Index: branches/2011/dev_INGV_2011/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90
===================================================================
--- branches/2011/dev_INGV_2011/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90	(revision 2991)
+++ branches/2011/dev_INGV_2011/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90	(revision 2992)
@@ -7,8 +7,9 @@
    !! History :  OPA  ! 1987-09  (P. Andrich)  Original code
    !!            4.0  ! 1991-11  (G. Madec)
-   !!            7.0  ! 1996-01  (G. Madec)  complet rewriting of multitasking suppression of common work arrays
-   !!            8.0  ! 1997-06 (G. Madec)  complete rewriting of zdfmix
+   !!            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
@@ -20,10 +21,11 @@
    !!   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 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 eosbn2, ONLY : nn_eos
 
    IMPLICIT NONE
@@ -39,4 +41,10 @@
    REAL(wp) ::   rn_avmri = 100.e-4_wp   ! maximum value of the vertical eddy viscosity
    REAL(wp) ::   rn_alp   =   5._wp      ! coefficient of the parameterization
+   REAL(wp) ::   rn_ekmfc =   0.7_wp     ! Ekman Factor Coeff
+   REAL(wp) ::   rn_mldmin=   1.0_wp     ! minimum mixed layer (ML) depth    
+   REAL(wp) ::   rn_mldmax=1000.0_wp     ! maximum mixed layer depth
+   REAL(wp) ::   rn_wtmix =  10.0_wp     ! Vertical eddy Diff. in the ML
+   REAL(wp) ::   rn_wvmix =  10.0_wp     ! Vertical eddy Visc. in the ML
+   LOGICAL  ::   ln_mldw  = .TRUE.       ! Use or not the MLD parameters
 
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tmric   !: coef. for the horizontal mean at t-point
@@ -67,14 +75,14 @@
       !!                    
       !! ** Purpose :   Compute the before eddy viscosity and diffusivity as
-      !!              a function of the local richardson number.
+      !!                a function of the local richardson number.
       !!
       !! ** Method  :   Local richardson number dependent formulation of the 
-      !!              vertical eddy viscosity and diffusivity coefficients. 
+      !!                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
+      !!                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 
@@ -84,19 +92,35 @@
       !!      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 wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
-      USE wrk_nemo, ONLY:   zwx => wrk_2d_1     ! 2D workspace
-      !!
-      INTEGER, INTENT( in ) ::   kt         ! ocean time-step indexocean time step
-      !!
-      INTEGER  ::   ji, jj, jk               ! dummy loop indices
-      REAL(wp) ::   zcoef, zdku, zdkv, zri, z05alp     ! temporary scalars
-      !!----------------------------------------------------------------------
-
-      IF( wrk_in_use(2, 1) ) THEN
+      USE wrk_nemo, ONLY:   zwx => wrk_2d_1         ! 2D workspace
+      USE wrk_nemo, ONLY:   ekm_dep => wrk_2d_2     ! 2D workspace
+      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
+      !!----------------------------------------------------------------------
+
+      IF( wrk_in_use(2, 1,2) ) THEN
          CALL ctl_stop('zdf_ric : requested workspace array unavailable')   ;   RETURN
       ENDIF
@@ -145,8 +169,60 @@
       !                                                ! ===============
       !
+      IF( ln_mldw ) THEN
+
+      !  Compute Ekman depth from wind stress forcing.
+      ! -------------------------------------------------------
+      zflageos = ( 0.5 + SIGN( 0.5, nn_eos - 1. ) ) * rau0
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zrhos          = rhop(ji,jj,1) + zflageos * ( 1. - tmask(ji,jj,1) )
+            zustar         = SQRT( taum(ji,jj) / ( zrhos +  rsmall ) )
+            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. )
       !
-      IF( wrk_not_released(2, 1) )   CALL ctl_stop('zdf_ric: failed to release workspace array')
+      IF( wrk_not_released(2, 1,2) )   CALL ctl_stop('zdf_ric: failed to release workspace array')
       !
    END SUBROUTINE zdf_ric
@@ -168,5 +244,6 @@
       INTEGER :: ji, jj, jk   ! dummy loop indices
       !!
-      NAMELIST/namzdf_ric/ rn_avmri, rn_alp, nn_ric
+      NAMELIST/namzdf_ric/ rn_avmri, rn_alp   , nn_ric  , rn_ekmfc,  &
+         &                rn_mldmin, rn_mldmax, rn_wtmix, rn_wvmix, ln_mldw
       !!----------------------------------------------------------------------
       !
@@ -179,7 +256,13 @@
          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,*) '      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
       !
