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
   !!            4.0  !  2017-04  (G. Madec)  remove CPP ddm key & avm at t-point only
   !!----------------------------------------------------------------------

   !!----------------------------------------------------------------------
   !!   zdf_ric_init  : initialization, namelist read, & parameters control
   !!   zdf_ric       : update momentum and tracer Kz from the Richardson number
   !!   ric_rst       : read/write RIC restart in ocean restart file
   !!----------------------------------------------------------------------
   USE oce            ! ocean dynamics and tracers variables
   USE dom_oce        ! ocean space and time domain variables
   USE zdf_oce        ! vertical physics: variables
   USE phycst         ! physical constants
   USE sbc_oce,  ONLY :   taum
   !
   USE in_out_manager ! I/O manager
   USE iom            ! I/O manager library
   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)


   IMPLICIT NONE
   PRIVATE

   PUBLIC   zdf_ric         ! called by zdfphy.F90
   PUBLIC   ric_rst         ! called by zdfphy.F90
   PUBLIC   zdf_ric_init    ! called by nemogcm.F90

   !                        !!* 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

   !! * Substitutions
#  include "do_loop_substitute.h90"
#  include "domzgr_substitute.h90"
   !!----------------------------------------------------------------------
   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
   !! $Id: zdfric.F90 15277 2021-09-22 13:19:40Z clem $
   !! Software governed by the CeCILL license (see ./LICENSE)
   !!----------------------------------------------------------------------
CONTAINS

   SUBROUTINE zdf_ric_init
      !!----------------------------------------------------------------------
      !!                 ***  ROUTINE zdf_ric_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
      !!----------------------------------------------------------------------
      !
      READ  ( numnam_ref, namzdf_ric, IOSTAT = ios, ERR = 901)
901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_ric in reference namelist' )

      READ  ( numnam_cfg, namzdf_ric, IOSTAT = ios, ERR = 902 )
902   IF( ios >  0 ) CALL ctl_nam ( ios , 'namzdf_ric in configuration namelist' )
      IF(lwm) WRITE ( numond, namzdf_ric )
      !
      IF(lwp) THEN                   ! Control print
         WRITE(numout,*)
         WRITE(numout,*) 'zdf_ric_init : Ri depend vertical mixing scheme'
         WRITE(numout,*) '~~~~~~~~~~~~'
         WRITE(numout,*) '   Namelist namzdf_ric : set Kz=F(Ri) parameters'
         WRITE(numout,*) '      maximum vertical viscosity        rn_avmri  = ', rn_avmri
         WRITE(numout,*) '      coefficient                       rn_alp    = ', rn_alp
         WRITE(numout,*) '      exponent                          nn_ric    = ', nn_ric
         WRITE(numout,*) '      Ekman layer enhanced mixing       ln_mldw   = ', ln_mldw
         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
      ENDIF
      !
      CALL ric_rst( nit000, 'READ' )  !* read or initialize all required files
      !
   END SUBROUTINE zdf_ric_init


   SUBROUTINE zdf_ric( kt, Kmm, p_sh2, p_avm, p_avt )
      !!----------------------------------------------------------------------
      !!                 ***  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(uu(:,:,:,Kbb)) )+mj( dk(vv(:,:,:,Kbb)) ) ]
      !!                    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
      !!                and rn_alp, nn_ric are adjustable parameters.
      !!                Typical values : rn_alp=5. and nn_ric=2.
      !!
      !!      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.24) 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 zdfphy.F90
      !!
      !! ** Action  :   avm, avt  mixing coeff (inner domain values only)
      !!
      !! References : Pacanowski & Philander 1981, JPO, 1441-1451.
      !!              PFJ Lermusiaux 2001.
      !!----------------------------------------------------------------------
      INTEGER                         , INTENT(in   ) ::   kt             ! ocean time-step
      INTEGER                         , INTENT(in   ) ::   Kmm            ! ocean time level index
      REAL(wp), DIMENSION(A2D(0) ,jpk), INTENT(in   ) ::   p_sh2          ! shear production term
      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   p_avm          ! vertical eddy viscosity (w-points)
      REAL(wp), DIMENSION(A2D(0) ,jpk), INTENT(inout) ::   p_avt          ! vertical eddy diffusivity (w-points)
      !!
      INTEGER  ::   ji, jj, jk                  ! dummy loop indices
      REAL(wp) ::   zcfRi, zav, zustar, zhek    ! local scalars
      REAL(wp), DIMENSION(T2D(0)) ::   zh_ekm  ! 2D workspace
      !!----------------------------------------------------------------------
      !
      !                       !==  avm and avt = F(Richardson number)  ==!
      DO_3D( 0, 0, 0, 0, 2, jpkm1 )       ! coefficient = F(richardson number) (avm-weighted Ri)
         zcfRi = 1._wp / (  1._wp + rn_alp * MAX(  0._wp , avm(ji,jj,jk) * rn2(ji,jj,jk) / ( p_sh2(ji,jj,jk) + 1.e-20 ) )  )
         zav   = rn_avmri * zcfRi**nn_ric
         !                          ! avm and avt coefficients
         p_avm(ji,jj,jk) = MAX(  zav         , avmb(jk)  ) * wmask(ji,jj,jk)
         p_avt(ji,jj,jk) = MAX(  zav * zcfRi , avtb(jk)  ) * wmask(ji,jj,jk)
      END_3D
      !
!!gm BUG <<<<====  This param can't work at low latitude
!!gm               it provides there much to thick mixed layer ( summer 150m in GYRE configuration !!! )
      !
      IF( ln_mldw ) THEN      !==  set a minimum value in the Ekman layer  ==!
         !
         DO_2D( 0, 0, 0, 0 )
            zustar = SQRT( taum(ji,jj) * r1_rho0 )
            zhek   = rn_ekmfc * zustar / ( ABS( ff_t(ji,jj) ) + rsmall )   ! Ekman depth
            zh_ekm(ji,jj) = MAX(  rn_mldmin , MIN( zhek , rn_mldmax )  )   ! set allowed range
         END_2D
         DO_3D( 0, 0, 0, 0, 2, jpkm1 )   !* minimum mixing coeff. within the Ekman layer
            IF( gdept(ji,jj,jk,Kmm) < zh_ekm(ji,jj) ) THEN
               p_avm(ji,jj,jk) = MAX(  p_avm(ji,jj,jk), rn_wvmix  ) * wmask(ji,jj,jk)
               p_avt(ji,jj,jk) = MAX(  p_avt(ji,jj,jk), rn_wtmix  ) * wmask(ji,jj,jk)
            ENDIF
         END_3D
      ENDIF
      !
   END SUBROUTINE zdf_ric


   SUBROUTINE ric_rst( kt, cdrw )
      !!---------------------------------------------------------------------
      !!                   ***  ROUTINE ric_rst  ***
      !!
      !! ** Purpose :   Read or write TKE file (en) in restart file
      !!
      !! ** Method  :   use of IOM library
      !!                if the restart does not contain TKE, en is either
      !!                set to rn_emin or recomputed
      !!----------------------------------------------------------------------
      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
      !
      INTEGER ::   jit, jk    ! dummy loop indices
      INTEGER ::   id1, id2   ! local integers
      !!----------------------------------------------------------------------
      !
      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
         !                                   ! ---------------
         !           !* Read the restart file
         IF( ln_rstart ) THEN
            id1 = iom_varid( numror, 'avt_k', ldstop = .FALSE. )
            id2 = iom_varid( numror, 'avm_k', ldstop = .FALSE. )
            !
            IF( MIN( id1, id2 ) > 0 ) THEN         ! restart exists => read it
               CALL iom_get( numror, jpdom_auto, 'avt_k', avt_k )
               CALL iom_get( numror, jpdom_auto, 'avm_k', avm_k )
            ENDIF
         ENDIF
         !           !* otherwise Kz already set to the background value in zdf_phy_init
         !
      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
         !                                   ! -------------------
         IF(lwp) WRITE(numout,*) '---- ric-rst ----'
         CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k )
         CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k)
         !
      ENDIF
      !
   END SUBROUTINE ric_rst

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