Index: branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
===================================================================
--- branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 5260)
+++ branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 5989)
@@ -28,5 +28,5 @@
    !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability
    !!----------------------------------------------------------------------
-#if defined key_zdftke   ||   defined key_esopa
+#if defined key_zdftke
    !!----------------------------------------------------------------------
    !!   'key_zdftke'                                   TKE vertical physics
@@ -53,4 +53,8 @@
    USE timing         ! Timing
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+#if defined key_agrif
+   USE agrif_opa_interp
+   USE agrif_opa_update
+#endif
 
    IMPLICIT NONE
@@ -85,9 +89,7 @@
    REAL(wp) ::   rhftau_scl = 1.0_wp       ! scale factor applied to HF part of taum  (nn_etau=3)
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   en             !: now turbulent kinetic energy   [m2/s2]
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htau           ! depth of tke penetration (nn_htau)
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dissl          ! now mixing lenght of dissipation
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avt_k , avm_k  ! not enhanced Kz
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avmu_k, avmv_k ! not enhanced Kz
+   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   apdlr          ! now mixing lenght of dissipation
 #if defined key_c1d
    !                                                                        !!** 1D cfg only  **   ('key_c1d')
@@ -100,5 +102,5 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
+   !! NEMO/OPA 3.7 , NEMO Consortium (2015)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -115,7 +117,6 @@
          &      e_pdl(jpi,jpj,jpk) , e_ric(jpi,jpj,jpk) ,                          &
 #endif
-         &      en    (jpi,jpj,jpk) , htau  (jpi,jpj)    , dissl(jpi,jpj,jpk) ,     & 
-         &      avt_k (jpi,jpj,jpk) , avm_k (jpi,jpj,jpk),                          &
-         &      avmu_k(jpi,jpj,jpk) , avmv_k(jpi,jpj,jpk), STAT= zdf_tke_alloc      )
+         &      htau  (jpi,jpj)    , dissl(jpi,jpj,jpk) ,     & 
+         &      apdlr(jpi,jpj,jpk) ,                                           STAT= zdf_tke_alloc      )
          !
       IF( lk_mpp             )   CALL mpp_sum ( zdf_tke_alloc )
@@ -173,4 +174,9 @@
       !!----------------------------------------------------------------------
       !
+#if defined key_agrif 
+      ! interpolation parent grid => child grid for avm_k ( ex : at west border: update column 1 and 2)
+      IF( .NOT.Agrif_Root() )   CALL Agrif_Tke
+#endif
+      !
       IF( kt /= nit000 ) THEN   ! restore before value to compute tke
          avt (:,:,:) = avt_k (:,:,:) 
@@ -189,5 +195,10 @@
       avmv_k(:,:,:) = avmv(:,:,:) 
       !
-   END SUBROUTINE zdf_tke
+#if defined key_agrif
+      ! Update child grid f => parent grid 
+      IF( .NOT.Agrif_Root() )   CALL Agrif_Update_Tke( kt )      ! children only
+#endif      
+     ! 
+  END SUBROUTINE zdf_tke
 
 
@@ -221,14 +232,15 @@
       REAL(wp) ::   zzd_up, zzd_lw                  !    -         -
 !!bfr      REAL(wp) ::   zebot                           !    -         -
-      INTEGER , POINTER, DIMENSION(:,:  ) :: imlc
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: zhlc
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zpelc, zdiag, zd_up, zd_lw
+      INTEGER , POINTER, DIMENSION(:,:  ) ::   imlc
+      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zhlc
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zpelc, zdiag, zd_up, zd_lw, z3du, z3dv
+      REAL(wp)                            ::   zri  !   local Richardson number
       !!--------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('tke_tke')
       !
-      CALL wrk_alloc( jpi,jpj, imlc )    ! integer
-      CALL wrk_alloc( jpi,jpj, zhlc ) 
-      CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 
+      CALL wrk_alloc( jpi,jpj,       imlc )    ! integer
+      CALL wrk_alloc( jpi,jpj,       zhlc ) 
+      CALL wrk_alloc( jpi,jpj,jpk,   zpelc, zdiag, zd_up, zd_lw, z3du, z3dv ) 
       !
       zbbrau = rn_ebb / rau0       ! Local constant initialisation
@@ -244,5 +256,5 @@
          DO jj = 2, jpjm1            ! en(mikt(ji,jj))   = rn_emin
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               en(ji,jj,mikt(ji,jj))=rn_emin * tmask(ji,jj,1)
+               en(ji,jj,mikt(ji,jj)) = rn_emin * tmask(ji,jj,1)
             END DO
          END DO
@@ -265,7 +277,5 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !                     en(bot)   = (rn_ebb0/rau0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin)
-!CDIR NOVERRCHK
 !!    DO jj = 2, jpjm1
-!CDIR NOVERRCHK
 !!       DO ji = fs_2, fs_jpim1   ! vector opt.
 !!          ztx2 = bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) + &
@@ -306,9 +316,6 @@
          END DO
          zcof = 0.016 / SQRT( zrhoa * zcdrag )
-!CDIR NOVERRCHK
          DO jk = 2, jpkm1         !* TKE Langmuir circulation source term added to en
-!CDIR NOVERRCHK
-            DO jj = 2, jpjm1
-!CDIR NOVERRCHK
+            DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   zus  = zcof * SQRT( taum(ji,jj) )           ! Stokes drift
@@ -332,18 +339,37 @@
       !
       DO jk = 2, jpkm1           !* Shear production at uw- and vw-points (energy conserving form)
-         DO jj = 1, jpj                 ! here avmu, avmv used as workspace
-            DO ji = 1, jpi
-               avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) )   &
-                  &                            * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) )   & 
-                  &                            / (  fse3uw_n(ji,jj,jk)               &
-                  &                              *  fse3uw_b(ji,jj,jk)  )
-               avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) )   &
-                  &                            * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) )   &
-                  &                            / (  fse3vw_n(ji,jj,jk)               &
-                  &                              *  fse3vw_b(ji,jj,jk)  )
-            END DO
-         END DO
-      END DO
-      !
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               z3du(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk  ) + avm(ji+1,jj,jk) )   &
+                  &                 * (  un(ji,jj,jk-1) -  un(ji  ,jj,jk) )   &
+                  &                 * (  ub(ji,jj,jk-1) -  ub(ji  ,jj,jk) ) * wumask(ji,jj,jk) &
+                  &                 / (  fse3uw_n(ji,jj,jk) * fse3uw_b(ji,jj,jk) )
+               z3dv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk  ) + avm(ji,jj+1,jk) )   &
+                  &                 * (  vn(ji,jj,jk-1) -  vn(ji,jj  ,jk) )   &
+                  &                 * (  vb(ji,jj,jk-1) -  vb(ji,jj  ,jk) ) * wvmask(ji,jj,jk) &
+                  &                 / (  fse3vw_n(ji,jj,jk) * fse3vw_b(ji,jj,jk) )
+            END DO
+         END DO
+      END DO
+      !
+      IF( nn_pdl == 1 ) THEN      !* Prandtl number case: compute apdlr
+         ! Note that zesh2 is also computed in the next loop.
+         ! We decided to compute it twice to keep code readability and avoid an IF case in the DO loops
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  !                                          ! shear prod. at w-point weightened by mask
+                  zesh2  =  ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
+                     &    + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )    
+                  !                                          ! local Richardson number
+                  zri   = MAX( rn2b(ji,jj,jk), 0._wp ) * avm(ji,jj,jk) / ( zesh2 + rn_bshear )
+                  apdlr(ji,jj,jk) = MAX(  0.1_wp,  ri_cri / MAX( ri_cri , zri )  )
+                  
+               END DO
+            END DO
+         END DO
+         !
+      ENDIF
+      !         
       DO jk = 2, jpkm1           !* Matrix and right hand side in en
          DO jj = 2, jpjm1
@@ -354,8 +380,8 @@
                zzd_lw = zcof * ( avm  (ji,jj,jk  ) + avm  (ji,jj,jk-1) )   &  ! lower diagonal
                   &          / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) )
-                  !                                                           ! shear prod. at w-point weightened by mask
-               zesh2  =  ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
-                  &    + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )    
-                  !
+               !                                   ! shear prod. at w-point weightened by mask
+               zesh2  =  ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
+                  &    + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )    
+               !
                zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw)
                zd_lw(ji,jj,jk) = zzd_lw
@@ -377,7 +403,5 @@
          END DO
       END DO
-      !
-      ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
-      DO jj = 2, jpjm1
+      DO jj = 2, jpjm1                             ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
          DO ji = fs_2, fs_jpim1   ! vector opt.
             zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
@@ -391,7 +415,5 @@
          END DO
       END DO
-      !
-      ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
-      DO jj = 2, jpjm1
+      DO jj = 2, jpjm1                             ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
          DO ji = fs_2, fs_jpim1   ! vector opt.
             en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1)
@@ -434,9 +456,6 @@
          END DO
       ELSEIF( nn_etau == 3 ) THEN       !* penetration belox the mixed layer (HF variability)
-!CDIR NOVERRCHK
          DO jk = 2, jpkm1
-!CDIR NOVERRCHK
-            DO jj = 2, jpjm1
-!CDIR NOVERRCHK
+            DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
@@ -453,7 +472,7 @@
       CALL lbc_lnk( en, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
       !
-      CALL wrk_dealloc( jpi,jpj, imlc )    ! integer
-      CALL wrk_dealloc( jpi,jpj, zhlc ) 
-      CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 
+      CALL wrk_dealloc( jpi,jpj,       imlc )    ! integer
+      CALL wrk_dealloc( jpi,jpj,       zhlc ) 
+      CALL wrk_dealloc( jpi,jpj,jpk,   zpelc, zdiag, zd_up, zd_lw, z3du, z3dv ) 
       !
       IF( nn_timing == 1 )  CALL timing_stop('tke_tke')
@@ -499,5 +518,5 @@
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       REAL(wp) ::   zrn2, zraug, zcoef, zav     ! local scalars
-      REAL(wp) ::   zdku, zpdlr, zri, zsqen     !   -      -
+      REAL(wp) ::   zdku, zri, zsqen            !   -      -
       REAL(wp) ::   zdkv, zemxl, zemlm, zemlp   !   -      -
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zmpdl, zmxlm, zmxld
@@ -529,12 +548,9 @@
       ENDIF
       !
-!CDIR NOVERRCHK
       DO jk = 2, jpkm1              ! interior value : l=sqrt(2*e/n^2)
-!CDIR NOVERRCHK
-         DO jj = 2, jpjm1
-!CDIR NOVERRCHK
+         DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
                zrn2 = MAX( rn2(ji,jj,jk), rsmall )
-               zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) )
+               zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 )  )
             END DO
          END DO
@@ -543,9 +559,10 @@
       !                     !* Physical limits for the mixing length
       !
-      zmxld(:,:,1  ) = zmxlm(:,:,1)   ! surface set to the minimum value 
+      zmxld(:,:, 1 ) = zmxlm(:,:,1)   ! surface set to the minimum value 
       zmxld(:,:,jpk) = rmxl_min       ! last level  set to the minimum value
       !
       SELECT CASE ( nn_mxl )
       !
+ !!gm Not sure of that coding for ISF....
       ! where wmask = 0 set zmxlm == fse3w
       CASE ( 0 )           ! bounded by the distance to surface and bottom
@@ -606,9 +623,6 @@
             END DO
          END DO
-!CDIR NOVERRCHK
          DO jk = 2, jpkm1
-!CDIR NOVERRCHK
-            DO jj = 2, jpjm1
-!CDIR NOVERRCHK
+            DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   zemlm = MIN ( zmxld(ji,jj,jk),  zmxlm(ji,jj,jk) )
@@ -630,9 +644,6 @@
       !                     !  Vertical eddy viscosity and diffusivity  (avmu, avmv, avt)
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-!CDIR NOVERRCHK
       DO jk = 1, jpkm1            !* vertical eddy viscosity & diffivity at w-points
-!CDIR NOVERRCHK
-         DO jj = 2, jpjm1
-!CDIR NOVERRCHK
+         DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
                zsqen = SQRT( en(ji,jj,jk) )
@@ -660,19 +671,10 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zcoef = avm(ji,jj,jk) * 2._wp * fse3w(ji,jj,jk) * fse3w(ji,jj,jk)
-                  !                                          ! shear
-                  zdku = avmu(ji-1,jj,jk) * ( un(ji-1,jj,jk-1) - un(ji-1,jj,jk) ) * ( ub(ji-1,jj,jk-1) - ub(ji-1,jj,jk) )   &
-                    &  + avmu(ji  ,jj,jk) * ( un(ji  ,jj,jk-1) - un(ji  ,jj,jk) ) * ( ub(ji  ,jj,jk-1) - ub(ji  ,jj,jk) )
-                  zdkv = avmv(ji,jj-1,jk) * ( vn(ji,jj-1,jk-1) - vn(ji,jj-1,jk) ) * ( vb(ji,jj-1,jk-1) - vb(ji,jj-1,jk) )   &
-                    &  + avmv(ji,jj  ,jk) * ( vn(ji,jj  ,jk-1) - vn(ji,jj  ,jk) ) * ( vb(ji,jj  ,jk-1) - vb(ji,jj  ,jk) )
-                  !                                          ! local Richardson number
-                  zri   = MAX( rn2b(ji,jj,jk), 0._wp ) * zcoef / (zdku + zdkv + rn_bshear )
-                  zpdlr = MAX(  0.1_wp,  0.2 / MAX( 0.2 , zri )  )
-!!gm and even better with the use of the "true" ri_crit=0.22222...  (this change the results!)
-!!gm              zpdlr = MAX(  0.1_wp,  ri_crit / MAX( ri_crit , zri )  )
-                  avt(ji,jj,jk)   = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk)
+                  avt(ji,jj,jk)   = MAX( apdlr(ji,jj,jk) * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
 # if defined key_c1d
-                  e_pdl(ji,jj,jk) = zpdlr * wmask(ji,jj,jk)  ! c1d configuration : save masked Prandlt number
-                  e_ric(ji,jj,jk) = zri   * wmask(ji,jj,jk)  ! c1d config. : save Ri
+                  e_pdl(ji,jj,jk) = apdlr(ji,jj,jk) * wmask(ji,jj,jk)    ! c1d configuration : save masked Prandlt number
+!!gm bug NO zri here....
+!!gm remove the specific diag for c1d !
+                  e_ric(ji,jj,jk) = zri * wmask(ji,jj,jk)                            ! c1d config. : save Ri
 # endif
               END DO
@@ -761,5 +763,5 @@
       IF( nn_pdl  < 0   .OR.  nn_pdl  > 1 )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' )
       IF( nn_htau < 0   .OR.  nn_htau > 1 )   CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' )
-      IF( nn_etau == 3 .AND. .NOT. lk_cpl )   CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' )
+      IF( nn_etau == 3 .AND. .NOT. ln_cpl )   CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' )
 
       IF( ln_mxl0 ) THEN
