ecearth3/sources/nemo-3.6/NEMO/OPA_SRC/ZDF/zdfric.F90

MODULE zdfric
   !!======================================================================
   !!                       ***  MODULE  zdfric  ***
   !! Ocean physics:  vertical mixing coefficient compute from the local
   !!                 Richardson number dependent formulation
   !!======================================================================
   !! History :  OPA  ! 1987-09  (P. Andrich)  Original code
   !!            4.0  ! 1991-11  (G. Madec)
   !!            7.0  ! 1996-01  (G. Madec)  complete rewriting of multitasking suppression of common work arrays
   !!            8.0  ! 1997-06  (G. Madec)  complete rewriting of zdfmix
   !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module
   !!            3.3  ! 2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
   !!            3.3.1! 2011-09  (P. Oddo) Mixed layer depth parameterization
   !!----------------------------------------------------------------------
#if defined key_zdfric   ||   defined key_esopa
   !!----------------------------------------------------------------------
   !!   'key_zdfric'                                             Kz = f(Ri)
   !!----------------------------------------------------------------------
   !!   zdf_ric       : update momentum and tracer Kz from the Richardson
   !!                  number computation
   !!   zdf_ric_init  : initialization, namelist read, & parameters control
   !!----------------------------------------------------------------------
   USE oce            ! ocean dynamics and tracers variables
   USE dom_oce        ! ocean space and time domain variables
   USE zdf_oce        ! ocean vertical physics
   USE in_out_manager ! I/O manager
   USE lbclnk         ! ocean lateral boundary condition (or mpp link)
   USE lib_mpp        ! MPP library
   USE wrk_nemo       ! work arrays
   USE timing         ! Timing
   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  

   USE eosbn2, ONLY : nn_eos

   IMPLICIT NONE
   PRIVATE

   PUBLIC   zdf_ric         ! called by step.F90
   PUBLIC   zdf_ric_init    ! called by opa.F90

   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfric = .TRUE.   !: Richardson vertical mixing flag

   !                           !!* Namelist namzdf_ric : Richardson number dependent Kz *
   INTEGER  ::   nn_ric         ! coefficient of the parameterization
   REAL(wp) ::   rn_avmri       ! maximum value of the vertical eddy viscosity
   REAL(wp) ::   rn_alp         ! coefficient of the parameterization
   REAL(wp) ::   rn_ekmfc       ! Ekman Factor Coeff
   REAL(wp) ::   rn_mldmin      ! minimum mixed layer (ML) depth    
   REAL(wp) ::   rn_mldmax      ! maximum mixed layer depth
   REAL(wp) ::   rn_wtmix       ! Vertical eddy Diff. in the ML
   REAL(wp) ::   rn_wvmix       ! Vertical eddy Visc. in the ML
   LOGICAL  ::   ln_mldw        ! Use or not the MLD parameters

   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tmric   !: coef. for the horizontal mean at t-point

   !! * Substitutions
#  include "domzgr_substitute.h90"
   !!----------------------------------------------------------------------
   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
   !! $Id: zdfric.F90 4624 2014-04-28 12:09:03Z acc $
   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------
CONTAINS

   INTEGER FUNCTION zdf_ric_alloc()
      !!----------------------------------------------------------------------
      !!                 ***  FUNCTION zdf_ric_alloc  ***
      !!----------------------------------------------------------------------
      ALLOCATE( tmric(jpi,jpj,jpk)   , STAT= zdf_ric_alloc )
      !
      IF( lk_mpp             )   CALL mpp_sum ( zdf_ric_alloc )
      IF( zdf_ric_alloc /= 0 )   CALL ctl_warn('zdf_ric_alloc: failed to allocate arrays')
   END FUNCTION zdf_ric_alloc


   SUBROUTINE zdf_ric( kt )
      !!----------------------------------------------------------------------
      !!                 ***  ROUTINE zdfric  ***
      !!                    
      !! ** Purpose :   Compute the before eddy viscosity and diffusivity as
      !!                a function of the local richardson number.
      !!
      !! ** Method  :   Local richardson number dependent formulation of the 
      !!                vertical eddy viscosity and diffusivity coefficients. 
      !!                The eddy coefficients are given by:
      !!                    avm = avm0 + avmb
      !!                    avt = avm0 / (1 + rn_alp*ri)
      !!                with ri  = N^2 / dz(u)**2
      !!                         = e3w**2 * rn2/[ mi( dk(ub) )+mj( dk(vb) ) ]
      !!                    avm0= rn_avmri / (1 + rn_alp*ri)**nn_ric
      !!      Where ri is the before local Richardson number,
      !!            rn_avmri is the maximum value reaches by avm and avt 
      !!            avmb and avtb are the background (or minimum) values
      !!            and rn_alp, nn_ric are adjustable parameters.
      !!      Typical values used are : avm0=1.e-2 m2/s, avmb=1.e-6 m2/s
      !!      avtb=1.e-7 m2/s, rn_alp=5. and nn_ric=2.
      !!      a numerical threshold is impose on the vertical shear (1.e-20)
      !!      As second step compute Ekman depth from wind stress forcing
      !!      and apply namelist provided vertical coeff within this depth.
      !!      The Ekman depth is:
      !!              Ustar = SQRT(Taum/rho0)
      !!              ekd= rn_ekmfc * Ustar / f0
      !!      Large et al. (1994, eq.29) suggest rn_ekmfc=0.7; however, the derivation
      !!      of the above equation indicates the value is somewhat arbitrary; therefore
      !!      we allow the freedom to increase or decrease this value, if the
      !!      Ekman depth estimate appears too shallow or too deep, respectively.
      !!      Ekd is then limited by rn_mldmin and rn_mldmax provided in the
      !!      namelist
      !!        N.B. the mask are required for implicit scheme, and surface
      !!      and bottom value already set in zdfini.F90
      !!
      !! References : Pacanowski & Philander 1981, JPO, 1441-1451.
      !!              PFJ Lermusiaux 2001.
      !!----------------------------------------------------------------------
      USE phycst,   ONLY:   rsmall,rau0
      USE sbc_oce,  ONLY:   taum
      !!
      INTEGER, INTENT( in ) ::   kt                           ! ocean time-step
      !!
      INTEGER  ::   ji, jj, jk                                ! dummy loop indices
      REAL(wp) ::   zcoef, zdku, zdkv, zri, z05alp, zflageos  ! temporary scalars
      REAL(wp) ::   zrhos, zustar
      REAL(wp), POINTER, DIMENSION(:,:) ::   zwx, ekm_dep  
      !!----------------------------------------------------------------------
      !
      IF( nn_timing == 1 )  CALL timing_start('zdf_ric')
      !
      CALL wrk_alloc( jpi,jpj, zwx, ekm_dep )
      !                                                ! ===============
      DO jk = 2, jpkm1                                 ! Horizontal slab
         !                                             ! ===============
         ! Richardson number (put in zwx(ji,jj))
         ! -----------------
         DO jj = 2, jpjm1
            DO ji = 2, jpim1
               zcoef = 0.5 / fse3w(ji,jj,jk)
               !                                            ! shear of horizontal velocity
               zdku = zcoef * (  ub(ji-1,jj,jk-1) + ub(ji,jj,jk-1)   &
                  &             -ub(ji-1,jj,jk  ) - ub(ji,jj,jk  )  )
               zdkv = zcoef * (  vb(ji,jj-1,jk-1) + vb(ji,jj,jk-1)   &
                  &             -vb(ji,jj-1,jk  ) - vb(ji,jj,jk  )  )
               !                                            ! richardson number (minimum value set to zero)
               zri = rn2(ji,jj,jk) / ( zdku*zdku + zdkv*zdkv + 1.e-20 )
               zwx(ji,jj) = MAX( zri, 0.e0 )
            END DO
         END DO
         CALL lbc_lnk( zwx, 'W', 1. )                       ! Boundary condition   (sign unchanged)

         ! Vertical eddy viscosity and diffusivity coefficients
         ! -------------------------------------------------------
         z05alp = 0.5_wp * rn_alp
         DO jj = 1, jpjm1                                   ! Eddy viscosity coefficients (avm)
            DO ji = 1, jpim1
               avmu(ji,jj,jk) = umask(ji,jj,jk) * rn_avmri / ( 1. + z05alp*( zwx(ji+1,jj)+zwx(ji,jj) ) )**nn_ric
               avmv(ji,jj,jk) = vmask(ji,jj,jk) * rn_avmri / ( 1. + z05alp*( zwx(ji,jj+1)+zwx(ji,jj) ) )**nn_ric
            END DO
         END DO
         DO jj = 2, jpjm1                                   ! Eddy diffusivity coefficients (avt)
            DO ji = 2, jpim1
               avt(ji,jj,jk) = tmric(ji,jj,jk) / ( 1._wp + rn_alp * zwx(ji,jj) )           &
                  &                            * (  avmu(ji,jj,jk) + avmu(ji-1,jj,jk)      &
                  &                               + avmv(ji,jj,jk) + avmv(ji,jj-1,jk)  )   &
                  &          + avtb(jk) * tmask(ji,jj,jk)
            END DO
         END DO
         DO jj = 2, jpjm1                                   ! Add the background coefficient on eddy viscosity
            DO ji = 2, jpim1
               avmu(ji,jj,jk) = avmu(ji,jj,jk) + avmb(jk) * umask(ji,jj,jk)
               avmv(ji,jj,jk) = avmv(ji,jj,jk) + avmb(jk) * vmask(ji,jj,jk)
            END DO
         END DO
         !                                             ! ===============
      END DO                                           !   End of slab
      !                                                ! ===============
      !
      IF( ln_mldw ) THEN

      !  Compute Ekman depth from wind stress forcing.
      ! -------------------------------------------------------
      !SF  zflageos = ( 0.5 + SIGN( 0.5, nn_eos - 1. ) ) * rau0
      !SF  DO jj = 1, jpj
      !SF     DO ji = 1, jpi
      !SF        zrhos          = rhop(ji,jj,1) + zflageos * ( 1. - tmask(ji,jj,1) )
      !SF        zustar         = SQRT( taum(ji,jj) / ( zrhos +  rsmall ) )
      !SF        ekm_dep(ji,jj) = rn_ekmfc * zustar / ( ABS( ff(ji,jj) ) + rsmall )
      !SF        ekm_dep(ji,jj) = MAX(ekm_dep(ji,jj),rn_mldmin) ! Minimun allowed
      !SF        ekm_dep(ji,jj) = MIN(ekm_dep(ji,jj),rn_mldmax) ! Maximum allowed
      !SF     END DO
      !SF  END DO

      DO jj = 1, jpj
         DO ji = 1, jpi
            zustar         = SQRT( taum(ji,jj) * r1_rau0 )
            ekm_dep(ji,jj) = rn_ekmfc * zustar / ( ABS( ff(ji,jj) ) + rsmall )
            ekm_dep(ji,jj) = MAX(ekm_dep(ji,jj),rn_mldmin) ! Minimun allowed
            ekm_dep(ji,jj) = MIN(ekm_dep(ji,jj),rn_mldmax) ! Maximum allowed
         END DO
      END DO

      ! In the first model level vertical diff/visc coeff.s 
      ! are always equal to the namelist values rn_wtmix/rn_wvmix
      ! -------------------------------------------------------
      DO jj = 1, jpj
         DO ji = 1, jpi
            avmv(ji,jj,1) = MAX( avmv(ji,jj,1), rn_wvmix )
            avmu(ji,jj,1) = MAX( avmu(ji,jj,1), rn_wvmix )
            avt( ji,jj,1) = MAX(  avt(ji,jj,1), rn_wtmix )
         END DO
      END DO

      !  Force the vertical mixing coef within the Ekman depth
      ! -------------------------------------------------------
      DO jk = 2, jpkm1
         DO jj = 1, jpj
            DO ji = 1, jpi
               IF( fsdept(ji,jj,jk) < ekm_dep(ji,jj) ) THEN
                  avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk), rn_wvmix )
                  avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk), rn_wvmix )
                  avt( ji,jj,jk) = MAX(  avt(ji,jj,jk), rn_wtmix )
               ENDIF
            END DO
         END DO
      END DO

      DO jk = 1, jpkm1                
         DO jj = 1, jpj
            DO ji = 1, jpi
               avmv(ji,jj,jk) = avmv(ji,jj,jk) * vmask(ji,jj,jk)
               avmu(ji,jj,jk) = avmu(ji,jj,jk) * umask(ji,jj,jk)
               avt( ji,jj,jk) = avt( ji,jj,jk) * tmask(ji,jj,jk)
            END DO
         END DO
      END DO

     ENDIF

      CALL lbc_lnk( avt , 'W', 1. )                         ! Boundary conditions   (unchanged sign)
      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )
      !
      CALL wrk_dealloc( jpi,jpj, zwx, ekm_dep )
      !
      IF( nn_timing == 1 )  CALL timing_stop('zdf_ric')
      !
   END SUBROUTINE zdf_ric


   SUBROUTINE zdf_ric_init
      !!----------------------------------------------------------------------
      !!                 ***  ROUTINE zdfbfr_init  ***
      !!                    
      !! ** Purpose :   Initialization of the vertical eddy diffusivity and
      !!      viscosity coef. for the Richardson number dependent formulation.
      !!
      !! ** Method  :   Read the namzdf_ric namelist and check the parameter values
      !!
      !! ** input   :   Namelist namzdf_ric
      !!
      !! ** Action  :   increase by 1 the nstop flag is setting problem encounter
      !!----------------------------------------------------------------------
      INTEGER :: ji, jj, jk   ! dummy loop indices
      INTEGER ::   ios        ! Local integer output status for namelist read
      !!
      NAMELIST/namzdf_ric/ rn_avmri, rn_alp   , nn_ric  , rn_ekmfc,  &
         &                rn_mldmin, rn_mldmax, rn_wtmix, rn_wvmix, ln_mldw
      !!----------------------------------------------------------------------
      !
      REWIND( numnam_ref )              ! Namelist namzdf_ric in reference namelist : Vertical diffusion Kz depends on Richardson number
      READ  ( numnam_ref, namzdf_ric, IOSTAT = ios, ERR = 901)
901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_ric in reference namelist', lwp )

      REWIND( numnam_cfg )              ! Namelist namzdf_ric in configuration namelist : Vertical diffusion Kz depends on Richardson number
      READ  ( numnam_cfg, namzdf_ric, IOSTAT = ios, ERR = 902 )
902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_ric in configuration namelist', lwp )
      IF(lwm) WRITE ( numond, namzdf_ric )
      !
      IF(lwp) THEN                   ! Control print
         WRITE(numout,*)
         WRITE(numout,*) 'zdf_ric : Ri depend vertical mixing scheme'
         WRITE(numout,*) '~~~~~~~'
         WRITE(numout,*) '   Namelist namzdf_ric : set Kz(Ri) parameters'
         WRITE(numout,*) '      maximum vertical viscosity     rn_avmri  = ', rn_avmri
         WRITE(numout,*) '      coefficient                    rn_alp    = ', rn_alp
         WRITE(numout,*) '      coefficient                    nn_ric    = ', nn_ric
         WRITE(numout,*) '      Ekman Factor Coeff             rn_ekmfc  = ', rn_ekmfc
         WRITE(numout,*) '      minimum mixed layer depth      rn_mldmin = ', rn_mldmin
         WRITE(numout,*) '      maximum mixed layer depth      rn_mldmax = ', rn_mldmax
         WRITE(numout,*) '      Vertical eddy Diff. in the ML  rn_wtmix  = ', rn_wtmix
         WRITE(numout,*) '      Vertical eddy Visc. in the ML  rn_wvmix  = ', rn_wvmix
         WRITE(numout,*) '      Use the MLD parameterization   ln_mldw   = ', ln_mldw
      ENDIF
      !
      !                              ! allocate zdfric arrays
      IF( zdf_ric_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_ric_init : unable to allocate arrays' )
      !
      DO jk = 1, jpk                 ! weighting mean array tmric for 4 T-points
         DO jj = 2, jpj              ! which accounts for coastal boundary conditions            
            DO ji = 2, jpi
               tmric(ji,jj,jk) =  tmask(ji,jj,jk)                                  &
                  &            / MAX( 1.,  umask(ji-1,jj  ,jk) + umask(ji,jj,jk)   &
                  &                      + vmask(ji  ,jj-1,jk) + vmask(ji,jj,jk)  )
            END DO
         END DO
      END DO
      tmric(:,1,:) = 0._wp
      !
      DO jk = 1, jpk                 ! Initialization of vertical eddy coef. to the background value
         avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
         avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
         avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
      END DO
      !
   END SUBROUTINE zdf_ric_init

#else
   !!----------------------------------------------------------------------
   !!   Dummy module :              NO Richardson dependent vertical mixing
   !!----------------------------------------------------------------------
   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfric = .FALSE.   !: Richardson mixing flag
CONTAINS
   SUBROUTINE zdf_ric_init         ! Dummy routine
   END SUBROUTINE zdf_ric_init
   SUBROUTINE zdf_ric( kt )        ! Dummy routine
      WRITE(*,*) 'zdf_ric: You should not have seen this print! error?', kt
   END SUBROUTINE zdf_ric
#endif

   !!======================================================================
END MODULE zdfric