Index: /trunk/NEMO/OPA_SRC/ZDF/zdftke.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 1049)
+++ /trunk/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 1050)
@@ -18,4 +18,9 @@
    !!             9.0  !  02-08  (G. Madec)  autotasking optimization
    !!             9.0  !  06-07  (S. Masson)  distributed restart using iom
+   !!              -   !  2005-07  (C. Ethe,  G.Madec) : update TKE physics:
+   !!                              - tke penetration (wind steering)
+   !!                              - suface condition for tke & mixing length
+   !!                              - Langmuir cells
+   !!             3.0  !  2007-11  (G. Madec)  avtb_2d, nn_avtb_2d 
    !!----------------------------------------------------------------------
 #if defined key_zdftke   ||   defined key_esopa
@@ -33,4 +38,5 @@
    USE sbc_oce         ! surface boundary condition: ocean
    USE phycst          ! physical constants
+   USE zdfmxl          ! mixed layer
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE prtctl          ! Print control
@@ -67,4 +73,15 @@
       &                  emin0 = 1.e-4_wp     ,  &  !: surface minimum value of tke (m2/s2)
       &                  ri_c  = 2._wp / 9._wp      !: critic Richardson number
+   !                                      !!! ** Namelist (namtke) **
+   INTEGER  ::   nn_avtb_2d = 1            !: = 0/1 horizontal shape for avtb
+   INTEGER  ::   nn_etau    = 0            !: = 0/1/2 tke depth penetration
+   INTEGER  ::   nn_htau    = 0            !: = 0/1/2 type of tke profile of penetration
+   REAL(wp) ::   rn_lmin    = 0.4_wp       !: surface min value of mixing turbulent length scale
+   REAL(wp) ::   rn_efr     = 1.0_wp       !: fraction of TKE surface value which penetrates in the ocean
+   LOGICAL  ::   ln_lsfc    = .FALSE.      !: mixing length scale surface value as function of wind stress or not
+   LOGICAL  ::   ln_lc      = .FALSE.      !: Langmuir cells (LC) as a source term of TKE or not
+   REAL(wp) ::   rn_lc      = 0.15_wp      !: coef to compute vertical velocity of LC
+
+   REAL(wp), DIMENSION (jpi,jpj)  :: avtb_2d ! set in tke_init, for other modif than ice
 
 # if defined key_c1d
@@ -148,4 +165,5 @@
          &          zmxlm => ta,  &  ! use ta as workspace
          &          zmxld => sa      ! use sa as workspace
+      USE oce     , ztkelc => va     ! use va as workspace
       !
       INTEGER, INTENT(in) ::   kt ! ocean time step
@@ -160,4 +178,9 @@
          &          zsqen, zesh2,            &  !
          &          zemxl, zemlm, zemlp
+      INTEGER , DIMENSION(jpi,jpj)     ::   imlc    ! 2D workspace
+      REAL(wp) ::   zraug, zus, zwlc, zind     ! temporary scalar
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zhtau   ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zhlc    ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zpelc   ! 3D workspace
       !!--------------------------------------------------------------------
 
@@ -165,5 +188,5 @@
 
       !                                            ! Local constant initialization
-      zmlmin = 1.e-8
+      zmlmin = 0.4
       zbbrau =  .5 * ebb / rau0
       zfact1 = -.5 * rdt * efave
@@ -177,6 +200,21 @@
       ! Buoyancy length scale: l=sqrt(2*e/n**2)
       ! ---------------------
-      zmxlm(:,:, 1 ) = zmlmin   ! surface set to the minimum value
-      zmxlm(:,:,jpk) = zmlmin   ! bottom  set to the minimum value
+      IF( ln_lsfc ) THEN      ! lsurf is function of wind stress : l=2.e5*sqrt(utau^2 + vtau^2)/(rau0*g)
+         zmxlm(:,:,1) = 0.e0
+         zraug = 0.5 * 2.e5 / ( rau0 * grav )
+         DO jj = 2, jpjm1
+!CDIR NOVERRCHK
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
+               zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
+               zmxlm(ji,jj,1  ) = zraug * SQRT( ztx2 * ztx2 + zty2 * zty2 )
+               ! set the surface minimum value to rn_lmin
+               zmxlm(ji,jj,1  ) = MAX( zmxlm(ji,jj,1) , rn_lmin )
+            END DO
+         END DO
+      ELSE                    ! surface set to the minimum value
+         zmxlm(:,:,1) = rn_lmin  
+      ENDIF
+      zmxlm(:,:,jpk) = rn_lmin   ! bottom  set to the minimum value
 !CDIR NOVERRCHK
       DO jk = 2, jpkm1
@@ -193,6 +231,6 @@
       ! Physical limits for the mixing length
       ! -------------------------------------
-      zmxld(:,:, 1 ) = zmlmin   ! surface set to the minimum value
-      zmxld(:,:,jpk) = zmlmin   ! bottom  set to the minimum value
+      zmxld(:,:, 1 ) = zmxlm(:,:,1)   ! surface set to the minimum value
+      zmxld(:,:,jpk) = rn_lmin           ! bottom  set to the minimum value
 
       SELECT CASE ( nmxl )
@@ -277,5 +315,60 @@
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! II  Tubulent kinetic energy time stepping
+      ! II  TKE Langmuir circulation source term
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( ln_lc ) THEN
+         !
+         ! Computation of total energy produce by LC : cumulative sum over jk
+         zpelc(:,:,1) =  MAX( rn2(:,:,1), 0. ) * fsdepw(:,:,1) * fse3w(:,:,1)
+         DO jk = 2, jpk
+            zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2(:,:,jk), 0. ) * fsdepw(:,:,jk) * fse3w(:,:,jk)
+         END DO
+         !
+         ! Computation of finite Langmuir Circulation depth
+         ! Initialization to the number of w ocean point mbathy
+         imlc(:,:) = mbathy(:,:)
+         DO jk = jpkm1, 2, -1
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ! Last w-level at which zpelc>=0.5*us*us with us=0.016*wind(starting from jpk-1)
+                  zus  = 0.000128 * wndm(ji,jj) * wndm(ji,jj)
+                  IF( zpelc(ji,jj,jk) > zus )   imlc(ji,jj) = jk
+               END DO
+            END DO
+         END DO
+         !
+         ! finite LC depth
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zhlc(ji,jj) = fsdepw(ji,jj,imlc(ji,jj))
+            END DO
+         END DO
+         !
+# if defined key_c1d
+         hlc(:,:) = zhlc(:,:) * tmask(:,:,1)      ! save finite Langmuir Circulation  depth
+# endif
+         !
+         ! TKE Langmuir circulation source term
+!CDIR NOVERRCHK
+         DO jk = 2, jpkm1
+!CDIR NOVERRCHK
+            DO jj = 2, jpjm1
+!CDIR NOVERRCHK
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  ! Stokes drift
+                  zus  = 0.016 * wndm(ji,jj)
+                  ! computation of vertical velocity due to LC
+                  zind = 0.5 - SIGN( 0.5, fsdepw(ji,jj,jk) - zhlc(ji,jj) )
+                  zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) )
+                  ! TKE Langmuir circulation source term
+                  ztkelc(ji,jj,jk) = ( zwlc * zwlc * zwlc ) / zhlc(ji,jj)
+               END DO
+            END DO
+         END DO
+         !
+      ENDIF
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! III  Tubulent kinetic energy time stepping
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
@@ -290,5 +383,5 @@
                zsqen = SQRT( en(ji,jj,jk) )
                zav   = ediff * zmxlm(ji,jj,jk) * zsqen
-               avt  (ji,jj,jk) = MAX( zav, avtb(jk) ) * tmask(ji,jj,jk)
+               avt  (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
                zmxlm(ji,jj,jk) = MAX( zav, avmb(jk) ) * tmask(ji,jj,jk)
                zmxld(ji,jj,jk) = zsqen / zmxld(ji,jj,jk)
@@ -309,5 +402,7 @@
             zesurf = zbbrau * SQRT( ztx2 * ztx2 + zty2 * zty2 )
             en (ji,jj,1) = MAX( zesurf, emin0 ) * tmask(ji,jj,1)
-            zmxlm(ji,jj,1  ) = avmb(1) * tmask(ji,jj,1)
+            zav  =  ediff * zmxlm(ji,jj,1) * SQRT( en(ji,jj,1) )* tmask(ji,jj,1)
+            zmxlm(ji,jj,1) = MAX( zav, avmb(1) ) * tmask(ji,jj,1)
+            avt  (ji,jj,1) = MAX( zav, avtb(1) * avtb_2d(ji,jj) ) * tmask(ji,jj,1)
             zmxlm(ji,jj,jpk) = 0.e0
          END DO
@@ -347,6 +442,12 @@
                   zwd(ji,jj,jk) = 1. - avmv(ji,jj,jk) - avmu(ji,jj,jk) + zfact2 * zmxld(ji,jj,jk)
                   ! right hand side in en 
-                  en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld(ji,jj,jk) * en   (ji,jj,jk)   &
-                  &                           +   rdt  * zmxlm(ji,jj,jk) * zesh2
+                  IF( .NOT. ln_lc ) THEN
+                     en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld (ji,jj,jk) * en(ji,jj,jk)   &
+                        &                        +   rdt  * zmxlm (ji,jj,jk) * zesh2
+                  ELSE
+                     en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld (ji,jj,jk) * en(ji,jj,jk)   &
+                        &                        +   rdt  * zmxlm (ji,jj,jk) * zesh2          &
+                        &                        +   rdt  * ztkelc(ji,jj,jk)
+                  ENDIF
                END DO
             END DO
@@ -390,7 +491,13 @@
                   zwd(ji,jj,jk) = 1. - avmv(ji,jj,jk) - avmu(ji,jj,jk) + zfact2 * zmxld(ji,jj,jk)
                   ! right hand side in en 
-                  en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld(ji,jj,jk) * en   (ji,jj,jk)   &
-                  &                           +   rdt  * zmxlm(ji,jj,jk) * zesh2
-                  ! save masked Prandlt number in zmxlm array
+                  IF( .NOT. ln_lc ) THEN
+                     en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld (ji,jj,jk) * en   (ji,jj,jk)   &
+                        &                        +   rdt  * zmxlm (ji,jj,jk) * zesh2
+                  ELSE
+                     en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld (ji,jj,jk) * en   (ji,jj,jk)   &
+                        &                        +   rdt  * zmxlm (ji,jj,jk) * zesh2             &
+                        &                        +   rdt  * ztkelc(ji,jj,jk)
+                  ENDIF
+                  ! save masked Prandlt number in zmxld array
                   zmxld(ji,jj,jk) = zpdl * tmask(ji,jj,jk)
                END DO
@@ -454,4 +561,54 @@
          END DO
       END DO
+     
+      ! Modify the value of TKE by an exponential assumption
+      ! en(ji,jj,jk)=en(ji,jj,jk)+en(ji,jj,1)*EXP(-fsdepw(ji,jj,jk)/ zhtau(ji,jj) )
+
+      SELECT CASE( nn_htau )      ! Choice of H value
+      !
+      CASE( 0 )
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zhtau(ji,jj) = 10.
+            END DO
+         END DO
+         !
+      CASE( 1 )
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zhtau(ji,jj) = MAX( 5., MIN( 40., 60.*ABS( SIN( rpi/180. * gphit(ji,jj) ) ) )   )
+            END DO
+         END DO
+         !
+      CASE( 2 )
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zhtau(ji,jj) =  MAX( 5.,6./4.* MIN( 40., 60.*ABS( SIN( rpi/180. * gphit(ji,jj) ) ) )   )
+            END DO
+         END DO
+         !
+      END SELECT
+
+      IF( nn_etau == 1 ) THEN
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  en(ji,jj,jk) = en(ji,jj,jk)   &
+                     &         + rn_efr * en(ji,jj,1)*EXP( -fsdepw(ji,jj,jk) / zhtau(ji,jj) ) &
+                     &                  * ( 1.e0 - fr_i(ji,jj) )  * tmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+      ELSEIF( nn_etau == 2 ) THEN     !  only at the base of the mixed layer
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               jk = nmln(ji,jj)
+               en(ji,jj,jk) = en(ji,jj,jk)    &
+                  &         + rn_efr * en(ji,jj,1)*EXP( -fsdepw(ji,jj,jk) / zhtau(ji,jj) ) &
+                  &                  * ( 1.e0 - fr_i(ji,jj) )  * tmask(ji,jj,jk)
+            END DO
+         END DO
+      ENDIF
+
 
       ! Lateral boundary conditions on ( avt, en )   (sign unchanged)
@@ -459,5 +616,5 @@
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! III.  Before vertical eddy vicosity and diffusivity coefficients
+      ! IV.  Before vertical eddy vicosity and diffusivity coefficients
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
@@ -563,5 +720,5 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   zpdl = zmxld(ji,jj,jk)
-                  avt(ji,jj,jk) = MAX( zpdl * avt(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk)
+                  avt(ji,jj,jk) = MAX( zpdl * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
                END DO
             END DO
@@ -621,7 +778,9 @@
       INTEGER ::           jk   ! dummy loop indices
 # endif
-
+      !!
       NAMELIST/namtke/ ln_rstke, ediff, ediss, ebb, efave, emin, emin0,   &
-         &             ri_c, nitke, nmxl, npdl, nave, navb
+                       ri_c, nitke, nmxl, npdl, nave, navb,    &
+                       ln_lsfc, rn_lmin, nn_avtb_2d, nn_etau, nn_htau, rn_efr, &
+                       ln_lc , rn_lc 
       !!----------------------------------------------------------------------
 
@@ -659,4 +818,12 @@
          WRITE(numout,*) '                and its associated coeff.   eboost = ', eboost
          WRITE(numout,*) '             constant background or profile navb   = ', navb
+         WRITE(numout,*) '             flag for compu.of bls as fun. of wind     ln_lsfc    = ', ln_lsfc
+         WRITE(numout,*) '             buoyancy lenght scale surface min value   rn_lmin    = ', rn_lmin
+         WRITE(numout,*) '             horizontal variation for avtb             nn_avtb_2d = ', nn_avtb_2d
+         WRITE(numout,*) '             test param. to add tke induced by wind    nn_etau    = ', nn_etau
+         WRITE(numout,*) '             flag for computation of exp. tke profile  nn_htau    = ', nn_htau
+         WRITE(numout,*) '             % of emin which pene. the thermocline     rn_efr     = ', rn_efr
+         WRITE(numout,*) '             flag to take into acc.  Langmuir circ.    ln_lc      = ', ln_lc
+         WRITE(numout,*) '             coef to compute verticla velocity of LC   rn_lc      = ', rn_lc
          WRITE(numout,*)
       ENDIF
@@ -665,4 +832,8 @@
       IF( nmxl < 0 .OR. nmxl > 3 ) CALL ctl_stop( '          bad flag: nmxl is < 0 or > 3 ' )
       IF( npdl < 0 .OR. npdl > 1 ) CALL ctl_stop( '          bad flag: npdl is < 0 or > 1 ' )
+      IF( nn_htau < 0    .OR. nn_htau > 2   )   CALL ctl_stop( 'bad flag: nn_htau is < 0 or > 3 ' )
+      IF( rn_lc   < 0.15 .OR. rn_lc   > 0.2 )   CALL ctl_stop( 'bad value: rn_lc must be > 0.15 and < 0.2 ' )
+
+      IF( nn_etau == 2 )   CALL zdf_mxl( nit000 )      ! Initialization of nmln 
 
       ! Horizontal average : initialization of weighting arrays 
@@ -743,19 +914,31 @@
       ! Background eddy viscosity and diffusivity profil
       ! ------------------------------------------------
-      IF( navb == 0 ) THEN
-         !   Define avmb, avtb from namelist parameter
+      IF( navb == 0 ) THEN      ! Define avmb, avtb from namelist parameter
          avmb(:) = avm0
          avtb(:) = avt0
-      ELSE
-         !   Background profile of avt (fit a theoretical/observational profile (Krauss 1990) 
+      ELSE                      ! Background profile of avt (fit a theoretical/observational profile (Krauss 1990) 
          avmb(:) = avm0
 !!bug  this is not valide neither in scoord
          IF(ln_sco .AND. lwp) WRITE(numout,cform_war)
-         IF(ln_sco .AND. lwp) WRITE(numout,*) '          avtb profile nort valid in sco'
+         IF(ln_sco .AND. lwp) WRITE(numout,*) '          avtb profile not valid in sco'
 
          avtb(:) = avt0 + ( 3.0e-4 - 2 * avt0 ) * 1.0e-4 * gdepw_0(:)   ! m2/s
       ENDIF
 
+      !                         ! 2D shape of the avtb
+      avtb_2d(:,:) = 1.e0             ! uniform 
+      !
+      IF( nn_avtb_2d == 1 ) THEN      ! decrease avtb in the equatorial band
+           !  -15S -5S : linear decrease from avt0 to avt0/10.
+           !  -5S  +5N : cst value avt0/10.
+           !   5N  15N : linear increase from avt0/10, to avt0
+           WHERE(-15. <= gphit .AND. gphit < -5 )   avtb_2d = (1.  - 0.09 * (gphit + 15.))
+           WHERE( -5. <= gphit .AND. gphit <  5 )   avtb_2d =  0.1
+           WHERE(  5. <= gphit .AND. gphit < 15 )   avtb_2d = (0.1 + 0.09 * (gphit -  5.))
+      ENDIF
+
       !   Increase the background in the surface layers
+!!gm  the increase is only required when using cen2 advection scheme 
+!!gm  for the equatorial upwelling.
       avmb(1) = 10.  * avmb(1)      ;      avtb(1) = 10.  * avtb(1)
       avmb(2) = 10.  * avmb(2)      ;      avtb(2) = 10.  * avtb(2)
