Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/AMM12/EXP00/namelist_cfg
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/AMM12/EXP00/namelist_cfg	(revision 7989)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/AMM12/EXP00/namelist_cfg	(revision 7990)
@@ -324,7 +324,7 @@
       rn_hsbfr =    1.6       !  heat/salt buoyancy flux ratio
    !
-   ln_zdftmx   = .false.   ! tidal mixing parameterization              (T =>   fill namzdf_tmx)
-   !
-   ln_zdfqiao  = .false.   ! surface wave-induced mixing                (T => ln_wave=ln_sdw=T )
+   !                       ! gravity wave-driven vertical mixing
+   ln_zdfiwm   = .false.      ! internal wave-induced mixing            (T =>   fill namzdf_iwm)
+   ln_zdfswm   = .false.      ! surface  wave-induced mixing            (T => ln_wave=ln_sdw=T )
    !
    !                       ! time-stepping
@@ -353,5 +353,5 @@
 /
 !-----------------------------------------------------------------------
-&namzdf_tmx    !   internal wave-driven mixing parameterization         (ln_zdftmx =T)
+&namzdf_iwm    !   internal wave-driven mixing parameterization         (ln_zdfiwm =T)
 !-----------------------------------------------------------------------
 /
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/C1D_PAPA/EXP00/namelist_cfg
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/C1D_PAPA/EXP00/namelist_cfg	(revision 7989)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/C1D_PAPA/EXP00/namelist_cfg	(revision 7990)
@@ -275,10 +275,11 @@
       rn_hsbfr =    1.6       !  heat/salt buoyancy flux ratio
    !
-   ln_zdftmx   = .false.   ! tidal mixing parameterization              (T =>   fill namzdf_tmx)
-   !
-   ln_zdfqiao  = .false.   ! surface wave-induced mixing                (T => ln_wave=ln_sdw=T )
+   !                       ! gravity wave-driven vertical mixing
+   ln_zdfiwm   = .false.      ! internal wave-induced mixing            (T =>   fill namzdf_iwm)
+   ln_zdfswm   = .false.      ! surface  wave-induced mixing            (T => ln_wave=ln_sdw=T )
    !
    !                       ! time-stepping
    ln_zdfexp   = .false.      ! split-explicit (T) or implicit (F) scheme
+      nn_zdfexp=    3            !  number of sub-timestep for ln_zdfexp=T
    !
    !                       ! coefficients
@@ -301,5 +302,5 @@
 /
 !-----------------------------------------------------------------------
-&namzdf_tmx    !    internal wave-driven mixing parameterization        (ln_zdftmx =T)
+&namzdf_iwm    !    internal wave-driven mixing parameterization        (ln_zdfiwm =T)
 !-----------------------------------------------------------------------
 /
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/GYRE_BFM/EXP00/namelist_cfg
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/GYRE_BFM/EXP00/namelist_cfg	(revision 7989)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/GYRE_BFM/EXP00/namelist_cfg	(revision 7990)
@@ -264,7 +264,7 @@
       rn_hsbfr =    1.6       !  heat/salt buoyancy flux ratio
    !
-   ln_zdftmx   = .false.   ! tidal mixing parameterization              (T =>   fill namzdf_tmx)
-   !
-   ln_zdfqiao  = .false.   ! surface wave-induced mixing                (T => ln_wave=ln_sdw=T )
+   !                       ! gravity wave-driven vertical mixing
+   ln_zdfiwm   = .false.      ! internal wave-induced mixing            (T =>   fill namzdf_iwm)
+   ln_zdfswm   = .false.      ! surface  wave-induced mixing            (T => ln_wave=ln_sdw=T )
    !
    !                       ! time-stepping
@@ -292,5 +292,5 @@
 /
 !-----------------------------------------------------------------------
-&namzdf_tmx    !    internal wave-driven mixing parameterization        (ln_zdftmx =T)
+&namzdf_iwm    !    internal wave-driven mixing parameterization        (ln_zdfiwm =T)
 !-----------------------------------------------------------------------
 /
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/GYRE_PISCES/EXP00/namelist_cfg
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/GYRE_PISCES/EXP00/namelist_cfg	(revision 7989)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/GYRE_PISCES/EXP00/namelist_cfg	(revision 7990)
@@ -218,11 +218,11 @@
       rn_hsbfr =    1.6       !  heat/salt buoyancy flux ratio
    !
-   ln_zdftmx   = .false.   ! tidal mixing parameterization              (T =>   fill namzdf_tmx)
-   !
-   ln_zdfqiao  = .false.   ! enhanced wave vertical mixing Qiao (2010) (T => ln_wave=T & ln_sdw=T & fill namsbc_wave)
+   !                       ! gravity wave-driven vertical mixing
+   ln_zdfiwm   = .false.      ! internal wave-induced mixing            (T =>   fill namzdf_iwm)
+   ln_zdfswm   = .false.      ! surface  wave-induced mixing            (T => ln_wave=ln_sdw=T )
    !
    !                       ! time-stepping
-   ln_zdfexp   = .false.   ! split-explicit (T) or implicit (F) time stepping scheme
-      nn_zdfexp=    3         !  number of sub-timestep for ln_zdfexp=T
+   ln_zdfexp   = .false.      ! split-explicit (T) or implicit (F) scheme
+      nn_zdfexp=    3            !  number of sub-timestep for ln_zdfexp=T
    !
    !                       !  Coefficients
@@ -251,9 +251,4 @@
 /
 !-----------------------------------------------------------------------
-&namzdf_tmx    !   tidal mixing parameterization                        (ln_zdftmx=T)
-!-----------------------------------------------------------------------
-   ln_tmx_itf  = .false.   !  ITF specific parameterisation
-/
-!-----------------------------------------------------------------------
 &nammpp        !   Massively Parallel Processing                        ("key_mpp_mpi)
 !-----------------------------------------------------------------------
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/ORCA2_LIM3_PISCES/EXP00/1_namelist_cfg
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/ORCA2_LIM3_PISCES/EXP00/1_namelist_cfg	(revision 7989)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/ORCA2_LIM3_PISCES/EXP00/1_namelist_cfg	(revision 7990)
@@ -226,11 +226,11 @@
       rn_hsbfr =    1.6       !  heat/salt buoyancy flux ratio
    !
-   ln_zdftmx   = .true.    ! tidal mixing parameterization              (T =>   fill namzdf_tmx)
-   !
-   ln_zdfqiao  = .false.   ! enhanced wave vertical mixing Qiao (2010) (T => ln_wave=T & ln_sdw=T & fill namsbc_wave)
+   !                       ! gravity wave-driven vertical mixing
+   ln_zdfiwm   = .true.       ! internal wave-induced mixing            (T =>   fill namzdf_iwm)
+   ln_zdfswm   = .false.      ! surface  wave-induced mixing            (T => ln_wave=ln_sdw=T )
    !
    !                       ! time-stepping
-   ln_zdfexp   = .false.   ! split-explicit (T) or implicit (F) time stepping scheme
-      nn_zdfexp=    3         !  number of sub-timestep for ln_zdfexp=T
+   ln_zdfexp   = .false.      ! split-explicit (T) or implicit (F) scheme
+      nn_zdfexp=    3            !  number of sub-timestep for ln_zdfexp=T
    !
    !                       !  Coefficients
@@ -245,5 +245,5 @@
 /
 !-----------------------------------------------------------------------
-&namzdf_tmx    !    internal wave-driven mixing parameterization        (ln_zdftmx =T)
+&namzdf_iwm    !    internal wave-driven mixing parameterization        (ln_zdfiwm =T)
 !-----------------------------------------------------------------------
    nn_zpyc     = 1         !  pycnocline-intensified dissipation scales as N (=1) or N^2 (=2)
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/ORCA2_LIM3_PISCES/EXP00/namelist_cfg
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/ORCA2_LIM3_PISCES/EXP00/namelist_cfg	(revision 7989)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/ORCA2_LIM3_PISCES/EXP00/namelist_cfg	(revision 7990)
@@ -27,6 +27,5 @@
 /
 !-----------------------------------------------------------------------
-&namcrs        !   Grid coarsening for dynamics output and/or
-               !   passive tracer coarsened online simulations
+&namcrs        !   coarsened grid (for outputs and/or TOP)              (ln_crs =T)
 !-----------------------------------------------------------------------
 /
@@ -255,11 +254,11 @@
       rn_hsbfr =    1.6       !  heat/salt buoyancy flux ratio
    !
-   ln_zdftmx   = .true.    ! tidal mixing parameterization              (T =>   fill namzdf_tmx)
-   !
-   ln_zdfqiao  = .false.   ! enhanced wave vertical mixing Qiao (2010) (T => ln_wave=T & ln_sdw=T & fill namsbc_wave)
+   !                       ! gravity wave-driven vertical mixing
+   ln_zdfiwm   = .false.      ! internal wave-induced mixing            (T =>   fill namzdf_iwm)
+   ln_zdfswm   = .false.      ! surface  wave-induced mixing            (T => ln_wave=ln_sdw=T )
    !
    !                       ! time-stepping
-   ln_zdfexp   = .false.   ! split-explicit (T) or implicit (F) time stepping scheme
-      nn_zdfexp=    3         !  number of sub-timestep for ln_zdfexp=T
+   ln_zdfexp   = .false.      ! split-explicit (T) or implicit (F) scheme
+      nn_zdfexp=    3            !  number of sub-timestep for ln_zdfexp=T
    !
    !                       !  Coefficients
@@ -274,5 +273,5 @@
 /
 !-----------------------------------------------------------------------
-&namzdf_tmx    !   tidal mixing parameterization                        (ln_zdftmx =T)
+&namzdf_iwm    !   tidal mixing parameterization                        (ln_zdfiwm =T)
 !-----------------------------------------------------------------------
    nn_zpyc     = 2         !  pycnocline-intensified dissipation scales as N (=1) or N^2 (=2)
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/SHARED/field_def_nemo-opa.xml
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/SHARED/field_def_nemo-opa.xml	(revision 7989)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/SHARED/field_def_nemo-opa.xml	(revision 7990)
@@ -398,10 +398,10 @@
          <field id="avm_evd"      long_name="convective enhancement of vertical viscosity"     standard_name="ocean_vertical_momentum_diffusivity_due_to_convection"   unit="m2/s" />
 
-         <!-- variables available with ln_zdftmx =T -->
+         <!-- variables available with ln_zdfiwm =T -->
          <field id="av_ratio"     long_name="S over T diffusivity ratio"            standard_name="salinity_over_temperature_diffusivity_ratio"                     unit="1"    />
          <field id="av_wave"      long_name="wave-induced vertical diffusivity"     standard_name="ocean_vertical_tracer_diffusivity_due_to_internal_waves"         unit="m2/s" />
-         <field id="bflx_tmx"     long_name="wave-induced buoyancy flux"            standard_name="buoyancy_flux_due_to_internal_waves"                             unit="W/kg" />
-         <field id="pcmap_tmx"    long_name="power consumed by wave-driven mixing"  standard_name="vertically_integrated_power_consumption_by_wave_driven_mixing"   unit="W/m2"      grid_ref="grid_W_2D" />
-         <field id="emix_tmx"     long_name="power density available for mixing"    standard_name="power_available_for_mixing_from_breaking_internal_waves"         unit="W/kg" />
+         <field id="bflx_iwm"     long_name="wave-induced buoyancy flux"            standard_name="buoyancy_flux_due_to_internal_waves"                             unit="W/kg" />
+         <field id="pcmap_iwm"    long_name="power consumed by wave-driven mixing"  standard_name="vertically_integrated_power_consumption_by_wave_driven_mixing"   unit="W/m2"      grid_ref="grid_W_2D" />
+         <field id="emix_iwm"     long_name="power density available for mixing"    standard_name="power_available_for_mixing_from_breaking_internal_waves"         unit="W/kg" />
 
          <!-- variables available with diaar5 -->   
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 7989)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 7990)
@@ -11,5 +11,5 @@
 !!              6 - Tracer           (nameos, namtra_adv, namtra_ldf, namtra_ldfeiv, namtra_dmp)
 !!              7 - dynamics         (namdyn_adv, namdyn_vor, namdyn_hpg, namdyn_spg, namdyn_ldf)
-!!              8 - Verical physics  (namzdf, namzdf_ric, namzdf_tke, namzdf_gls, namzdf_tmx, namzdf_tmx_new)
+!!              8 - Verical physics  (namzdf, namzdf_ric, namzdf_tke, namzdf_gls, namzdf_iwm)
 !!              9 - miscellaneous    (nammpp, namctl)
 !!             10 - diagnostics      (namnc4, namtrd, namspr, namflo, namhsb, namsto)
@@ -875,5 +875,5 @@
 !!    namzdf_tke    TKE vertical mixing                                 (ln_zdftke=T)
 !!    namzdf_gls    GLS vertical mixing                                 (ln_zdfgls=T)
-!!    namzdf_tmx    tidal mixing parameterization                       (ln_zdftmx=T)
+!!    namzdf_iwm    tidal mixing parameterization                       (ln_zdfiwm=T)
 !!======================================================================
 !
@@ -899,7 +899,7 @@
       rn_hsbfr =    1.6       !  heat/salt buoyancy flux ratio
    !
-   ln_zdftmx   = .false.   ! tidal mixing parameterization              (T =>   fill namzdf_tmx)
-   !
-   ln_zdfqiao  = .false.   ! surface wave-induced mixing                (T => ln_wave=ln_sdw=T )
+   !                       ! gravity wave-driven vertical mixing
+   ln_zdfiwm   = .false.      ! internal wave-induced mixing            (T =>   fill namzdf_iwm)
+   ln_zdfswm   = .false.      ! surface  wave-induced mixing            (T => ln_wave=ln_sdw=T )
    !
    !                       ! time-stepping
@@ -919,10 +919,10 @@
    rn_alp      =    5.     !  coefficient of the parameterization
    nn_ric      =    2      !  coefficient of the parameterization
-   rn_ekmfc    =    0.7    !  Factor in the Ekman depth Equation
-   rn_mldmin   =    1.0    !  minimum allowable mixed-layer depth estimate (m)
-   rn_mldmax   = 1000.0    !  maximum allowable mixed-layer depth estimate (m)
-   rn_wtmix    =   10.0    !  vertical eddy viscosity coeff [m2/s] in the mixed-layer
-   rn_wvmix    =   10.0    !  vertical eddy diffusion coeff [m2/s] in the mixed-layer
-   ln_mldw     =  .true.   !  Flag to use or not the mixed layer depth param.
+   ln_mldw     =  .false.  !  enhanced mixing in the Ekman layer
+      rn_ekmfc    =    0.7    !  Factor in the Ekman depth Equation
+      rn_mldmin   =    1.0    !  minimum allowable mixed-layer depth estimate (m)
+      rn_mldmax   = 1000.0    !  maximum allowable mixed-layer depth estimate (m)
+      rn_wtmix    =   10.0    !  vertical eddy viscosity coeff [m2/s] in the mixed-layer
+      rn_wvmix    =   10.0    !  vertical eddy diffusion coeff [m2/s] in the mixed-layer
 /
 !-----------------------------------------------------------------------
@@ -974,5 +974,5 @@
 /
 !-----------------------------------------------------------------------
-&namzdf_tmx    !    internal wave-driven mixing parameterization        (ln_zdftmx =T)
+&namzdf_iwm    !    internal wave-driven mixing parameterization        (ln_zdfiwm =T)
 !-----------------------------------------------------------------------
    nn_zpyc     = 1         !  pycnocline-intensified dissipation scales as N (=1) or N^2 (=2)
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/ISOMIP/EXP00/namelist_cfg
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/ISOMIP/EXP00/namelist_cfg	(revision 7989)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/ISOMIP/EXP00/namelist_cfg	(revision 7990)
@@ -344,7 +344,7 @@
    ln_zdfddm   = .false.   ! double diffusive mixing
    !
-   ln_zdftmx   = .false.   ! tidal mixing parameterization              (T =>   fill namzdf_tmx)
-   !
-   ln_zdfqiao  = .false.   ! surface wave-induced mixing                (T => ln_wave=ln_sdw=T )
+   !                       ! gravity wave-driven vertical mixing
+   ln_zdfiwm   = .false.      ! internal wave-induced mixing            (T =>   fill namzdf_iwm)
+   ln_zdfswm   = .false.      ! surface  wave-induced mixing            (T => ln_wave=ln_sdw=T )
    !
    !                       ! time-stepping
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/SAS_BIPER/EXP00/namelist_cfg
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/SAS_BIPER/EXP00/namelist_cfg	(revision 7989)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/SAS_BIPER/EXP00/namelist_cfg	(revision 7990)
@@ -223,5 +223,5 @@
 /
 !-----------------------------------------------------------------------
-&namzdf_tmx    !   tidal mixing parameterization                        (ln_zdftmx =T)
+&namzdf_iwm    !   tidal mixing parameterization                        (ln_zdfiwm =T)
 !-----------------------------------------------------------------------
 /
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/WAD/EXP00/namelist_cfg
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/WAD/EXP00/namelist_cfg	(revision 7989)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/WAD/EXP00/namelist_cfg	(revision 7990)
@@ -328,5 +328,5 @@
 !!    namzdf_tke    TKE vertical mixing                                 (ln_zdftke=T)
 !!    namzdf_gls    GLS vertical mixing                                 (ln_zdfgls=T)
-!!    namzdf_tmx    tidal mixing parameterization                       (ln_zdftmx=T)
+!!    namzdf_iwm    tidal mixing parameterization                       (ln_zdfiwm=T)
 !!======================================================================
 !-----------------------------------------------------------------------
@@ -351,10 +351,11 @@
       rn_hsbfr =    1.6       !  heat/salt buoyancy flux ratio
    !
-   ln_zdftmx   = .false.   ! tidal mixing parameterization              (T =>   fill namzdf_tmx)
-   !
-   ln_zdfqiao  = .false.   ! surface wave-induced mixing (Qiao et al. 2010) (T =>   ln_wave=ln_sdw=T. & fill namsbc_wave)
+   !                       ! gravity wave-driven vertical mixing
+   ln_zdfiwm   = .false.      ! internal wave-induced mixing            (T =>   fill namzdf_iwm)
+   ln_zdfswm   = .false.      ! surface  wave-induced mixing            (T => ln_wave=ln_sdw=T )
    !
    !                       ! time-stepping
-   ln_zdfexp   = .false.   ! split-explicit (T) or implicit (F) scheme
+   ln_zdfexp   = .false.      ! split-explicit (T) or implicit (F) scheme
+      nn_zdfexp=    3            !  number of sub-timestep for ln_zdfexp=T
    !
    !                       ! coefficients
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/cfg.txt
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/cfg.txt	(revision 7989)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/cfg.txt	(revision 7990)
@@ -6,4 +6,4 @@
 ORCA2_OFF_TRC OPA_SRC OFF_SRC TOP_SRC
 ORCA2_LIM3_PISCES OPA_SRC LIM_SRC_3 TOP_SRC NST_SRC
+GYRE_PISCES OPA_SRC TOP_SRC
 GYRE_PISCES_XIOS OPA_SRC TOP_SRC
-GYRE_PISCES OPA_SRC TOP_SRC
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90	(revision 7989)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90	(revision 7990)
@@ -9,5 +9,5 @@
    USE dom_oce         ! ocean space and time domain
    USE zdf_oce         ! ocean vertical physics    
-   USE zdfgls, ONLY: mxln
+   USE zdfgls   , ONLY : hmxn
    USE in_out_manager  ! I/O units
    USE iom             ! I/0 library
@@ -107,5 +107,5 @@
       IF( ln_zdfgls ) THEN
          en_25h(:,:,:) = en(:,:,:)
-         rmxln_25h(:,:,:) = mxln(:,:,:)
+         rmxln_25h(:,:,:) = hmxn(:,:,:)
       ENDIF
 #if defined key_lim3 || defined key_lim2
@@ -169,5 +169,5 @@
          IF( ln_zdfgls ) THEN
             en_25h(:,:,:)        = en_25h(:,:,:) + en(:,:,:)
-            rmxln_25h(:,:,:)      = rmxln_25h(:,:,:) + mxln(:,:,:)
+            rmxln_25h(:,:,:)      = rmxln_25h(:,:,:) + hmxn(:,:,:)
          ENDIF
          cnt_25h = cnt_25h + 1
@@ -217,19 +217,19 @@
          zw3d(:,:,:) = vn_25h(:,:,:)*vmask(:,:,:) + zmdi*(1.0-vmask(:,:,:))
          CALL iom_put("vomecrty25h", zw3d  )   ! j-current
-         zw3d(:,:,:) = wn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+         zw3d(:,:,:) = wn_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
          CALL iom_put("vomecrtz25h", zw3d )   ! k-current
          ! Write vertical physics
-         zw3d(:,:,:) = avt_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+         zw3d(:,:,:) = avt_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
          CALL iom_put("avt25h", zw3d )   ! diffusivity
-         zw3d(:,:,:) = avm_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+         zw3d(:,:,:) = avm_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
          CALL iom_put("avm25h", zw3d)   ! viscosity
          IF( ln_zdftke ) THEN
-            zw3d(:,:,:) = en_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+            zw3d(:,:,:) = en_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
             CALL iom_put("tke25h", zw3d)   ! tke
          ENDIF
          IF( ln_zdfgls ) THEN
-            zw3d(:,:,:) = en_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+            zw3d(:,:,:) = en_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
             CALL iom_put("tke25h", zw3d)   ! tke
-            zw3d(:,:,:) = rmxln_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+            zw3d(:,:,:) = rmxln_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
             CALL iom_put( "mxln25h",zw3d)
          ENDIF
@@ -249,5 +249,5 @@
          IF( ln_zdfgls ) THEN
             en_25h(:,:,:) = en(:,:,:)
-            rmxln_25h(:,:,:) = mxln(:,:,:)
+            rmxln_25h(:,:,:) = hmxn(:,:,:)
          ENDIF
          cnt_25h = 1
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90	(revision 7989)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90	(revision 7990)
@@ -42,5 +42,4 @@
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
    
@@ -51,5 +50,4 @@
       !! ** Purpose :   compute the vertical ocean dynamics physics.
       !!---------------------------------------------------------------------
-      !!
       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
       !
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 7989)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 7990)
@@ -12,5 +12,5 @@
 
    !!----------------------------------------------------------------------
-   !!   dyn_zdf_imp   : compute the vertical diffusion using a implicit scheme
+   !!   dyn_zdf_imp   : compute the vertical diffusion using a implicit time-stepping
    !!                   together with the Leap-Frog time integration.
    !!----------------------------------------------------------------------
@@ -39,5 +39,5 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 4.0 , NEMO Consortium (2017)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 7989)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 7990)
@@ -9,4 +9,5 @@
    !!            3.4  ! 2011_11  (C. Harris) more flexibility + multi-category fields
    !!----------------------------------------------------------------------
+
    !!----------------------------------------------------------------------
    !!   namsbc_cpl      : coupled formulation namlist
@@ -974,5 +975,5 @@
       !!                        emp          upward mass flux [evap. - precip. (- runoffs) (- calving)] (ocean only case)
       !!----------------------------------------------------------------------
-      USE zdf_oce,  ONLY : ln_zdfqiao
+      USE zdf_oce,  ONLY : ln_zdfswm
 
       IMPLICIT NONE
@@ -1159,15 +1160,15 @@
       !                                                      !      Wave mean period     !
       !                                                      ! ========================= ! 
-         IF( srcv(jpr_wper)%laction ) wmp(:,:) = frcv(jpr_wper)%z3(:,:,1)
+         IF( srcv(jpr_wper)%laction )   wmp(:,:) = frcv(jpr_wper)%z3(:,:,1)
       !
       !                                                      ! ========================= ! 
       !                                                      !  Significant wave height  !
       !                                                      ! ========================= ! 
-         IF( srcv(jpr_hsig)%laction ) hsw(:,:) = frcv(jpr_hsig)%z3(:,:,1)
+         IF( srcv(jpr_hsig)%laction )   hsw(:,:) = frcv(jpr_hsig)%z3(:,:,1)
       !
       !                                                      ! ========================= ! 
-      !                                                      !    Vertical mixing Qiao   !
+      !                                                      !    surface wave mixing    !
       !                                                      ! ========================= ! 
-         IF( srcv(jpr_wnum)%laction .AND. ln_zdfqiao ) wnum(:,:) = frcv(jpr_wnum)%z3(:,:,1)
+         IF( srcv(jpr_wnum)%laction .AND. ln_zdfswm )   wnum(:,:) = frcv(jpr_wnum)%z3(:,:,1)
 
          ! Calculate the 3D Stokes drift both in coupled and not fully uncoupled mode
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 7989)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 7990)
@@ -19,5 +19,5 @@
    USE oce            ! ocean variables
    USE sbc_oce	       ! Surface boundary condition: ocean fields
-   USE zdf_oce,  ONLY : ln_zdfqiao
+   USE zdf_oce,  ONLY : ln_zdfswm
    USE bdy_oce        ! open boundary condition variables
    USE domvvl         ! domain: variable volume layers
@@ -227,5 +227,5 @@
          !
          ! Read also wave number if needed, so that it is available in coupling routines
-         IF( ln_zdfqiao .AND. .NOT.cpl_wnum ) THEN
+         IF( ln_zdfswm .AND. .NOT.cpl_wnum ) THEN
             CALL fld_read( kt, nn_fsbc, sf_wn )          ! read wave parameters from external forcing
             wnum(:,:) = sf_wn(1)%fnow(:,:,1)
@@ -345,5 +345,5 @@
          vsd(:,:,:) = 0._wp
          wsd(:,:,:) = 0._wp
-         ! Wave number needed only if ln_zdfqiao=T
+         ! Wave number needed only if ln_zdfswm=T
          IF( .NOT. cpl_wnum ) THEN
             ALLOCATE( sf_wn(1), STAT=ierror )           !* allocate and fill sf_wave with sn_wnum
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90	(revision 7989)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90	(revision 7990)
@@ -16,6 +16,5 @@
    PUBLIC  zdf_oce_alloc    ! Called in nemogcm.F90
 
-   !                                 !!* namelist namzdf: vertical diffusion *
-!!gm
+   !                            !!* namelist namzdf: vertical physics *
    !                             ! vertical closure scheme flags
    LOGICAL , PUBLIC ::   ln_zdfcst   !: constant coefficients
@@ -23,13 +22,5 @@
    LOGICAL , PUBLIC ::   ln_zdftke   !: Turbulent Kinetic Energy closure
    LOGICAL , PUBLIC ::   ln_zdfgls   !: Generic Length Sclare closure
-   !                             ! tidal induced mixing   
-   LOGICAL , PUBLIC ::   ln_zdftmx   !: tidal mixing parameterization flag
-   !                             ! double diffusion
-   LOGICAL , PUBLIC ::   ln_zdfddm   !: double diffusive mixing flag
-   REAL(wp), PUBLIC ::      rn_avts     !: maximum value of avs for salt fingering
-   REAL(wp), PUBLIC ::      rn_hsbfr    !: heat/salt buoyancy flux ratio
-!!gm
-   LOGICAL , PUBLIC ::   ln_zdfexp   !: explicit vertical diffusion scheme flag
-   INTEGER , PUBLIC ::   nn_zdfexp   !: number of sub-time step (explicit time stepping)
+   !                             ! convection
    LOGICAL , PUBLIC ::   ln_zdfevd   !: convection: enhanced vertical diffusion flag
    INTEGER , PUBLIC ::      nn_evdm     !: =0/1 flag to apply enhanced avm or not
@@ -38,6 +29,14 @@
    INTEGER , PUBLIC ::      nn_npc      !: non penetrative convective scheme call  frequency
    INTEGER , PUBLIC ::      nn_npcp     !: non penetrative convective scheme print frequency
-   !                             ! Surface wave-induced mixing
-   LOGICAL , PUBLIC ::   ln_zdfqiao  !: Enhanced wave vertical mixing Qiao(2010) formulation flag
+   !                             ! double diffusion
+   LOGICAL , PUBLIC ::   ln_zdfddm   !: double diffusive mixing flag
+   REAL(wp), PUBLIC ::      rn_avts     !: maximum value of avs for salt fingering
+   REAL(wp), PUBLIC ::      rn_hsbfr    !: heat/salt buoyancy flux ratio
+   !                             ! gravity wave-induced vertical mixing
+   LOGICAL , PUBLIC ::   ln_zdfswm   !: surface  wave-induced mixing flag
+   LOGICAL , PUBLIC ::   ln_zdfiwm   !: internal wave-induced mixing flag
+   !                             ! time-stepping
+   LOGICAL , PUBLIC ::   ln_zdfexp   !: explicit vertical diffusion scheme flag
+   INTEGER , PUBLIC ::      nn_zdfexp   !: number of sub-time step (explicit time stepping)
    !                             ! coefficients 
    REAL(wp), PUBLIC ::   rn_avm0     !: vertical eddy viscosity (m2/s)
@@ -47,15 +46,13 @@
 
 
-   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avt  , avs     !: vertical eddy diffusivity coef at w-pt [m2/s]
-
+   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avt  , avs     !: vertical T & S  diffusivity coef at w-pt [m2/s]
+   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avm            !: vertical momentum viscosity coef at w-pt [m2/s]
+   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avmu , avmv    !: vertical viscosity coef at uw- & vw-pts  [m2/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avt_k , avm_k  !: avt, avs computed by turbulent closure alone
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   en             !: now turbulent kinetic energy   [m2/s2]
    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::   avmb , avtb    !: background profile of avm and avt
    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   avtb_2d        !: horizontal shape of background Kz profile
-   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   bfrua, bfrva   !: Bottom friction coefficients set in zdfbfr
-   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   tfrua, tfrva   !: top friction coefficients set in zdfbfr
-   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avmu , avmv    !: vertical viscosity coef at uw- & vw-pts       [m2/s]
-   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avm            !: vertical viscosity & diffusivity coef at w-pt [m2/s]
-   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), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   en              !: now turbulent kinetic energy   [m2/s2]
+   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   bfrua, bfrva   !: bottom friction coefficients
+   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   tfrua, tfrva   !: top    friction coefficients
 
    !!----------------------------------------------------------------------
@@ -71,12 +68,9 @@
       !!----------------------------------------------------------------------
       !
-      ALLOCATE(avmb(jpk) , bfrua(jpi,jpj) ,                         &
-         &     avtb(jpk) , bfrva(jpi,jpj) , avtb_2d(jpi,jpj) ,      &
-         &     tfrua(jpi, jpj), tfrva(jpi, jpj)              ,      &
-         &     avmu  (jpi,jpj,jpk), avm   (jpi,jpj,jpk)      ,      &
-         &     avmv  (jpi,jpj,jpk), avt   (jpi,jpj,jpk)      , avs   (jpi,jpj,jpk),          &
-         &     avt_k (jpi,jpj,jpk), avm_k (jpi,jpj,jpk)      ,      & 
-         &     avmu_k(jpi,jpj,jpk), avmv_k(jpi,jpj,jpk)      ,      &
-         &     en    (jpi,jpj,jpk), STAT = zdf_oce_alloc )
+      ALLOCATE( avm  (jpi,jpj,jpk) , avt  (jpi,jpj,jpk) , avs(jpi,jpj,jpk) ,        &
+         &      avm_k(jpi,jpj,jpk) , avt_k(jpi,jpj,jpk) , en (jpi,jpj,jpk) ,        & 
+         &      avmb(jpk) , bfrua(jpi,jpj) , tfrua(jpi, jpj) ,                      &
+         &      avtb(jpk) , bfrva(jpi,jpj) , tfrva(jpi, jpj) , avtb_2d(jpi,jpj) ,   &
+         &      avmu(jpi,jpj,jpk), avmv(jpi,jpj,jpk) ,  STAT = zdf_oce_alloc )
          !
       IF( zdf_oce_alloc /= 0 )   CALL ctl_warn('zdf_oce_alloc: failed to allocate arrays')
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90	(revision 7989)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90	(revision 7990)
@@ -8,4 +8,5 @@
    !!            3.3  ! 2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
    !!            3.6  ! 2013-04  (G. Madec, F. Roquet) zrau compute locally using interpolation of alpha & beta
+   !!            4.0  !  2017-04  (G. Madec)  remove CPP ddm key & avm at t-point only 
    !!----------------------------------------------------------------------
 
@@ -22,5 +23,4 @@
    USE prtctl         ! Print control
    USE lib_mpp        ! MPP library
-   USE wrk_nemo       ! work arrays
    USE timing         ! Timing
 
@@ -63,5 +63,5 @@
       !!      avt = avt + zavft + zavdt
       !!      avs = avs + zavfs + zavds
-      !!      avmu, avmv are required to remain at least above avt and avs.
+      !!      avm is required to remain at least above avt and avs.
       !!      
       !! ** Action  :   avt, avs : updated vertical eddy diffusivity coef. for T & S
@@ -77,10 +77,8 @@
       REAL(wp) ::   zavft, zavfs    !   -      -
       REAL(wp) ::   zavdt, zavds    !   -      -
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3
+      REAL(wp), DIMENSION(jpi,jpj) ::   zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('zdf_ddm')
-      !
-      CALL wrk_alloc( jpi,jpj, zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3 )
+      IF( nn_timing == 1 )   CALL timing_start('zdf_ddm')
       !
       !                                                ! ===============
@@ -89,7 +87,12 @@
          ! Define the mask 
          ! ---------------
+!!gm  WORK to be done:   change the code from vector optimisation to scalar one.
+!!gm                     ==>>>  test in the loop instead of use of mask arrays
+!!gm                            and many acces in memory
+         
          DO jj = 1, jpj                                ! R=zrau = (alpha / beta) (dk[t] / dk[s])
             DO ji = 1, jpi
                zrw =   ( gdepw_n(ji,jj,jk  ) - gdept_n(ji,jj,jk) )   &
+!!gm please, use e3w_n below 
                   &  / ( gdept_n(ji,jj,jk-1) - gdept_n(ji,jj,jk) ) 
                !
@@ -156,36 +159,11 @@
             END DO
          END DO
-
-
-         ! Increase avmu, avmv if necessary
-         ! --------------------------------
-!!gm to be changed following the definition of avm.
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk),    &
-                  &                  avt(ji,jj,jk), avt(ji+1,jj,jk),   &
-                  &                  avs(ji,jj,jk), avs(ji+1,jj,jk) )  * wumask(ji,jj,jk)
-               avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk),    &
-                  &                  avt(ji,jj,jk), avt(ji,jj+1,jk),   &
-                  &                  avs(ji,jj,jk), avs(ji,jj+1,jk) )  * wvmask(ji,jj,jk)
-            END DO
-         END DO
          !                                                ! ===============
       END DO                                              !   End of slab
       !                                                   ! ===============
       !
-      CALL lbc_lnk( avt , 'W', 1._wp )     ! Lateral boundary conditions   (unchanged sign)
-      CALL lbc_lnk( avs , 'W', 1._wp )
-      CALL lbc_lnk( avm , 'W', 1._wp )
-      CALL lbc_lnk( avmu, 'U', 1._wp ) 
-      CALL lbc_lnk( avmv, 'V', 1._wp )
-
       IF(ln_ctl) THEN
          CALL prt_ctl(tab3d_1=avt , clinfo1=' ddm  - t: ', tab3d_2=avs , clinfo2=' s: ', ovlap=1, kdim=jpk)
-         CALL prt_ctl(tab3d_1=avmu, clinfo1=' ddm  - u: ', mask1=umask, &
-            &         tab3d_2=avmv, clinfo2=       ' v: ', mask2=vmask, ovlap=1, kdim=jpk)
       ENDIF
-      !
-      CALL wrk_dealloc( jpi,jpj, zrau, zmsks, zmskf, zmskd1, zmskd2, zmskd3 )
       !
       IF( nn_timing == 1 )  CALL timing_stop('zdf_ddm')
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90	(revision 7989)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90	(revision 7990)
@@ -8,4 +8,5 @@
    !!   NEMO     1.0  !  2002-06  (G. Madec)  F90: Free form and module
    !!            3.2  !  2009-03  (M. Leclair, G. Madec, R. Benshila) test on both before & after
+   !!            4.0  !  2017-04  (G. Madec)  evd applied on avm (at t-point) 
    !!----------------------------------------------------------------------
 
@@ -23,5 +24,4 @@
    USE iom             ! for iom_put
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE wrk_nemo        ! work arrays
    USE timing          ! Timing
 
@@ -45,16 +45,18 @@
       !!      sivity coefficients when a static instability is encountered.
       !!
-      !! ** Method  :   avt, avm, and the 4 neighbouring avmu, avmv coefficients
-      !!      are set to avevd (namelist parameter) if the water column is 
-      !!      statically unstable (i.e. if rn2 < -1.e-12 )
+      !! ** Method  :   tracer (and momentum if nn_evdm=1) vertical mixing 
+      !!              coefficients are set to rn_evd (namelist parameter) 
+      !!              if the water column is statically unstable.
+      !!                The test of static instability is performed using
+      !!              Brunt-Vaisala frequency (rn2 < -1.e-12) of to successive
+      !!              time-step (Leap-Frog environnement): before and
+      !!              now time-step.
       !!
-      !! ** Action  :   avt, avm, avmu, avmv updted in static instability cases
-      !!
-      !! References :   Lazar, A., these de l'universite Paris VI, France, 1997
+      !! ** Action  :   avt, avm   enhanced where static instability occurs
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt   ! ocean time-step indexocean time step
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step indexocean time step
       !
       INTEGER ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zavt_evd, zavm_evd
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zavt_evd, zavm_evd
       !!----------------------------------------------------------------------
       !
@@ -68,5 +70,4 @@
       ENDIF
       !
-      CALL wrk_alloc( jpi,jpj,jpk,   zavt_evd, zavm_evd ) 
       !
       zavt_evd(:,:,:) = avt(:,:,:)           ! set avt prior to evd application
@@ -74,35 +75,39 @@
       SELECT CASE ( nn_evdm )
       !
-      CASE ( 1 )           ! enhance vertical eddy viscosity and diffusivity (if rn2<-1.e-12)
+      CASE ( 1 )           !==  enhance tracer & momentum Kz  ==!   (if rn2<-1.e-12)
          !
          zavm_evd(:,:,:) = avm(:,:,:)           ! set avm prior to evd application
          !
+!! change last digits results
+!         WHERE( MAX( rn2(2:jpi,2:jpj,2:jpkm1), rn2b(2:jpi,2:jpj,2:jpkm1) )  <= -1.e-12 ) THEN
+!            avt(2:jpi,2:jpj,2:jpkm1) = rn_evd * wmask(2:jpi,2:jpj,2:jpkm1)
+!            avm(2:jpi,2:jpj,2:jpkm1) = rn_evd * wmask(2:jpi,2:jpj,2:jpkm1)
+!         END WHERE
+         !
          DO jk = 1, jpkm1 
-            DO jj = 2, jpj             ! no vector opt.
-               DO ji = 2, jpi
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1
                   IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) THEN
-                     avt (ji  ,jj  ,jk) = rn_evd * tmask(ji  ,jj  ,jk)
-                     avm (ji  ,jj  ,jk) = rn_evd * tmask(ji  ,jj  ,jk)
-                     avmu(ji  ,jj  ,jk) = rn_evd * umask(ji  ,jj  ,jk)
-                     avmu(ji-1,jj  ,jk) = rn_evd * umask(ji-1,jj  ,jk)
-                     avmv(ji  ,jj  ,jk) = rn_evd * vmask(ji  ,jj  ,jk)
-                     avmv(ji  ,jj-1,jk) = rn_evd * vmask(ji  ,jj-1,jk)
+                     avt(ji,jj,jk) = rn_evd * wmask(ji,jj,jk)
+                     avm(ji,jj,jk) = rn_evd * wmask(ji,jj,jk)
                   ENDIF
                END DO
             END DO
          END DO 
-         CALL lbc_lnk( avt , 'W', 1. )   ;   CALL lbc_lnk( avm , 'W', 1. )   ! Lateral boundary conditions
-         CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )
          !
          zavm_evd(:,:,:) = avm(:,:,:) - zavm_evd(:,:,:)   ! change in avm due to evd
          CALL iom_put( "avm_evd", zavm_evd )              ! output this change
          !
-      CASE DEFAULT         ! enhance vertical eddy diffusivity only (if rn2<-1.e-12) 
+      CASE DEFAULT         !==  enhance tracer Kz  ==!   (if rn2<-1.e-12) 
+!! change last digits results
+!         WHERE( MAX( rn2(2:jpi,2:jpj,2:jpkm1), rn2b(2:jpi,2:jpj,2:jpkm1) )  <= -1.e-12 )
+!            avt(2:jpi,2:jpj,2:jpkm1) = rn_evd * wmask(2:jpi,2:jpj,2:jpkm1)
+!         END WHERE
+
          DO jk = 1, jpkm1
-!!!         WHERE( rn2(:,:,jk) <= -1.e-12 ) avt(:,:,jk) = tmask(:,:,jk) * avevd   ! agissant sur T SEUL! 
-            DO jj = 1, jpj             ! loop over the whole domain (no lbc_lnk call)
-               DO ji = 1, jpi
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1
                   IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 )   &
-                     avt(ji,jj,jk) = rn_evd * tmask(ji,jj,jk)
+                     avt(ji,jj,jk) = rn_evd * wmask(ji,jj,jk)
                END DO
             END DO
@@ -110,10 +115,8 @@
          !
       END SELECT 
-
+      !
       zavt_evd(:,:,:) = avt(:,:,:) - zavt_evd(:,:,:)   ! change in avt due to evd
       CALL iom_put( "avt_evd", zavt_evd )              ! output this change
       IF( l_trdtra ) CALL trd_tra( kt, 'TRA', jp_tem, jptra_evd, zavt_evd )
-      !
-      CALL wrk_dealloc( jpi,jpj,jpk,   zavt_evd, zavm_evd ) 
       !
       IF( nn_timing == 1 )  CALL timing_stop('zdf_evd')
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 7989)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 7990)
@@ -5,6 +5,7 @@
    !!                 turbulent closure parameterization
    !!======================================================================
-   !! History :   3.0  !  2009-09  (G. Reffray)  Original code
-   !!             3.3  !  2010-10  (C. Bricaud)  Add in the reference
+   !! History :  3.0  !  2009-09  (G. Reffray)  Original code
+   !!            3.3  !  2010-10  (C. Bricaud)  Add in the reference
+   !!            4.0  !  2017-04  (G. Madec)  remove CPP keys & avm at t-point only 
    !!----------------------------------------------------------------------
 
@@ -36,10 +37,10 @@
    PRIVATE
 
-   PUBLIC   zdf_gls        ! routine called in step module
-   PUBLIC   zdf_gls_init   ! routine called in zdfphy module
-   PUBLIC   gls_rst        ! routine called in zdfphy module
+   PUBLIC   zdf_gls        ! called in zdfphy
+   PUBLIC   zdf_gls_init   ! called in zdfphy
+   PUBLIC   gls_rst        ! called in zdfphy
 
    !
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   mxln    !: now mixing length
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hmxn    !: now mixing length
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zwall   !: wall function
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ustars2 !: Squared surface velocity scale at T-points
@@ -113,5 +114,5 @@
       !!                ***  FUNCTION zdf_gls_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( mxln(jpi,jpj,jpk), zwall(jpi,jpj,jpk) ,     &
+      ALLOCATE( hmxn(jpi,jpj,jpk), zwall(jpi,jpj,jpk) ,     &
          &      ustars2(jpi,jpj) , ustarb2(jpi,jpj)   , STAT= zdf_gls_alloc )
          !
@@ -148,21 +149,20 @@
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_elem_b    ! element of the second matrix diagonal
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_elem_c    ! element of the third  matrix diagonal
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z3du, z3dv  ! u- and v-shears
       !!--------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('zdf_gls')
+      IF( nn_timing == 1 )   CALL timing_start('zdf_gls')
       !
       CALL wrk_alloc( jpi,jpj,       zdep, zkar, zflxs, zhsro )
-      CALL wrk_alloc( jpi,jpj,jpk,   eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi  )
-      
+      CALL wrk_alloc( jpi,jpj,jpk,   eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi )
+      CALL wrk_alloc( jpi,jpj,jpk,   z3du, z3dv )
+
       ! Preliminary computing
 
       ustars2 = 0._wp   ;   ustarb2 = 0._wp   ;   psi  = 0._wp   ;   zwall_psi = 0._wp
 
-      IF( kt /= nit000 ) THEN   ! restore before value to compute tke
-         avt (:,:,:) = avt_k (:,:,:)
-         avm (:,:,:) = avm_k (:,:,:)
-         avmu(:,:,:) = avmu_k(:,:,:)
-         avmv(:,:,:) = avmv_k(:,:,:) 
-      ENDIF
+      ! restore before values computed GLS alone
+      avt(:,:,:) = avt_k(:,:,:)
+      avm(:,:,:) = avm_k(:,:,:)
 
       ! Compute surface and bottom friction at T-points
@@ -183,39 +183,37 @@
       END DO    
 
-      ! Set surface roughness length
-      SELECT CASE ( nn_z0_met )
-      !
-      CASE ( 0 )             ! Constant roughness          
+      SELECT CASE ( nn_z0_met )      !==  Set surface roughness length  ==!
+      CASE ( 0 )                          ! Constant roughness          
          zhsro(:,:) = rn_hsro
       CASE ( 1 )             ! Standard Charnock formula
          zhsro(:,:) = MAX(rsbc_zs1 * ustars2(:,:), rn_hsro)
       CASE ( 2 )             ! Roughness formulae according to Rascle et al., Ocean Modelling (2008)
+!!gm         zcof = 2._wp * 0.6_wp / 28._wp
+!!gm         zdep(:,:)  = 30._wp * TANH(  zcof/ SQRT( MAX(ustars2(:,:),rsmall) )  )      ! Wave age (eq. 10)
          zdep(:,:)  = 30.*TANH(2.*0.3/(28.*SQRT(MAX(ustars2(:,:),rsmall))))             ! Wave age (eq. 10)
          zhsro(:,:) = MAX(rsbc_zs2 * ustars2(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11)
       CASE ( 3 )             ! Roughness given by the wave model (coupled or read in file)
+!!gm  BUG missing a multiplicative coefficient....
          zhsro(:,:) = hsw(:,:)
       END SELECT
 
-      ! Compute shear and dissipation rate
-      DO jk = 2, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               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) )   &
-                  &                            / (  e3uw_n(ji,jj,jk) * e3uw_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) )   &
-                  &                            / (  e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) )
-               eps(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / mxln(ji,jj,jk)
+      DO jk = 2, jpkm1              !==  Compute shear and dissipation rate  ==!
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               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) ) / (  e3uw_n(ji,jj,jk) * e3uw_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) ) / (  e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) )
+                  !
+               eps(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / hmxn(ji,jj,jk)
             END DO
          END DO
       END DO
-      !
-      ! Lateral boundary conditions (avmu,avmv) (sign unchanged)
-      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )
 
       ! Save tke at before time step
       eb  (:,:,:) = en  (:,:,:)
-      mxlb(:,:,:) = mxln(:,:,:)
+      mxlb(:,:,:) = hmxn(:,:,:)
 
       IF( nn_clos == 0 ) THEN    ! Mellor-Yamada
@@ -223,5 +221,5 @@
             DO jj = 2, jpjm1 
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zup   = mxln(ji,jj,jk) * gdepw_n(ji,jj,mbkt(ji,jj)+1)
+                  zup   = hmxn(ji,jj,jk) * gdepw_n(ji,jj,mbkt(ji,jj)+1)
                   zdown = vkarmn * gdepw_n(ji,jj,jk) * ( -gdepw_n(ji,jj,jk) + gdepw_n(ji,jj,mbkt(ji,jj)+1) )
                   zcoef = ( zup / MAX( zdown, rsmall ) )
@@ -247,9 +245,9 @@
       DO jk = 2, jpkm1
          DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+            DO ji = 2, jpim1
                !
                ! shear prod. at w-point weightened by mask
-               shear(ji,jj,jk) =  ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1.e0 , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
-                  &             + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1.e0 , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )
+               shear(ji,jj,jk) =  ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1.e0 , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
+                  &             + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1.e0 , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )
                !
                ! stratif. destruction
@@ -278,19 +276,12 @@
                ! building the matrix
                zcof = rfact_tke * tmask(ji,jj,jk)
-               !
-               ! lower diagonal
-               z_elem_a(ji,jj,jk) = zcof * ( avm  (ji,jj,jk  ) + avm  (ji,jj,jk-1) )   &
-                  &                      / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk  ) )
-               !
-               ! upper diagonal
-               z_elem_c(ji,jj,jk) = zcof * ( avm  (ji,jj,jk+1) + avm  (ji,jj,jk  ) )   &
-                  &                      / ( e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk) )
-               !
-               ! diagonal
-               z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk)  &
-                  &                       + rdt * zdiss * tmask(ji,jj,jk) 
-               !
-               ! right hand side in en
-               en(ji,jj,jk) = en(ji,jj,jk) + rdt * zesh2 * tmask(ji,jj,jk)
+               !                                               ! lower diagonal
+               z_elem_a(ji,jj,jk) = zcof * ( avm(ji,jj,jk  ) + avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) )
+               !                                               ! upper diagonal
+               z_elem_c(ji,jj,jk) = zcof * ( avm(ji,jj,jk+1) + avm(ji,jj,jk  ) ) / ( e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk) )
+               !                                               ! diagonal
+               z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk)  + rdt * zdiss * wmask(ji,jj,jk) 
+               !                                               ! right hand side in en
+               en(ji,jj,jk) = en(ji,jj,jk) + rdt * zesh2 * wmask(ji,jj,jk)
             END DO
          END DO
@@ -300,6 +291,6 @@
       !
       ! Set surface condition on zwall_psi (1 at the bottom)
-      zwall_psi(:,:,1) = zwall_psi(:,:,2)
-      zwall_psi(:,:,jpk) = 1.
+      zwall_psi(:,:, 1 ) = zwall_psi(:,:,2)
+      zwall_psi(:,:,jpk) = 1._wp
       !
       ! Surface boundary condition on tke
@@ -353,5 +344,5 @@
       !
       CASE ( 0 )             ! Dirichlet 
-         !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = rn_lmin
+         !                      ! en(ibot) = u*^2 / Co2 and hmxn(ibot) = rn_lmin
          !                      ! Balance between the production and the dissipation terms
          DO jj = 2, jpjm1
@@ -503,16 +494,12 @@
                ! building the matrix
                zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk)
-               ! lower diagonal
-               z_elem_a(ji,jj,jk) = zcof * ( avm  (ji,jj,jk  ) + avm  (ji,jj,jk-1) )   &
-                  &                      / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk  ) )
-               ! upper diagonal
-               z_elem_c(ji,jj,jk) = zcof * ( avm  (ji,jj,jk+1) + avm  (ji,jj,jk  ) )   &
-                  &                      / ( e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk) )
-               ! diagonal
-               z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk)  &
-                  &                       + rdt * zdiss * tmask(ji,jj,jk)
-               !
-               ! right hand side in psi
-               psi(ji,jj,jk) = psi(ji,jj,jk) + rdt * zesh2 * tmask(ji,jj,jk)
+               !                                               ! lower diagonal
+               z_elem_a(ji,jj,jk) = zcof * ( avm(ji,jj,jk  ) + avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) )
+               !                                               ! upper diagonal
+               z_elem_c(ji,jj,jk) = zcof * ( avm(ji,jj,jk+1) + avm(ji,jj,jk  ) ) / ( e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk) )
+               !                                               ! diagonal
+               z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk)
+               !                                               ! right hand side in psi
+               psi(ji,jj,jk) = psi(ji,jj,jk) + rdt * zesh2 * wmask(ji,jj,jk)
             END DO
          END DO
@@ -565,5 +552,4 @@
       zflxs(:,:) = zdep(:,:) * zflxs(:,:)
       psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w_n(:,:,2)
-
       !   
       !
@@ -577,5 +563,5 @@
       !
       CASE ( 0 )             ! Dirichlet 
-         !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * rn_bfrz0
+         !                      ! en(ibot) = u*^2 / Co2 and hmxn(ibot) = vkarmn * rn_bfrz0
          !                      ! Balance between the production and the dissipation terms
          DO jj = 2, jpjm1
@@ -700,13 +686,13 @@
       ! Limit dissipation rate under stable stratification
       ! --------------------------------------------------
-      DO jk = 1, jpkm1 ! Note that this set boundary conditions on mxln at the same time
+      DO jk = 1, jpkm1 ! Note that this set boundary conditions on hmxn at the same time
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1    ! vector opt.
                ! limitation
                eps(ji,jj,jk)  = MAX( eps(ji,jj,jk), rn_epsmin )
-               mxln(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk)
+               hmxn(ji,jj,jk)  = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk)
                ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 
                zrn2 = MAX( rn2(ji,jj,jk), rsmall )
-               IF (ln_length_lim) mxln(ji,jj,jk) = MIN(  rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), mxln(ji,jj,jk) )
+               IF (ln_length_lim) hmxn(ji,jj,jk) = MIN(  rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), hmxn(ji,jj,jk) )
             END DO
          END DO
@@ -733,7 +719,7 @@
                   sm = ( rb1**(-1._wp/3._wp) + ( 18._wp*ra1*ra1 + 9._wp*ra1*ra2*(1._wp-rc2) )*sh*gh ) / (1._wp-9._wp*ra1*ra2*gh)
                   !
-                  ! Store stability function in avmu and avmv
-                  avmu(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)
-                  avmv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)
+                  ! Store stability function in z3du and z3dv
+                  z3du(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)
+                  z3dv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)
                END DO
             END DO
@@ -761,7 +747,7 @@
                   sh = (rs4 - rs5*gh + rs6*gm) / rcff
                   !
-                  ! Store stability function in avmu and avmv
-                  avmu(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)
-                  avmv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)
+                  ! Store stability function in z3du and z3dv
+                  z3du(ji,jj,jk) = rc_diff * sh * tmask(ji,jj,jk)
+                  z3dv(ji,jj,jk) = rc_diff * sm * tmask(ji,jj,jk)
                END DO
             END DO
@@ -773,9 +759,9 @@
       ! Lines below are useless if GOTM style Dirichlet conditions are used
 
-      avmv(:,:,1) = avmv(:,:,2)
+      z3dv(:,:,1) = z3dv(:,:,2)
 
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            avmv(ji,jj,mbkt(ji,jj)+1) = avmv(ji,jj,mbkt(ji,jj))
+            z3dv(ji,jj,mbkt(ji,jj)+1) = z3dv(ji,jj,mbkt(ji,jj))
          END DO
       END DO
@@ -786,44 +772,31 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zsqen         = SQRT( 2._wp * en(ji,jj,jk) ) * mxln(ji,jj,jk)
-               zav           = zsqen * avmu(ji,jj,jk)
-               avt(ji,jj,jk) = MAX( zav, avtb(jk) )*tmask(ji,jj,jk) ! apply mask for zdfmxl routine
-               zav           = zsqen * avmv(ji,jj,jk)
-               avm(ji,jj,jk) = MAX( zav, avmb(jk) ) ! Note that avm is not masked at the surface and the bottom
+               zsqen         = SQRT( 2._wp * en(ji,jj,jk) ) * hmxn(ji,jj,jk)
+               zav           = zsqen * z3du(ji,jj,jk)
+               avt(ji,jj,jk) = MAX( zav, avtb(jk) ) * wmask(ji,jj,jk) ! apply mask for zdfmxl routine
+               zav           = zsqen * z3dv(ji,jj,jk)
+               avm(ji,jj,jk) = MAX( zav, avmb(jk) )                   ! Note that avm is not masked at the surface and the bottom
             END DO
          END DO
       END DO
-      !
-      ! Lateral boundary conditions (sign unchanged)
       avt(:,:,1)  = 0._wp
+!!gm I'm sure this is needed to compute the shear term at next time-step
       CALL lbc_lnk( avm, 'W', 1. )   ;   CALL lbc_lnk( avt, 'W', 1. )
-
-      DO jk = 2, jpkm1            !* vertical eddy viscosity at u- and v-points
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * umask(ji,jj,jk)
-               avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * vmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-      avmu(:,:,1) = 0._wp             ;   avmv(:,:,1) = 0._wp                 ! set surface to zero
-      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )       ! Lateral boundary conditions
-
+      !
       IF(ln_ctl) THEN
-         CALL prt_ctl( tab3d_1=en  , clinfo1=' gls  - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk)
-         CALL prt_ctl( tab3d_1=avmu, clinfo1=' gls  - u: ', mask1=umask,                   &
-            &          tab3d_2=avmv, clinfo2=       ' v: ', mask2=vmask, ovlap=1, kdim=jpk )
+         CALL prt_ctl( tab3d_1=en , clinfo1=' gls  - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk)
+         CALL prt_ctl( tab3d_1=avm, clinfo1=' gls  - m: ', ovlap=1, kdim=jpk )
       ENDIF
       !
-      avt_k (:,:,:) = avt (:,:,:)
-      avm_k (:,:,:) = avm (:,:,:)
-      avmu_k(:,:,:) = avmu(:,:,:)
-      avmv_k(:,:,:) = avmv(:,:,:)
+      avt_k(:,:,:) = avt(:,:,:)      !==  store avt, avm values computed by GLS only  ==!
+      avm_k(:,:,:) = avm(:,:,:)
+      !
+      IF( lrst_oce )   CALL gls_rst( kt, 'WRITE' )     ! write the GLS info in the restart file
       !
       CALL wrk_dealloc( jpi,jpj,       zdep, zkar, zflxs, zhsro )
       CALL wrk_dealloc( jpi,jpj,jpk,   eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('zdf_gls')
-      !
+      CALL wrk_dealloc( jpi,jpj,jpk,   z3du, z3dv )
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('zdf_gls')
       !
    END SUBROUTINE zdf_gls
@@ -835,8 +808,7 @@
       !!                     
       !! ** Purpose :   Initialization of the vertical eddy diffivity and 
-      !!      viscosity when using a gls turbulent closure scheme
+      !!              viscosity computed using a GLS turbulent closure scheme
       !!
       !! ** Method  :   Read the namzdf_gls namelist and check the parameters
-      !!      called at the first timestep (nit000)
       !!
       !! ** input   :   Namlist namzdf_gls
@@ -892,5 +864,5 @@
       ENDIF
 
-      !                                !* allocate gls arrays
+      !                                !* allocate GLS arrays
       IF( zdf_gls_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_gls_init : unable to allocate arrays' )
 
@@ -1120,20 +1092,13 @@
       rsbc_psi1 = -0.5_wp * rdt * rc0**(rpp-2._wp*rmm) / rsc_psi
       rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking 
-
+      !
       rfact_tke = -0.5_wp / rsc_tke * rdt                                ! Cst used for the Diffusion term of tke
       rfact_psi = -0.5_wp / rsc_psi * rdt                                ! Cst used for the Diffusion term of tke
-
+      !
       !                                !* Wall proximity function
       zwall (:,:,:) = 1._wp * tmask(:,:,:)
 
-      !                                !* set vertical eddy coef. to the background value
-      DO jk = 1, jpk
-         avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
-         avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)
-         avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
-         avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
-      END DO
-      !                              
-      CALL gls_rst( nit000, 'READ' )   !* read or initialize all required files
+      !                                !* read or initialize all required files  
+      CALL gls_rst( nit000, 'READ' )      ! (en, avt_k, avm_k, hmxn)
       !
       IF( nn_timing == 1 )  CALL timing_stop('zdf_gls_init')
@@ -1152,9 +1117,9 @@
       !!                set to rn_emin or recomputed (nn_igls/=0)
       !!----------------------------------------------------------------------
-      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
-      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
+      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
+      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
       !
       INTEGER ::   jit, jk   ! dummy loop indices
-      INTEGER ::   id1, id2, id3, id4, id5, id6
+      INTEGER ::   id1, id2, id3, id4
       INTEGER ::   ji, jj, ikbu, ikbv
       REAL(wp)::   cbx, cby
@@ -1165,31 +1130,29 @@
          IF( ln_rstart ) THEN                   !* Read the restart file
             id1 = iom_varid( numror, 'en'   , ldstop = .FALSE. )
-            id2 = iom_varid( numror, 'avt'  , ldstop = .FALSE. )
-            id3 = iom_varid( numror, 'avm'  , ldstop = .FALSE. )
-            id4 = iom_varid( numror, 'avmu' , ldstop = .FALSE. )
-            id5 = iom_varid( numror, 'avmv' , ldstop = .FALSE. )
-            id6 = iom_varid( numror, 'mxln' , ldstop = .FALSE. )
+            id2 = iom_varid( numror, 'avt_k', ldstop = .FALSE. )
+            id3 = iom_varid( numror, 'avm_k', ldstop = .FALSE. )
+            id4 = iom_varid( numror, 'hmxn' , ldstop = .FALSE. )
             !
-            IF( MIN( id1, id2, id3, id4, id5, id6 ) > 0 ) THEN        ! all required arrays exist
-               CALL iom_get( numror, jpdom_autoglo, 'en'    , en     )
-               CALL iom_get( numror, jpdom_autoglo, 'avt'   , avt    )
-               CALL iom_get( numror, jpdom_autoglo, 'avm'   , avm    )
-               CALL iom_get( numror, jpdom_autoglo, 'avmu'  , avmu   )
-               CALL iom_get( numror, jpdom_autoglo, 'avmv'  , avmv   )
-               CALL iom_get( numror, jpdom_autoglo, 'mxln'  , mxln   )
+            IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN        ! all required arrays exist
+               CALL iom_get( numror, jpdom_autoglo, 'en'   , en    )
+               CALL iom_get( numror, jpdom_autoglo, 'avt_k', avt_k )
+               CALL iom_get( numror, jpdom_autoglo, 'avm_k', avm_k )
+               CALL iom_get( numror, jpdom_autoglo, 'hmxn' , hmxn  )
             ELSE                        
-               IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and mxln computed by iterative loop'
-               en  (:,:,:) = rn_emin
-               mxln(:,:,:) = 0.05        
-               avt_k (:,:,:) = avt (:,:,:)
-               avm_k (:,:,:) = avm (:,:,:)
-               avmu_k(:,:,:) = avmu(:,:,:)
-               avmv_k(:,:,:) = avmv(:,:,:)
+               IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and hmxn computed by iterative loop'
+               en   (:,:,:) = rn_emin
+               hmxn (:,:,:) = 0.05_wp
+               avt_k(:,:,:) = avt(:,:,:)
+               avm_k(:,:,:) = avm(:,:,:)
                DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_gls( jit )   ;   END DO
             ENDIF
          ELSE                                   !* Start from rest
-            IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and mxln by background values'
+            IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and hmxn by background values'
             en  (:,:,:) = rn_emin
-            mxln(:,:,:) = 0.05       
+            hmxn(:,:,:) = 0.05_wp
+            DO jk = 1, jpk
+               avt_k(:,:,jk) = avtb(jk) * wmask(:,:,jk)
+               avm_k(:,:,jk) = avmb(jk) * wmask(:,:,jk)
+            END DO
          ENDIF
          !
@@ -1197,10 +1160,8 @@
          !                                   ! -------------------
          IF(lwp) WRITE(numout,*) '---- gls-rst ----'
-         CALL iom_rstput( kt, nitrst, numrow, 'en'   , en     ) 
-         CALL iom_rstput( kt, nitrst, numrow, 'avt'  , avt_k  )
-         CALL iom_rstput( kt, nitrst, numrow, 'avm'  , avm_k  )
-         CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k ) 
-         CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv_k )
-         CALL iom_rstput( kt, nitrst, numrow, 'mxln' , mxln   )
+         CALL iom_rstput( kt, nitrst, numrow, 'en'   , en    ) 
+         CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k )
+         CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k )
+         CALL iom_rstput( kt, nitrst, numrow, 'hmxn' , hmxn  )
          !
       ENDIF
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfiwm.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfiwm.F90	(revision 7990)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfiwm.F90	(revision 7990)
@@ -0,0 +1,482 @@
+MODULE zdfiwm
+   !!========================================================================
+   !!                       ***  MODULE  zdfiwm  ***
+   !! Ocean physics: Internal gravity wave-driven vertical mixing
+   !!========================================================================
+   !! History :  1.0  !  2004-04  (L. Bessieres, G. Madec)  Original code
+   !!             -   !  2006-08  (A. Koch-Larrouy)  Indonesian strait
+   !!            3.3  !  2010-10  (C. Ethe, G. Madec)  reorganisation of initialisation phase
+   !!            3.6  !  2016-03  (C. de Lavergne)  New param: internal wave-driven mixing 
+   !!            4.0  !  2017-04  (G. Madec)  renamed module, remove the old param. and the CPP keys
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   zdf_iwm       : global     momentum & tracer Kz with wave induced Kz
+   !!   zdf_iwm_init  : global     momentum & tracer Kz with wave induced Kz
+   !!----------------------------------------------------------------------
+   USE oce            ! ocean dynamics and tracers variables
+   USE dom_oce        ! ocean space and time domain variables
+   USE zdf_oce        ! ocean vertical physics variables
+   USE zdfddm         ! ocean vertical physics: double diffusive mixing
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+   USE eosbn2         ! ocean equation of state
+   USE phycst         ! physical constants
+   USE prtctl         ! Print control
+   USE in_out_manager ! I/O manager
+   USE iom            ! I/O Manager
+   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)  
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   zdf_iwm         ! called in step module 
+   PUBLIC   zdf_iwm_init    ! called in nemogcm module 
+   PUBLIC   zdf_iwm_alloc   ! called in nemogcm module
+
+   !                       !!* Namelist  namzdf_iwm : internal wave-driven mixing *
+   INTEGER  ::  nn_zpyc     ! pycnocline-intensified mixing energy proportional to N (=1) or N^2 (=2)
+   LOGICAL  ::  ln_mevar    ! variable (=T) or constant (=F) mixing efficiency
+   LOGICAL  ::  ln_tsdiff   ! account for differential T/S wave-driven mixing (=T) or not (=F)
+
+   REAL(wp) ::  r1_6 = 1._wp / 6._wp
+
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ebot_iwm     ! power available from high-mode wave breaking (W/m2)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   epyc_iwm     ! power available from low-mode, pycnocline-intensified wave breaking (W/m2)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ecri_iwm     ! power available from low-mode, critical slope wave breaking (W/m2)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hbot_iwm     ! WKB decay scale for high-mode energy dissipation (m)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hcri_iwm     ! decay scale for low-mode critical slope dissipation (m)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   emix_iwm     ! local energy density available for mixing (W/kg)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   bflx_iwm     ! buoyancy flux Kz * N^2 (W/kg)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   pcmap_iwm    ! vertically integrated buoyancy flux (W/m2)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zav_ratio    ! S/T diffusivity ratio (only for ln_tsdiff=T)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zav_wave     ! Internal wave-induced diffusivity
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , NEMO Consortium (2016)
+   !! $Id$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   INTEGER FUNCTION zdf_iwm_alloc()
+      !!----------------------------------------------------------------------
+      !!                ***  FUNCTION zdf_iwm_alloc  ***
+      !!----------------------------------------------------------------------
+      ALLOCATE(     ebot_iwm(jpi,jpj),  epyc_iwm(jpi,jpj),  ecri_iwm(jpi,jpj)    ,   &
+      &             hbot_iwm(jpi,jpj),  hcri_iwm(jpi,jpj),  emix_iwm(jpi,jpj,jpk),   &
+      &         bflx_iwm(jpi,jpj,jpk), pcmap_iwm(jpi,jpj), zav_ratio(jpi,jpj,jpk),   & 
+      &         zav_wave(jpi,jpj,jpk), STAT=zdf_iwm_alloc     )
+      !
+      IF( lk_mpp             )   CALL mpp_sum ( zdf_iwm_alloc )
+      IF( zdf_iwm_alloc /= 0 )   CALL ctl_warn('zdf_iwm_alloc: failed to allocate arrays')
+   END FUNCTION zdf_iwm_alloc
+
+
+   SUBROUTINE zdf_iwm( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE zdf_iwm  ***
+      !!                   
+      !! ** Purpose :   add to the vertical mixing coefficients the effect of
+      !!              breaking internal waves.
+      !!
+      !! ** Method  : - internal wave-driven vertical mixing is given by:
+      !!                  Kz_wave = min(  100 cm2/s, f(  Reb = emix_iwm /( Nu * N^2 )  )
+      !!              where emix_iwm is the 3D space distribution of the wave-breaking 
+      !!              energy and Nu the molecular kinematic viscosity.
+      !!              The function f(Reb) is linear (constant mixing efficiency)
+      !!              if the namelist parameter ln_mevar = F and nonlinear if ln_mevar = T.
+      !!
+      !!              - Compute emix_iwm, the 3D power density that allows to compute
+      !!              Reb and therefrom the wave-induced vertical diffusivity.
+      !!              This is divided into three components:
+      !!                 1. Bottom-intensified low-mode dissipation at critical slopes
+      !!                     emix_iwm(z) = ( ecri_iwm / rau0 ) * EXP( -(H-z)/hcri_iwm )
+      !!                                   / ( 1. - EXP( - H/hcri_iwm ) ) * hcri_iwm
+      !!              where hcri_iwm is the characteristic length scale of the bottom 
+      !!              intensification, ecri_iwm a map of available power, and H the ocean depth.
+      !!                 2. Pycnocline-intensified low-mode dissipation
+      !!                     emix_iwm(z) = ( epyc_iwm / rau0 ) * ( sqrt(rn2(z))^nn_zpyc )
+      !!                                   / SUM( sqrt(rn2(z))^nn_zpyc * e3w(z) )
+      !!              where epyc_iwm is a map of available power, and nn_zpyc
+      !!              is the chosen stratification-dependence of the internal wave
+      !!              energy dissipation.
+      !!                 3. WKB-height dependent high mode dissipation
+      !!                     emix_iwm(z) = ( ebot_iwm / rau0 ) * rn2(z) * EXP(-z_wkb(z)/hbot_iwm)
+      !!                                   / SUM( rn2(z) * EXP(-z_wkb(z)/hbot_iwm) * e3w(z) )
+      !!              where hbot_iwm is the characteristic length scale of the WKB bottom 
+      !!              intensification, ebot_iwm is a map of available power, and z_wkb is the
+      !!              WKB-stretched height above bottom defined as
+      !!                    z_wkb(z) = H * SUM( sqrt(rn2(z'>=z)) * e3w(z'>=z) )
+      !!                                 / SUM( sqrt(rn2(z'))    * e3w(z')    )
+      !!
+      !!              - update the model vertical eddy viscosity and diffusivity: 
+      !!                     avt  = avt  +    av_wave
+      !!                     avm  = avm  +    av_wave
+      !!
+      !!              - if namelist parameter ln_tsdiff = T, account for differential mixing:
+      !!                     avs  = avt  +    av_wave * diffusivity_ratio(Reb)
+      !!
+      !! ** Action  : - Define emix_iwm used to compute internal wave-induced mixing
+      !!              - avt, avs, avm, increased by internal wave-driven mixing    
+      !!
+      !! References :  de Lavergne et al. 2015, JPO; 2016, in prep.
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step 
+      !
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   ztpc         ! scalar workspace
+      REAL(wp), DIMENSION(:,:)  , POINTER ::  zfact     ! Used for vertical structure
+      REAL(wp), DIMENSION(:,:)  , POINTER ::  zhdep     ! Ocean depth
+      REAL(wp), DIMENSION(:,:,:), POINTER ::  zwkb      ! WKB-stretched height above bottom
+      REAL(wp), DIMENSION(:,:,:), POINTER ::  zweight   ! Weight for high mode vertical distribution
+      REAL(wp), DIMENSION(:,:,:), POINTER ::  znu_t     ! Molecular kinematic viscosity (T grid)
+      REAL(wp), DIMENSION(:,:,:), POINTER ::  znu_w     ! Molecular kinematic viscosity (W grid)
+      REAL(wp), DIMENSION(:,:,:), POINTER ::  zReb      ! Turbulence intensity parameter
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )   CALL timing_start('zdf_iwm')
+      !
+      CALL wrk_alloc( jpi,jpj,       zfact, zhdep )
+      CALL wrk_alloc( jpi,jpj,jpk,   zwkb, zweight, znu_t, znu_w, zReb )
+
+      !                          ! ----------------------------- !
+      !                          !  Internal wave-driven mixing  !  (compute zav_wave)
+      !                          ! ----------------------------- !
+      !                             
+      !                        !* Critical slope mixing: distribute energy over the time-varying ocean depth,
+      !                                                 using an exponential decay from the seafloor.
+      DO jj = 1, jpj                ! part independent of the level
+         DO ji = 1, jpi
+            zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1)       ! depth of the ocean
+            zfact(ji,jj) = rau0 * (  1._wp - EXP( -zhdep(ji,jj) / hcri_iwm(ji,jj) )  )
+            IF( zfact(ji,jj) /= 0._wp )   zfact(ji,jj) = ecri_iwm(ji,jj) / zfact(ji,jj)
+         END DO
+      END DO
+
+      DO jk = 2, jpkm1              ! complete with the level-dependent part
+         emix_iwm(:,:,jk) = zfact(:,:) * (  EXP( ( gde3w_n(:,:,jk  ) - zhdep(:,:) ) / hcri_iwm(:,:) )                      &
+            &                             - EXP( ( gde3w_n(:,:,jk-1) - zhdep(:,:) ) / hcri_iwm(:,:) )  ) * wmask(:,:,jk)   &
+!!gm delta(gde3w_n) = e3t_n  !!  Please verify the grid-point position w versus t-point
+            &                          / ( gde3w_n(:,:,jk) - gde3w_n(:,:,jk-1) )
+      END DO
+
+      !                        !* Pycnocline-intensified mixing: distribute energy over the time-varying 
+      !                        !* ocean depth as proportional to sqrt(rn2)^nn_zpyc
+      !
+      SELECT CASE ( nn_zpyc )
+      !
+      CASE ( 1 )               ! Dissipation scales as N (recommended)
+         !
+         zfact(:,:) = 0._wp
+         DO jk = 2, jpkm1              ! part independent of the level
+            zfact(:,:) = zfact(:,:) + e3w_n(:,:,jk) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk)
+         END DO
+         !
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rau0 * zfact(ji,jj) )
+            END DO
+         END DO
+         !
+         DO jk = 2, jpkm1              ! complete with the level-dependent part
+            emix_iwm(:,:,jk) = emix_iwm(:,:,jk) + zfact(:,:) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk)
+         END DO
+         !
+      CASE ( 2 )               ! Dissipation scales as N^2
+         !
+         zfact(:,:) = 0._wp
+         DO jk = 2, jpkm1              ! part independent of the level
+            zfact(:,:) = zfact(:,:) + e3w_n(:,:,jk) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk)
+         END DO
+         !
+         DO jj= 1, jpj
+            DO ji = 1, jpi
+               IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_iwm(ji,jj) / ( rau0 * zfact(ji,jj) )
+            END DO
+         END DO
+         !
+         DO jk = 2, jpkm1              ! complete with the level-dependent part
+            emix_iwm(:,:,jk) = emix_iwm(:,:,jk) + zfact(:,:) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk)
+         END DO
+         !
+      END SELECT
+
+      !                        !* WKB-height dependent mixing: distribute energy over the time-varying 
+      !                        !* ocean depth as proportional to rn2 * exp(-z_wkb/rn_hbot)
+      !
+      zwkb(:,:,:) = 0._wp
+      zfact(:,:) = 0._wp
+      DO jk = 2, jpkm1
+         zfact(:,:) = zfact(:,:) + e3w_n(:,:,jk) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk)
+         zwkb(:,:,jk) = zfact(:,:)
+      END DO
+      !
+      DO jk = 2, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( zfact(ji,jj) /= 0 )   zwkb(ji,jj,jk) = zhdep(ji,jj) * ( zfact(ji,jj) - zwkb(ji,jj,jk) )   &
+                  &                                     * tmask(ji,jj,jk) / zfact(ji,jj)
+            END DO
+         END DO
+      END DO
+      zwkb(:,:,1) = zhdep(:,:) * tmask(:,:,1)
+      !
+      zweight(:,:,:) = 0._wp
+      DO jk = 2, jpkm1
+         zweight(:,:,jk) = MAX( 0._wp, rn2(:,:,jk) ) * hbot_iwm(:,:) * wmask(:,:,jk)                    &
+            &   * (  EXP( -zwkb(:,:,jk) / hbot_iwm(:,:) ) - EXP( -zwkb(:,:,jk-1) / hbot_iwm(:,:) )  )
+      END DO
+      !
+      zfact(:,:) = 0._wp
+      DO jk = 2, jpkm1              ! part independent of the level
+         zfact(:,:) = zfact(:,:) + zweight(:,:,jk)
+      END DO
+      !
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ebot_iwm(ji,jj) / ( rau0 * zfact(ji,jj) )
+         END DO
+      END DO
+      !
+      DO jk = 2, jpkm1              ! complete with the level-dependent part
+         emix_iwm(:,:,jk) = emix_iwm(:,:,jk) + zweight(:,:,jk) * zfact(:,:) * wmask(:,:,jk)   &
+            &                                / ( gde3w_n(:,:,jk) - gde3w_n(:,:,jk-1) )
+      END DO
+      !
+      ! Calculate molecular kinematic viscosity
+      znu_t(:,:,:) = 1.e-4_wp * (  17.91_wp - 0.53810_wp * tsn(:,:,:,jp_tem) + 0.00694_wp * tsn(:,:,:,jp_tem) * tsn(:,:,:,jp_tem)  &
+         &                                  + 0.02305_wp * tsn(:,:,:,jp_sal)  ) * tmask(:,:,:) * r1_rau0
+      DO jk = 2, jpkm1
+         znu_w(:,:,jk) = 0.5_wp * ( znu_t(:,:,jk-1) + znu_t(:,:,jk) ) * wmask(:,:,jk)
+      END DO
+      !
+      ! Calculate turbulence intensity parameter Reb
+      DO jk = 2, jpkm1
+         zReb(:,:,jk) = emix_iwm(:,:,jk) / MAX( 1.e-20_wp, znu_w(:,:,jk) * rn2(:,:,jk) )
+      END DO
+      !
+      ! Define internal wave-induced diffusivity
+      DO jk = 2, jpkm1
+         zav_wave(:,:,jk) = znu_w(:,:,jk) * zReb(:,:,jk) * r1_6   ! This corresponds to a constant mixing efficiency of 1/6
+      END DO
+      !
+      IF( ln_mevar ) THEN              ! Variable mixing efficiency case : modify zav_wave in the
+         DO jk = 2, jpkm1              ! energetic (Reb > 480) and buoyancy-controlled (Reb <10.224 ) regimes
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  IF( zReb(ji,jj,jk) > 480.00_wp ) THEN
+                     zav_wave(ji,jj,jk) = 3.6515_wp * znu_w(ji,jj,jk) * SQRT( zReb(ji,jj,jk) )
+                  ELSEIF( zReb(ji,jj,jk) < 10.224_wp ) THEN
+                     zav_wave(ji,jj,jk) = 0.052125_wp * znu_w(ji,jj,jk) * zReb(ji,jj,jk) * SQRT( zReb(ji,jj,jk) )
+                  ENDIF
+               END DO
+            END DO
+         END DO
+      ENDIF
+      !
+      DO jk = 2, jpkm1                 ! Bound diffusivity by molecular value and 100 cm2/s
+         zav_wave(:,:,jk) = MIN(  MAX( 1.4e-7_wp, zav_wave(:,:,jk) ), 1.e-2_wp  ) * wmask(:,:,jk)
+      END DO
+      !
+      IF( kt == nit000 ) THEN        !* Control print at first time-step: diagnose the energy consumed by zav_wave
+         ztpc = 0._wp
+!!gm used of glosum 3D....
+         DO jk = 2, jpkm1
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ztpc = ztpc + e3w_n(ji,jj,jk) * e1e2t(ji,jj)   &
+                     &         * MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj)
+               END DO
+            END DO
+         END DO
+         IF( lk_mpp )   CALL mpp_sum( ztpc )
+         ztpc = rau0 * ztpc ! Global integral of rauo * Kz * N^2 = power contributing to mixing 
+         !
+         IF(lwp) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) 'zdf_iwm : Internal wave-driven mixing (iwm)'
+            WRITE(numout,*) '~~~~~~~ '
+            WRITE(numout,*)
+            WRITE(numout,*) '      Total power consumption by av_wave: ztpc =  ', ztpc * 1.e-12_wp, 'TW'
+         ENDIF
+      ENDIF
+
+      !                          ! ----------------------- !
+      !                          !   Update  mixing coefs  !                          
+      !                          ! ----------------------- !
+      !      
+      IF( ln_tsdiff ) THEN          !* Option for differential mixing of salinity and temperature
+         DO jk = 2, jpkm1              ! Calculate S/T diffusivity ratio as a function of Reb
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zav_ratio(ji,jj,jk) = ( 0.505_wp + 0.495_wp *                                                                  &
+                      &   TANH(    0.92_wp * (   LOG10(  MAX( 1.e-20_wp, zReb(ji,jj,jk) * 5._wp * r1_6 )  ) - 0.60_wp   )    )   &
+                      &                 ) * wmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         CALL iom_put( "av_ratio", zav_ratio )
+         DO jk = 2, jpkm1           !* update momentum & tracer diffusivity with wave-driven mixing
+            avs(:,:,jk) = avs(:,:,jk) + zav_wave(:,:,jk) * zav_ratio(:,:,jk)
+            avt(:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk)
+            avm(:,:,jk) = avm(:,:,jk) + zav_wave(:,:,jk)
+         END DO
+         !
+      ELSE                          !* update momentum & tracer diffusivity with wave-driven mixing
+         DO jk = 2, jpkm1
+            avs(:,:,jk) = avs(:,:,jk) + zav_wave(:,:,jk)
+            avt(:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk)
+            avm(:,:,jk) = avm(:,:,jk) + zav_wave(:,:,jk)
+         END DO
+      ENDIF
+
+      !                             !* output internal wave-driven mixing coefficient
+      CALL iom_put( "av_wave", zav_wave )
+                                    !* output useful diagnostics: N^2, Kz * N^2 (bflx_iwm), 
+                                    !  vertical integral of rau0 * Kz * N^2 (pcmap_iwm), energy density (emix_iwm)
+      IF( iom_use("bflx_iwm") .OR. iom_use("pcmap_iwm") ) THEN
+         bflx_iwm(:,:,:) = MAX( 0._wp, rn2(:,:,:) ) * zav_wave(:,:,:)
+         pcmap_iwm(:,:) = 0._wp
+         DO jk = 2, jpkm1
+            pcmap_iwm(:,:) = pcmap_iwm(:,:) + e3w_n(:,:,jk) * bflx_iwm(:,:,jk) * wmask(:,:,jk)
+         END DO
+         pcmap_iwm(:,:) = rau0 * pcmap_iwm(:,:)
+         CALL iom_put( "bflx_iwm", bflx_iwm )
+         CALL iom_put( "pcmap_iwm", pcmap_iwm )
+      ENDIF
+      CALL iom_put( "bn2", rn2 )
+      CALL iom_put( "emix_iwm", emix_iwm )
+      
+      CALL wrk_dealloc( jpi,jpj,       zfact, zhdep )
+      CALL wrk_dealloc( jpi,jpj,jpk,   zwkb, zweight, znu_t, znu_w, zReb )
+
+      IF(ln_ctl)   CALL prt_ctl(tab3d_1=zav_wave , clinfo1=' iwm - av_wave: ', tab3d_2=avt, clinfo2=' avt: ', ovlap=1, kdim=jpk)
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('zdf_iwm')
+      !
+   END SUBROUTINE zdf_iwm
+
+
+   SUBROUTINE zdf_iwm_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE zdf_iwm_init  ***
+      !!                     
+      !! ** Purpose :   Initialization of the wave-driven vertical mixing, reading
+      !!              of input power maps and decay length scales in netcdf files.
+      !!
+      !! ** Method  : - Read the namzdf_iwm namelist and check the parameters
+      !!
+      !!              - Read the input data in NetCDF files :
+      !!              power available from high-mode wave breaking (mixing_power_bot.nc)
+      !!              power available from pycnocline-intensified wave-breaking (mixing_power_pyc.nc)
+      !!              power available from critical slope wave-breaking (mixing_power_cri.nc)
+      !!              WKB decay scale for high-mode wave-breaking (decay_scale_bot.nc)
+      !!              decay scale for critical slope wave-breaking (decay_scale_cri.nc)
+      !!
+      !! ** input   : - Namlist namzdf_iwm
+      !!              - NetCDF files : mixing_power_bot.nc, mixing_power_pyc.nc, mixing_power_cri.nc,
+      !!              decay_scale_bot.nc decay_scale_cri.nc
+      !!
+      !! ** Action  : - Increase by 1 the nstop flag is setting problem encounter
+      !!              - Define ebot_iwm, epyc_iwm, ecri_iwm, hbot_iwm, hcri_iwm
+      !!
+      !! References : de Lavergne et al. JPO, 2015 ; de Lavergne PhD 2016
+      !!              de Lavergne et al. in prep., 2017
+      !!----------------------------------------------------------------------
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      INTEGER  ::   inum         ! local integer
+      INTEGER  ::   ios
+      REAL(wp) ::   zbot, zpyc, zcri   ! local scalars
+      !!
+      NAMELIST/namzdf_iwm_new/ nn_zpyc, ln_mevar, ln_tsdiff
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('zdf_iwm_init')
+      !
+      REWIND( numnam_ref )              ! Namelist namzdf_iwm in reference namelist : Wave-driven mixing
+      READ  ( numnam_ref, namzdf_iwm_new, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_iwm in reference namelist', lwp )
+      !
+      REWIND( numnam_cfg )              ! Namelist namzdf_iwm in configuration namelist : Wave-driven mixing
+      READ  ( numnam_cfg, namzdf_iwm_new, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_iwm in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, namzdf_iwm_new )
+      !
+      IF(lwp) THEN                  ! Control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'zdf_iwm_init : internal wave-driven mixing'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist namzdf_iwm_new : set wave-driven mixing parameters'
+         WRITE(numout,*) '      Pycnocline-intensified diss. scales as N (=1) or N^2 (=2) = ', nn_zpyc
+         WRITE(numout,*) '      Variable (T) or constant (F) mixing efficiency            = ', ln_mevar
+         WRITE(numout,*) '      Differential internal wave-driven mixing (T) or not (F)   = ', ln_tsdiff
+      ENDIF
+      
+      ! The new wave-driven mixing parameterization elevates avt and avm in the interior, and
+      ! ensures that avt remains larger than its molecular value (=1.4e-7). Therefore, avtb should 
+      ! be set here to a very small value, and avmb to its (uniform) molecular value (=1.4e-6).
+      avmb(:) = 1.4e-6_wp        ! viscous molecular value
+      avtb(:) = 1.e-10_wp        ! very small diffusive minimum (background avt is specified in zdf_iwm)    
+      avtb_2d(:,:) = 1.e0_wp     ! uniform 
+      IF(lwp) THEN                  ! Control print
+         WRITE(numout,*)
+         WRITE(numout,*) '   Force the background value applied to avm & avt in TKE to be everywhere ',   &
+            &               'the viscous molecular value & a very small diffusive value, resp.'
+      ENDIF
+            
+      !                             ! allocate iwm arrays
+      IF( zdf_iwm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_iwm_init : unable to allocate iwm arrays' )
+      !
+      !                             ! read necessary fields
+      CALL iom_open('mixing_power_bot',inum)       ! energy flux for high-mode wave breaking [W/m2]
+      CALL iom_get  (inum, jpdom_data, 'field', ebot_iwm, 1 ) 
+      CALL iom_close(inum)
+      !
+      CALL iom_open('mixing_power_pyc',inum)       ! energy flux for pynocline-intensified wave breaking [W/m2]
+      CALL iom_get  (inum, jpdom_data, 'field', epyc_iwm, 1 )
+      CALL iom_close(inum)
+      !
+      CALL iom_open('mixing_power_cri',inum)       ! energy flux for critical slope wave breaking [W/m2]
+      CALL iom_get  (inum, jpdom_data, 'field', ecri_iwm, 1 )
+      CALL iom_close(inum)
+      !
+      CALL iom_open('decay_scale_bot',inum)        ! spatially variable decay scale for high-mode wave breaking [m]
+      CALL iom_get  (inum, jpdom_data, 'field', hbot_iwm, 1 )
+      CALL iom_close(inum)
+      !
+      CALL iom_open('decay_scale_cri',inum)        ! spatially variable decay scale for critical slope wave breaking [m]
+      CALL iom_get  (inum, jpdom_data, 'field', hcri_iwm, 1 )
+      CALL iom_close(inum)
+
+      ebot_iwm(:,:) = ebot_iwm(:,:) * ssmask(:,:)
+      epyc_iwm(:,:) = epyc_iwm(:,:) * ssmask(:,:)
+      ecri_iwm(:,:) = ecri_iwm(:,:) * ssmask(:,:)
+
+      ! Set once for all to zero the first and last vertical levels of appropriate variables
+      emix_iwm (:,:, 1 ) = 0._wp
+      emix_iwm (:,:,jpk) = 0._wp
+      zav_ratio(:,:, 1 ) = 0._wp
+      zav_ratio(:,:,jpk) = 0._wp
+      zav_wave (:,:, 1 ) = 0._wp
+      zav_wave (:,:,jpk) = 0._wp
+
+      zbot = glob_sum( e1e2t(:,:) * ebot_iwm(:,:) )
+      zpyc = glob_sum( e1e2t(:,:) * epyc_iwm(:,:) )
+      zcri = glob_sum( e1e2t(:,:) * ecri_iwm(:,:) )
+      IF(lwp) THEN
+         WRITE(numout,*) '      High-mode wave-breaking energy:             ', zbot * 1.e-12_wp, 'TW'
+         WRITE(numout,*) '      Pycnocline-intensifed wave-breaking energy: ', zpyc * 1.e-12_wp, 'TW'
+         WRITE(numout,*) '      Critical slope wave-breaking energy:        ', zcri * 1.e-12_wp, 'TW'
+      ENDIF
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('zdf_iwm_init')
+      !
+   END SUBROUTINE zdf_iwm_init
+
+   !!======================================================================
+END MODULE zdfiwm
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfphy.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfphy.F90	(revision 7989)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfphy.F90	(revision 7990)
@@ -2,5 +2,5 @@
    !!======================================================================
    !!                      ***  MODULE  zdfphy  ***
-   !! Ocean physics :   manager of vertical mixing parametrizations
+   !! Vertical ocean physics :   manager of all vertical physics packages
    !!======================================================================
    !! History :  4.0  !  2017-04  (G. Madec)  original code
@@ -8,24 +8,26 @@
 
    !!----------------------------------------------------------------------
-   !!   zdf_phy_init  : initialization of all vertical physics pakages
+   !!   zdf_phy_init  : initialization of all vertical physics packages
    !!   zdf_phy       : upadate at each time-step the vertical mixing coeff. 
    !!----------------------------------------------------------------------
-   USE par_oce        ! mesh and scale factors
-   USE zdf_oce        ! TKE vertical mixing          
+   USE par_oce        ! ocean parameters
+   USE zdf_oce        ! vertical physics: shared variables         
+   USE zdfbfr         ! vertical physics: bottom friction
+   USE zdfric         ! vertical physics: Richardson vertical mixing   
+   USE zdftke         ! vertical physics: TKE vertical mixing
+   USE zdfgls         ! vertical physics: GLS vertical mixing
+   USE zdfddm         ! vertical physics: double diffusion mixing      
+   USE zdfevd         ! vertical physics: convection via enhanced vertical diffusion  
+   USE zdfiwm         ! vertical physics: internal wave-induced mixing  
+   USE zdfswm         ! vertical physics: surface  wave-induced mixing
+   USE zdfmxl         ! vertical physics: mixed layer
+   USE tranpc         ! convection: non penetrative adjustment
+   USE trc_oce        ! variables shared between passive tracer & ocean           
    USE sbc_oce        ! surface module (only for nn_isf in the option compatibility test)
-   USE zdfbfr         ! bottom friction
-   USE zdftke         ! TKE vertical mixing
-   USE zdfgls         ! GLS vertical mixing
-   USE zdfric         ! Richardson vertical mixing   
-   USE zdfddm         ! double diffusion mixing      
-   USE zdfevd         ! enhanced vertical diffusion
-   USE zdftmx         ! internal tide-induced mixing
-   USE zdfqiao          !Qiao module wave induced mixing   (zdf_qiao routine)
-   USE zdfmxl           ! Mixed-layer depth                (zdf_mxl routine)
-   USE tranpc         ! convection: non penetrative adjustment
    USE sbcrnf         ! surface boundary condition: runoff variables
    !
    USE in_out_manager ! I/O manager
    USE iom            ! IOM library
+   USE lbclnk         ! lateral boundary conditions
    USE lib_mpp        ! distribued memory computing
 
@@ -33,10 +35,16 @@
    PRIVATE
 
-   PUBLIC   zdf_phy_init   ! routine called by nemogcm.F90
-   PUBLIC   zdf_phy        ! routine called by step.F90
-   
-   
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
+   PUBLIC   zdf_phy_init  ! called by nemogcm.F90
+   PUBLIC   zdf_phy       ! called by step.F90
+
+   INTEGER ::   nzdf_phy   ! type of vertical closure used 
+   !                       ! associated indicators
+   INTEGER, PARAMETER ::   np_CST = 1   ! Constant Kz
+   INTEGER, PARAMETER ::   np_RIC = 2   ! Richardson number dependent Kz
+   INTEGER, PARAMETER ::   np_TKE = 3   ! Turbulente Kinetic Eenergy closure scheme for Kz
+   INTEGER, PARAMETER ::   np_GLS = 4   ! Generic Length Scale closure scheme for Kz
+      
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , NEMO Consortium (2017)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -51,4 +59,5 @@
       !!
       !! ** Method  :   Read namelist namzdf, control logicals 
+      !!                set horizontal shape and vertical profile of background mixing coef.
       !!----------------------------------------------------------------------
       INTEGER ::   ioptio, ios   ! local integers
@@ -58,28 +67,24 @@
          &             ln_zdfnpc, nn_npc , nn_npcp,                  &     ! convection : npc
          &             ln_zdfddm, rn_avts, rn_hsbfr,                 &     ! double diffusion
-         &             ln_zdftmx,                                    &     ! tidal mixing
-         &             ln_zdfqiao,                                   &     ! surface wave-induced mixing
+         &             ln_zdfswm,                                    &     ! surface  wave-induced mixing
+         &             ln_zdfiwm,                                    &     ! internal  -      -      -
          &             ln_zdfexp, nn_zdfexp,                         &     ! time-stepping
          &             rn_avm0, rn_avt0, nn_avb, nn_havtb                  ! coefficients
-
-
-!!org      NAMELIST/namzdf/ rn_avm0, rn_avt0, nn_avb, nn_havtb, ln_zdfexp, nn_zdfexp,  &
-!!org         &        ln_zdfevd, nn_evdm, rn_avevd, ln_zdfnpc, nn_npc, nn_npcp,       &
-!!org         &        ln_zdfqiao
-      !!----------------------------------------------------------------------
-
+      !!----------------------------------------------------------------------
+      !
+      !                           !==  Namelist  ==!
       REWIND( numnam_ref )              ! Namelist namzdf in reference namelist : Vertical mixing parameters
       READ  ( numnam_ref, namzdf, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf in reference namelist', lwp )
-
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namzdf in reference namelist', lwp )
+      !
       REWIND( numnam_cfg )              ! Namelist namzdf in reference namelist : Vertical mixing parameters
       READ  ( numnam_cfg, namzdf, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf in configuration namelist', lwp )
-      IF(lwm) WRITE ( numond, namzdf )
-
-      IF(lwp) THEN               !* Parameter print
+902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namzdf in configuration namelist', lwp )
+      IF(lwm)   WRITE ( numond, namzdf )
+      !
+      IF(lwp) THEN                      ! Parameter print
          WRITE(numout,*)
-         WRITE(numout,*) 'zdf_phy_init : vertical physics'
-         WRITE(numout,*) '~~~~~~~~'
+         WRITE(numout,*) 'zdf_phy_init: vertical physics'
+         WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namzdf : set vertical mixing mixing parameters'
          WRITE(numout,*) '      vertical closure scheme'
@@ -98,104 +103,33 @@
          WRITE(numout,*) '         maximum avs for dd mixing                 rn_avts = ', rn_avts
          WRITE(numout,*) '         heat/salt buoyancy flux ratio             rn_hsbfr= ', rn_hsbfr
-         WRITE(numout,*) '      surface wave-induced mixing                ln_zdfqiao= ', ln_zdfqiao                                          ! surface wave induced mixing
-         WRITE(numout,*) '      tidal mixing                               ln_zdftmx = ', ln_zdftmx
-         WRITE(numout,*) '      time splitting / backward scheme           ln_zdfexp = ', ln_zdfexp
+         WRITE(numout,*) '      gravity wave-induced mixing'
+         WRITE(numout,*) '         surface  wave (Qiao et al 2010)         ln_zdfswm = ', ln_zdfswm                                          ! surface wave induced mixing
+         WRITE(numout,*) '         internal wave (de Lavergne et al 2017)  ln_zdfiwm = ', ln_zdfiwm
+         WRITE(numout,*) '      time-steping scheme'
+         WRITE(numout,*) '         time splitting (T) / implicit (F)       ln_zdfexp = ', ln_zdfexp
          WRITE(numout,*) '         number of sub-time step (ln_zdfexp=T)   nn_zdfexp = ', nn_zdfexp
          WRITE(numout,*) '      coefficients : '
-         WRITE(numout,*) '      vertical eddy viscosity             rn_avm0   = ', rn_avm0
-         WRITE(numout,*) '      vertical eddy diffusivity           rn_avt0   = ', rn_avt0
-         WRITE(numout,*) '      constant background or profile      nn_avb    = ', nn_avb
-         WRITE(numout,*) '      horizontal variation for avtb       nn_havtb  = ', nn_havtb
-      ENDIF
-
-!!gm      IF(ln_zdfddm) THEN                    ! double diffusive mixing'
-!         avs(:,:,:) = rn_avt0 * wmask(:,:,:) 
-!!gm      ENDIF
-
-      !                          !* Parameter & logical controls
-      !                          !  ----------------------------
-      !
-      !                               ! ... check of vertical mixing scheme on tracers
-      !                                              ==> will be done in trazdf module
-      !
-      !                               ! ... check of mixing coefficient
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) '   vertical mixing option :'
-      ioptio = 0
-      IF( ln_zdfcst ) THEN
-         IF(lwp) WRITE(numout,*) '      constant eddy diffusion coefficients'
-         ioptio = ioptio+1
-      ENDIF
-      IF( ln_zdfric ) THEN
-         IF(lwp) WRITE(numout,*) '      Richardson dependent eddy coefficients'
-         ioptio = ioptio+1
-      ENDIF
-      IF( ln_zdftke ) THEN
-         IF(lwp) WRITE(numout,*) '      TKE dependent eddy coefficients'
-         ioptio = ioptio+1
-      ENDIF
-      IF( ln_zdfgls ) THEN
-         IF(lwp) WRITE(numout,*) '      GLS dependent eddy coefficients'
-         ioptio = ioptio+1
-      ENDIF
-      IF( ioptio == 0 .OR. ioptio > 1 )   &
-         &   CALL ctl_stop( ' one and only one vertical diffusion option has to be defined ' )
-      IF( ( ln_zdfric .OR. ln_zdfgls ) .AND. ln_isfcav )   &
-         &   CALL ctl_stop( ' only zdfcst and zdftke were tested with ice shelves cavities ' )
-      !
-      !                               ! ... Convection
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) '   convection :'
-      !
-#if defined key_top
-      IF( ln_zdfnpc )   CALL ctl_stop( ' zdf_phy_init: npc scheme is not working with key_top' )
-#endif
-      !
-      ioptio = 0
-      IF( ln_zdfnpc ) THEN
-         IF(lwp) WRITE(numout,*) '      use non penetrative convective scheme'
-         ioptio = ioptio+1
-      ENDIF
-      IF( ln_zdfevd ) THEN
-         IF(lwp) WRITE(numout,*) '      use enhanced vertical dif. scheme'
-         ioptio = ioptio+1
-      ENDIF
-      IF( ln_zdftke ) THEN
-         IF(lwp) WRITE(numout,*) '      use the 1.5 turbulent closure'
-      ENDIF
-      IF( ln_zdfgls ) THEN
-         IF(lwp) WRITE(numout,*) '      use the GLS closure scheme'
-      ENDIF
-      IF ( ioptio > 1 )   CALL ctl_stop( ' chose between ln_zdfnpc and ln_zdfevd' )
-      IF( ioptio == 0 .AND. .NOT.( ln_zdftke .OR. ln_zdfgls ) )           &
-         CALL ctl_stop( ' except for TKE or GLS physics, a convection scheme is',   &
-         &              ' required: ln_zdfevd or ln_zdfnpc logicals' )
-
-      !                               !* Background eddy viscosity and diffusivity profil
-      IF( nn_avb == 0 ) THEN                ! Define avmb, avtb from namelist parameter
+         WRITE(numout,*) '         vertical eddy viscosity                 rn_avm0   = ', rn_avm0
+         WRITE(numout,*) '         vertical eddy diffusivity               rn_avt0   = ', rn_avt0
+         WRITE(numout,*) '         constant background or profile          nn_avb    = ', nn_avb
+         WRITE(numout,*) '         horizontal variation for avtb           nn_havtb  = ', nn_havtb
+      ENDIF
+
+      !                          !==  Background eddy viscosity and diffusivity  ==!
+      IF( nn_avb == 0 ) THEN             ! Define avmb, avtb from namelist parameter
          avmb(:) = rn_avm0
          avtb(:) = rn_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(:) = rn_avm0
          avtb(:) = rn_avt0 + ( 3.e-4_wp - 2._wp * rn_avt0 ) * 1.e-4_wp * gdepw_1d(:)   ! m2/s
          IF(ln_sco .AND. lwp)   CALL ctl_warn( 'avtb profile not valid in sco' )
       ENDIF
-      !
-      IF( ln_rstart ) THEN                  !  Read avmb, avtb in restart (if exist)
-         ! if ln_traadv_cen, avmb, avtb have been modified in traadv_cen2 module. 
-         ! To ensure the restartability, avmb & avtb are written in the restart 
-         ! file in traadv_cen2 end read here. 
-         IF( iom_varid( numror, 'avmb', ldstop = .FALSE. ) > 0 ) THEN
-            CALL iom_get( numror, jpdom_unknown, 'avmb', avmb )
-            CALL iom_get( numror, jpdom_unknown, 'avtb', avtb )
-         ENDIF
-      ENDIF
-      !                                     ! 2D shape of the avtb
-      avtb_2d(:,:) = 1.e0                        ! uniform 
-      !
-      IF( nn_havtb == 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
+      !                                  ! 2D shape of the avtb
+      avtb_2d(:,:) = 1._wp                   ! uniform 
+      !
+      IF( nn_havtb == 1 ) THEN               ! decrease avtb by a factor of ten 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
@@ -203,28 +137,53 @@
       ENDIF
       !
-
-!!gm moved into zdf_phy_init
-!
-                            CALL zdf_bfr_init      ! bottom friction
-
-      ioptio = 0                 !==  type of vertical turbulent closure  ==!    (set nzdfphy)
-      !
-!      IF( ln_zdfcst ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_CST   ;   ENDIF
-!      IF( ln_zdfric ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_RIC   ;   CALL zdf_ric_init   ;   ENDIF
-!      IF( ln_zdftke ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_TKE   ;   CALL zdf_tke_init   ;   ENDIF
-!      IF( ln_zdfgls ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_GLS   ;   CALL zdf_gls_init   ;   ENDIF
-
-
-      !
-      IF( ln_zdfric     )   CALL zdf_ric_init      ! Richardson number dependent Kz
-      IF( ln_zdftke     )   CALL zdf_tke_init      ! TKE closure scheme
-      IF( ln_zdfgls     )   CALL zdf_gls_init      ! GLS closure scheme
-      IF( ln_zdftmx     )   CALL zdf_tmx_init      ! tidal vertical mixing
-!!gm
+      !                                   ! set 1st/last level Av to zero one for all
+      avt(:,:,1) = 0._wp   ;   avs(:,:,1) = 0._wp   ;   avm(:,:,1) = 0._wp
+
+      !                          !==  Convection  ==!
+      !
+      IF( ln_zdfnpc .AND. ln_zdfevd )   CALL ctl_stop( 'zdf_phy_init: chose between ln_zdfnpc and ln_zdfevd' )
+      IF( lk_top    .AND. ln_zdfnpc )   CALL ctl_stop( 'zdf_phy_init: npc scheme is not working with key_top' )
+      IF(lwp) THEN
+         WRITE(numout,*)
+         IF    ( ln_zdfnpc ) THEN  ;   WRITE(numout,*) '      convection: use non penetrative convective scheme'
+         ELSEIF( ln_zdfevd ) THEN  ;   WRITE(numout,*) '      convection: use enhanced vertical diffusion scheme'
+         ELSE                      ;   WRITE(numout,*) '      convection: no specific scheme used'
+         ENDIF
+      ENDIF
+
+      IF(lwp) THEN               !==  Double Diffusion Mixing parameterization  ==!   (ddm)
+         WRITE(numout,*)
+         IF( ln_zdfddm ) THEN   ;   WRITE(numout,*) '      use double diffusive mixing: avs /= avt'
+         ELSE                   ;   WRITE(numout,*) '      No  double diffusive mixing: avs = avt'
+         ENDIF
+      ENDIF
+
+      !                          !==  type of vertical turbulent closure  ==!    (set nzdf_phy)
+      ioptio = 0 
+      IF( ln_zdfcst ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_CST   ;   ENDIF
+      IF( ln_zdfric ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_RIC   ;   CALL zdf_ric_init   ;   ENDIF
+      IF( ln_zdftke ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_TKE   ;   CALL zdf_tke_init   ;   ENDIF
+      IF( ln_zdfgls ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_GLS   ;   CALL zdf_gls_init   ;   ENDIF
+      !
+      IF( ioptio /= 1 )    CALL ctl_stop( 'zdf_phy_init: one and only one vertical diffusion option has to be defined ' )
+      IF( ln_isfcav ) THEN
+      IF( ln_zdfric .OR. ln_zdfgls )    CALL ctl_stop( 'zdf_phy_init: zdfric and zdfgls never tested with ice shelves cavities ' )
+      ENDIF
+
+      !                          !== gravity wave-driven mixing  ==!
+      IF( ln_zdfiwm )   CALL zdf_iwm_init       ! internal wave-driven mixing
+      IF( ln_zdfswm )   CALL zdf_swm_init       ! surface  wave-driven mixing
+
+      !                          !== bottom friction  ==!
+      CALL zdf_bfr_init
+      !
+      !                          !== time-stepping  ==!
+      ! Check/update of time stepping done in dynzdf_init/trazdf_init
+      !!gm move it here ?
       !
    END SUBROUTINE zdf_phy_init
 
 
-   SUBROUTINE zdf_phy( kstp )
+   SUBROUTINE zdf_phy( kt )
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE zdf_phy  ***
@@ -238,43 +197,65 @@
       !!                bottom stress.....                               <<<<====verifier !
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kstp   ! ocean time-step index
-      !
-      INTEGER ::   ji, jj, jk    ! dummy loop indice
-      !!----------------------------------------------------------------------
-      !
-                         CALL zdf_bfr( kstp )         ! bottom friction (if quadratic)
-      !                                               ! Vertical eddy viscosity and diffusivity coefficients
-      IF( ln_zdfric  )   CALL zdf_ric ( kstp )             ! Richardson number dependent Kz
-      IF( ln_zdftke  )   CALL zdf_tke ( kstp )             ! TKE closure scheme for Kz
-      IF( ln_zdfgls  )   CALL zdf_gls ( kstp )             ! GLS closure scheme for Kz
-      IF( ln_zdfqiao )   CALL zdf_qiao( kstp )             ! Qiao vertical mixing 
-      !
-      IF( ln_zdfcst  ) THEN                                ! Constant Kz (reset avt, avm[uv] to the background value)
-         avt (:,:,:) = rn_avt0 * wmask (:,:,:)
-         avm (:,:,:) = rn_avm0 * wmask (:,:,:)
-         avmu(:,:,:) = rn_avm0 * wumask(:,:,:)
-         avmv(:,:,:) = rn_avm0 * wvmask(:,:,:)
-      ENDIF
-      !
-      IF( ln_rnf_mouth ) THEN                         ! increase diffusivity at rivers mouths
-         DO jk = 2, nkrnf   ;   avt(:,:,jk) = avt(:,:,jk) + 2._wp * rn_avt_rnf * rnfmsk(:,:) * tmask(:,:,jk)   ;   END DO
-      ENDIF
-      !
-      IF( ln_zdfevd  )   CALL zdf_evd( kstp )         ! enhanced vertical eddy diffusivity
-      !
-      IF( ln_zdfddm  ) THEN                           ! double diffusive mixing
-                         CALL zdf_ddm( kstp )
-      ELSE                                            ! avs=avt
-         DO jk = 2, jpkm1   ;   avs(:,:,jk) = avt(:,:,jk)   ;   END DO
-      ENDIF
-      !
-      IF( ln_zdftmx  )   CALL zdf_tmx( kstp )         ! tidal vertical mixing
-
-                         CALL zdf_mxl( kstp )         ! mixed layer depth
-
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      !
+      INTEGER ::   ji, jj, jk   ! dummy loop indice
+      !! ---------------------------------------------------------------------
+      !
+      CALL zdf_bfr( kt )                        !* bottom friction (if quadratic)
+      !                       
+      SELECT CASE ( nzdf_phy )                  !* Vertical eddy viscosity and diffusivity coefficients at w-points
+      CASE( np_RIC )   ;   CALL zdf_ric( kt )         ! Richardson number dependent Kz
+      CASE( np_TKE )   ;   CALL zdf_tke( kt )         ! TKE closure scheme for Kz
+      CASE( np_GLS )   ;   CALL zdf_gls( kt )         ! GLS closure scheme for Kz
+      CASE( np_CST )                                  ! Constant Kz (reset avt, avm to the background value)
+         avt(2:jpim1,2:jpjm1,1:jpkm1) = rn_avt0 * wmask(2:jpim1,2:jpjm1,1:jpkm1)
+         avm(2:jpim1,2:jpjm1,1:jpkm1) = rn_avm0 * wmask(2:jpim1,2:jpjm1,1:jpkm1)
+      END SELECT
+      !
+      IF( ln_rnf_mouth ) THEN                   !* increase diffusivity at rivers mouths
+         DO jk = 2, nkrnf
+            avt(:,:,jk) = avt(:,:,jk) + 2._wp * rn_avt_rnf * rnfmsk(:,:) * wmask(:,:,jk)
+         END DO
+      ENDIF
+      !
+      IF( ln_zdfevd )   CALL zdf_evd( kt )      !* convection: enhanced vertical eddy diffusivity
+      !
+      !                                         !* double diffusive mixing
+      IF( ln_zdfddm ) THEN                            ! update avt and compute avs
+                        CALL zdf_ddm( kt )
+      ELSE                                            ! same mixing on all tracers
+         avs(2:jpim1,2:jpjm1,1:jpkm1) = avt(2:jpim1,2:jpjm1,1:jpkm1)
+      ENDIF
+      !
+      !                                         !* wave-induced mixing 
+      IF( ln_zdfswm )   CALL zdf_swm( kt )            ! surface  wave (Qiao et al. 2004) 
+      IF( ln_zdfiwm )   CALL zdf_iwm( kt )            ! internal wave (de Lavergne et al 2017)
+
+
+      !                                         !* Lateral boundary conditions (sign unchanged)
+      CALL lbc_lnk( avm, 'W', 1. )
+      CALL lbc_lnk( avt, 'W', 1. )
+      CALL lbc_lnk( avs, 'W', 1. )
+      !
+!!gm  ===>>>   TO BE REMOVED 
+      DO jk = 1, jpkm1            !* vertical eddy viscosity at wu- and wv-points
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * wumask(ji,jj,jk)
+               avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * wvmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions
+!!gm end
+
+
+      CALL zdf_mxl( kt )                        !* mixed layer depth, and level
+
+!!gm  !==>> to be moved in zdftke & zdfgls    
                                                       ! write TKE or GLS information in the restart file
-      IF( lrst_oce .AND. ln_zdftke )   CALL tke_rst( kstp, 'WRITE' )
-      IF( lrst_oce .AND. ln_zdfgls )   CALL gls_rst( kstp, 'WRITE' )
-      !
+      IF( lrst_oce .AND. ln_zdftke )   CALL tke_rst( kt, 'WRITE' )
+      IF( lrst_oce .AND. ln_zdfgls )   CALL gls_rst( kt, 'WRITE' )
+      !      
    END SUBROUTINE zdf_phy
 
Index: anches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfqiao.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfqiao.F90	(revision 7989)
+++ 	(revision )
@@ -1,133 +1,0 @@
-MODULE zdfqiao
-   !!======================================================================
-   !!                       ***  MODULE  zdfqiao  ***
-   !! Qiao module      : vertical mixing enhancement due to surface waves 
-   !!======================================================================
-   !! History :  3.6  !  2014-10  (E. Clementi)  Original code
-   !!----------------------------------------------------------------------
-   !!   zdf_qiao        : compute Qiao parameters
-   !!----------------------------------------------------------------------
-
-   USE in_out_manager  ! I/O manager
-   USE lib_mpp         ! distribued memory computing library
-   USE sbc_oce         ! Surface boundary condition: ocean fields
-   USE zdf_oce
-   USE sbcwave         ! wave module
-   USE dom_oce
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)  
-   
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC zdf_qiao    ! routine called in step
-
-   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: qbv, qbvu, qbvv
-
-   !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
-   !! $Id: $
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
-CONTAINS
-
-   SUBROUTINE zdf_qiao( kt )
-      !!---------------------------------------------------------------------
-      !!                     ***  ROUTINE zdf_qiao ***
-      !!
-      !! ** Purpose :Compute the Qiao term (qbv) to be added to
-      !!             vertical viscosity and diffusivity coeffs.  
-      !!
-      !! ** Method  :qbv = alpha * A * Us(0) * exp (3 * k * z)
-      !!             
-      !! ** action  :Compute the Qiao wave dependent term 
-      !!             only if ln_wave=.true.
-      !!               
-      !!---------------------------------------------------------------------
-      INTEGER, INTENT( in  ) ::  kt   ! ocean time step
-      !
-      INTEGER :: jj, ji, jk   ! dummy loop indices
-      !!---------------------------------------------------------------------
-      !
-      IF( kt == nit000 ) THEN                   ! First call kt=nit000 !
-         IF( .NOT. ( ln_wave .AND. ln_sdw ) )   &
-            &   CALL ctl_stop ( 'Ask for wave Qiao enhanced turbulence but ln_wave   &
-            &                    and ln_sdw have to be activated')
-         IF( zdf_qiao_alloc() /= 0 )   &
-            &   CALL ctl_stop( 'STOP', 'zdf_qiao : unable to allocate arrays' )
-      ENDIF
-
-      !
-      ! Compute the Qiao term Bv (qbv) to be added to
-      ! vertical viscosity and diffusivity
-      ! qbv = alpha * A * Us(0) * exp (3 * k * z)
-      ! alpha here is set to 1
-      !---------------------------------------------------------------------------------
-      !
-!!gm Comment: I don't understand the use of min of 4 gdepw_n to define a quantity at w-point
-!!gm                       ==>> this is an error....
-      DO jk = 1, jpk
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1
-               qbv(ji,jj,jk) = 1.0 * 0.353553 * hsw(ji,jj) * tsd2d(ji,jj) *             &
-            &                  EXP(3.0 * wnum(ji,jj) *                                  &                     
-            &                  (-MIN( gdepw_n(ji  ,jj  ,jk), gdepw_n(ji+1,jj  ,jk),     &
-            &                         gdepw_n(ji  ,jj+1,jk), gdepw_n(ji+1,jj+1,jk))))   &
-            &                          * wmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-      !
-      CALL lbc_lnk( qbv, 'W', 1. )   ! Lateral boundary conditions
-         
-      !
-      ! Interpolate Qiao parameter qbv into the grid_U and grid_V
-      !----------------------------------------------------------
-      !
-      DO jk = 1, jpk
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1
-               qbvu(ji,jj,jk) = 0.5 * wumask(ji,jj,jk)  *              &  
-            &                  ( qbv(ji,jj,jk) + qbv(ji+1,jj  ,jk) )
-               qbvv(ji,jj,jk) = 0.5 * wvmask(ji,jj,jk)  *              &
-            &                  ( qbv(ji,jj,jk) + qbv(ji  ,jj+1,jk) )
-            END DO
-         END DO
-      END DO
-      ! 
-      CALL lbc_lnk( qbvu, 'U', 1. ) ; CALL lbc_lnk( qbvv, 'V', 1. )   ! Lateral boundary conditions
-
-      ! Enhance vertical mixing coeff.         
-      !-------------------------------
-      !
-!!gm  with double diffusion activated, avs is not updated...
-!!gm                    =====>>> BUG
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               avmu(ji,jj,jk) = ( avmu(ji,jj,jk) + qbvu(ji,jj,jk) ) * umask(ji,jj,jk)
-               avmv(ji,jj,jk) = ( avmv(ji,jj,jk) + qbvv(ji,jj,jk) ) * vmask(ji,jj,jk)
-               avt (ji,jj,jk) = ( avt (ji,jj,jk) + qbv (ji,jj,jk) ) * tmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-      !
-   END SUBROUTINE zdf_qiao
-
-
-   INTEGER FUNCTION zdf_qiao_alloc()
-      !!----------------------------------------------------------------------
-      !!                ***  FUNCTION zdf_qiao_alloc  ***
-      !!----------------------------------------------------------------------
-      ALLOCATE( qbv(jpi,jpj,jpk), qbvu(jpi,jpj,jpk), qbvv(jpi,jpj,jpk),   &
-         &      STAT = zdf_qiao_alloc )
-      !
-      IF( lk_mpp             )  CALL mpp_sum ( zdf_qiao_alloc )
-      IF( zdf_qiao_alloc > 0 )  CALL ctl_warn('zdf_qiao_alloc: allocation of arrays failed.')
-      !
-   END FUNCTION zdf_qiao_alloc
-      
-   !!======================================================================
-END MODULE zdfqiao
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90	(revision 7989)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90	(revision 7990)
@@ -5,11 +5,12 @@
    !!                 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
+   !! 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 
    !!----------------------------------------------------------------------
 
@@ -21,12 +22,14 @@
    USE dom_oce        ! ocean space and time domain variables
    USE zdf_oce        ! ocean vertical physics
+   USE phycst         ! physical constants
+   USE sbc_oce,  ONLY :   taum
+   USE eosbn2 ,  ONLY :   neos
+   !
    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 : neos
 
    IMPLICIT NONE
@@ -36,16 +39,14 @@
    PUBLIC   zdf_ric_init    ! called by opa.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
-
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tmric   !: coef. for the horizontal mean at t-point
+   !                        !!* 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
@@ -57,15 +58,4 @@
    !!----------------------------------------------------------------------
 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 )
@@ -105,126 +95,96 @@
       !!      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.
       !!----------------------------------------------------------------------
-      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))
-         ! -----------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step
+      !!
+      INTEGER  ::   ji, jj, jk                 ! dummy loop indices
+      LOGICAL  ::   ll_av_weight = .TRUE.      ! local logical
+      REAL(wp) ::   zsh2, zcfRi, zav, zustar   ! local scalars
+      REAL(wp), DIMENSION(jpi,jpj) ::   zdu2, zdv2, zh_ekm   ! 2D workspace
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )   CALL timing_start('zdf_ric')
+      !
+      DO jk = 2, jpkm1        !==  avm and avt = F(Richardson number)  ==!
+         !
+      IF( ll_av_weight ) THEN    !== viscosity weighted shear & stratification terms
+         !
+         DO jj = 1, jpjm1           !* viscosity weighted shear term
+            DO ji = 1, jpim1
+               zdu2(ji,jj) = 0.5 * ( avm(ji,jj,jk  ) + avm(ji+1,jj,jk) ) * wumask(ji,jj,jk)      &
+                  &              * (  un(ji,jj,jk-1) -  un(ji  ,jj,jk) )                         &
+                  &              * (  ub(ji,jj,jk-1) -  ub(ji  ,jj,jk) ) / (  e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) )
+               zdv2(ji,jj) = 0.5 * ( avm(ji,jj,jk  ) + avm(ji,jj+1,jk) ) * wumask(ji,jj,jk)      &
+                  &              * (  vn(ji,jj,jk-1) -  vn(ji,jj  ,jk) )                         &
+                  &              * (  vb(ji,jj,jk-1) -  vb(ji,jj  ,jk) ) / (  e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) )
+            END DO
+         END DO
          DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1
-               zcoef = 0.5 / e3w_n(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, fs_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 = fs_2, fs_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 = fs_2, fs_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.
-      ! -------------------------------------------------------
-      zflageos = ( 0.5 + SIGN( 0.5, neos - 1. ) ) * rau0
-      DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1
-            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_t(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
+            DO ji = 2, jpim1        !* Richardson number dependent coefficient
+               !                          ! shear of horizontal velocity
+               zsh2  =  ( zdu2(ji-1,jj) + zdu2(ji,jj) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
+                  &   + ( zdv2(ji,jj-1) + zdv2(ji,jj) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )    
+               !                          ! coefficient = F(richardson number)
+               zcfRi = 1._wp / (  1._wp + rn_alp * MAX(  0._wp , avt(ji,jj,jk) * rn2(ji,jj,jk) / ( zsh2 + 1.e-20 ) )  )
+               !
+               !                    !* Vertical eddy viscosity and diffusivity coefficients
+               zav = rn_avmri * zcfRi**nn_ric
+               avm(ji,jj,jk) = MAX(  zav         , avmb(jk)  ) * wmask(ji,jj,jk)
+               avt(ji,jj,jk) = MAX(  zav * zcfRi , avtb(jk)  ) * wmask(ji,jj,jk)
+            END DO
+         END DO
+      ELSE                    !==  classical Richardson number definition  ==!
+         DO jj = 1, jpjm1           !* shear term
+            DO ji = 1, jpim1
+               zdu2(ji,jj) = 0.5 * (  un(ji,jj,jk-1) -  un(ji  ,jj,jk) )                         &
+                  &              * (  ub(ji,jj,jk-1) -  ub(ji  ,jj,jk) ) / (  e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) )
+               zdv2(ji,jj) = 0.5 * (  vn(ji,jj,jk-1) -  vn(ji,jj  ,jk) )                         &
+                  &              * (  vb(ji,jj,jk-1) -  vb(ji,jj  ,jk) ) / (  e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) )
+            END DO
+         END DO
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1        !* Richardson number dependent coefficient
+               !                          ! shear of horizontal velocity
+               zsh2  =  ( zdu2(ji-1,jj) + zdu2(ji,jj) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
+                  &   + ( zdv2(ji,jj-1) + zdv2(ji,jj) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )    
+               !                          ! coefficient = F(richardson number)
+               zcfRi = 1._wp / (  1._wp + rn_alp * MAX(  0._wp , rn2(ji,jj,jk) / ( zsh2 + 1.e-20 ) )  )
+               !
+               !                    !* Vertical eddy viscosity and diffusivity coefficients
+               zav = rn_avmri * zcfRi**nn_ric
+               avm(ji,jj,jk) = MAX(  zav         , avmb(jk)  ) * wmask(ji,jj,jk)
+               avt(ji,jj,jk) = MAX(  zav * zcfRi , avtb(jk)  ) * wmask(ji,jj,jk)
+            END DO
+         END DO
+      ENDIF
+         !
       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 = 2, jpjm1
-         DO ji = fs_2, fs_jpim1
-            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 = 2, jpjm1
-            DO ji = fs_2, fs_jpim1
-               IF( gdept_n(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 = 2, jpjm1
-            DO ji = fs_2, fs_jpim1
-               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')
+      !
+      IF( ln_mldw ) THEN      !==  set a minimum value in the Ekman layer  ==!
+         !
+         DO jj = 2, jpjm1        !* Ekman depth
+            DO ji = 2, jpim1
+               zustar = SQRT( taum(ji,jj) * r1_rau0 )
+               zh_ekm(ji,jj) = rn_ekmfc * zustar / ( ABS( ff_t(ji,jj) ) + rsmall )     ! Ekman depth
+               zh_ekm(ji,jj) = MAX(  rn_mldmin , MIN( zh_ekm(ji,jj) , rn_mldmax )  )   ! set allowed rang
+            END DO
+         END DO
+         DO jk = 2, jpkm1        !* minimum mixing coeff. within the Ekman layer
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1
+                  IF( gdept_n(ji,jj,jk) < zh_ekm(ji,jj) ) THEN
+                     avm(ji,jj,jk) = MAX(  avm(ji,jj,jk), rn_wvmix  ) * wmask(ji,jj,jk)
+                     avt(ji,jj,jk) = MAX(  avt(ji,jj,jk), rn_wtmix  ) * wmask(ji,jj,jk)
+                  ENDIF
+               END DO
+            END DO
+         END DO
+      END IF
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('zdf_ric')
       !
    END SUBROUTINE zdf_ric
@@ -264,34 +224,19 @@
          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
+         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
       !
-      !                              ! 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)
+         avt(:,:,jk) = avtb(jk) * tmask(:,:,jk)
+         avm(:,:,jk) = avmb(jk) * umask(:,:,jk)
       END DO
       !
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfswm.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfswm.F90	(revision 7990)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfswm.F90	(revision 7990)
@@ -0,0 +1,95 @@
+MODULE zdfswm
+   !!======================================================================
+   !!                       ***  MODULE  zdfswm  ***
+   !! vertical physics :   surface wave-induced mixing 
+   !!======================================================================
+   !! History :  3.6  !  2014-10  (E. Clementi)  Original code
+   !!            4.0  !  2017-04  (G. Madec)  debug + simplifications
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   zdf_swm      : update Kz due to surface wave-induced mixing 
+   !!   zdf_swm_init : initilisation
+   !!----------------------------------------------------------------------
+   USE dom_oce        ! ocean domain variable
+   USE zdf_oce        ! vertical physics: mixing coefficients
+   USE sbc_oce        ! Surface boundary condition: ocean fields
+   USE sbcwave        ! wave module
+   !
+   USE in_out_manager ! I/O manager
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)  
+   USE lib_mpp        ! distribued memory computing library
+   
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC zdf_swm         ! routine called in zdp_phy
+   PUBLIC zdf_swm_init    ! routine called in zdf_phy_init
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , NEMO Consortium (2017) 
+   !! $Id:$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE zdf_swm( kt )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE zdf_swm ***
+      !!
+      !! ** Purpose :Compute the swm term (qbv) to be added to
+      !!             vertical viscosity and diffusivity coeffs.  
+      !!
+      !! ** Method  :   Compute the swm term Bv (zqb) and added it to
+      !!               vertical viscosity and diffusivity coefficients
+      !!                   zqb = alpha * A * Us(0) * exp (3 * k * z)
+      !!               where alpha is set here to 1
+      !!             
+      !! ** action  :    avt, avs, avm updated by the surface wave-induced mixing
+      !!                               (inner domain only)
+      !!               
+      !! reference : Qiao et al. GRL, 2004
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time step
+      !
+      INTEGER ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp)::   zcoef, zqb   ! local scalar
+      !!---------------------------------------------------------------------
+      !
+      zcoef = 1._wp * 0.353553_wp
+      DO jk = 2, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               zqb = zcoef * hsw(ji,jj) * tsd2d(ji,jj) * EXP( -3. * wnum(ji,jj) * gdepw_n(ji,jj,jk) ) * wmask(ji,jj,jk)
+               !
+               avt(ji,jj,jk) = avt(ji,jj,jk) + zqb
+               avs(ji,jj,jk) = avs(ji,jj,jk) + zqb
+               avm(ji,jj,jk) = avm(ji,jj,jk) + zqb
+            END DO
+         END DO
+      END DO
+      !
+   END SUBROUTINE zdf_swm
+   
+   
+   SUBROUTINE zdf_swm_init
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE zdf_swm_init ***
+      !!
+      !! ** Purpose :   surface wave-induced mixing initialisation  
+      !!
+      !! ** Method  :   check the availability of surface wave fields
+      !!---------------------------------------------------------------------
+      !
+      IF(lwp) THEN                  ! Control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'zdf_swm_init : surface wave-driven mixing'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+      ENDIF
+      IF(  .NOT.ln_wave .OR.   &
+         & .NOT.ln_sdw    )   CALL ctl_stop ( 'zdf_swm_init: ln_zdfswm=T but ln_wave and ln_sdw /= T')
+      !
+   END SUBROUTINE zdf_swm_init
+
+   !!======================================================================
+END MODULE zdfswm
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 7989)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 7990)
@@ -27,5 +27,5 @@
    !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
    !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability
-   !!            4.0  !  2017-04  (G. Madec)  Remove CPP keys
+   !!            4.0  !  2017-04  (G. Madec)  remove CPP ddm key & avm at t-point only 
    !!----------------------------------------------------------------------
 
@@ -159,5 +159,5 @@
       !!
       !! ** Action  :   compute en (now turbulent kinetic energy)
-      !!                update avt, avmu, avmv (before vertical eddy coef.)
+      !!                update avt, avm (before vertical eddy coef.)
       !!
       !! References : Gaspar et al., JGR, 1990,
@@ -175,19 +175,13 @@
 #endif
       !
-      IF( kt /= nit000 ) THEN   ! restore before value to compute tke
-         avt (:,:,:) = avt_k (:,:,:) 
-         avm (:,:,:) = avm_k (:,:,:) 
-         avmu(:,:,:) = avmu_k(:,:,:) 
-         avmv(:,:,:) = avmv_k(:,:,:) 
-      ENDIF 
-      !
-      CALL tke_tke      ! now tke (en)
-      !
-      CALL tke_avn      ! now avt, avm, avmu, avmv
-      !
-      avt_k (:,:,:) = avt (:,:,:) 
-      avm_k (:,:,:) = avm (:,:,:) 
-      avmu_k(:,:,:) = avmu(:,:,:) 
-      avmv_k(:,:,:) = avmv(:,:,:) 
+      avt(:,:,:) = avt_k(:,:,:)     ! restore before Kz computed with TKE alone
+      avm(:,:,:) = avm_k(:,:,:) 
+      !
+      CALL tke_tke                  ! now tke (en)
+      !
+      CALL tke_avn                  ! now avt, avm, dissl
+      !
+      avt_k(:,:,:) = avt(:,:,:)     ! store Kz computed with TKE alone
+      avm_k(:,:,:) = avm(:,:,:) 
       !
 #if defined key_agrif
@@ -195,5 +189,7 @@
       IF( .NOT.Agrif_Root() )   CALL Agrif_Update_Tke( kt )      ! children only
 #endif      
-     ! 
+      ! 
+      IF( lrst_oce )   CALL tke_rst( kt, 'WRITE' )
+      !
   END SUBROUTINE zdf_tke
 
@@ -207,5 +203,5 @@
       !! ** Method  : - TKE surface boundary condition
       !!              - source term due to Langmuir cells (Axell JGR 2002) (ln_lc=T)
-      !!              - source term due to shear (saved in avmu, avmv arrays)
+      !!              - source term due to shear (= Kz dz[Ub] * dz[Un] )
       !!              - Now TKE : resolution of the TKE equation by inverting
       !!                a tridiagonal linear system by a "methode de chasse"
@@ -213,6 +209,4 @@
       !!
       !! ** Action  : - en : now turbulent kinetic energy)
-      !!              - avmu, avmv : production of TKE by shear at u and v-points
-      !!                (= Kz dz[Ub] * dz[Un] )
       !! ---------------------------------------------------------------------
       INTEGER  ::   ji, jj, jk                      ! dummy loop arguments
@@ -228,5 +222,5 @@
       REAL(wp) ::   zzd_up, zzd_lw                  !    -         -
 !!bfr      REAL(wp) ::   zebot                           !    -         -
-      INTEGER , POINTER, DIMENSION(:,:  ) ::   imlc
+      INTEGER , DIMENSION(jpi,jpj) ::   imlc
       REAL(wp), POINTER, DIMENSION(:,:  ) ::   zhlc
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zpelc, zdiag, zd_up, zd_lw, z3du, z3dv
@@ -236,5 +230,4 @@
       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, z3du, z3dv ) 
@@ -473,7 +466,7 @@
          END DO
       ENDIF
+!!gm not sure we need this lbc ....
       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, z3du, z3dv ) 
@@ -516,11 +509,10 @@
       !!      with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise.
       !!
-      !! ** Action  : - avt : now vertical eddy diffusivity (w-point)
-      !!              - avmu, avmv : now vertical eddy viscosity at uw- and vw-points
+      !! ** Action  : - avt, avm : now vertical eddy diffusivity and viscosity (w-point)
       !!----------------------------------------------------------------------
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zrn2, zraug, zcoef, zav     ! local scalars
-      REAL(wp) ::   zdku, zri, zsqen            !   -      -
-      REAL(wp) ::   zdkv, zemxl, zemlm, zemlp   !   -      -
+      REAL(wp) ::   zrn2, zraug, zcoef, zav   ! local scalars
+      REAL(wp) ::   zdku,   zdkv, zsqen       !   -      -
+      REAL(wp) ::   zemxl, zemlm, zemlp       !   -      -
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zmpdl, zmxlm, zmxld
       !!--------------------------------------------------------------------
@@ -645,5 +637,5 @@
 
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      !                     !  Vertical eddy viscosity and diffusivity  (avmu, avmv, avt)
+      !                     !  Vertical eddy viscosity and diffusivity  (avm and avt)
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       DO jk = 1, jpkm1            !* vertical eddy viscosity & diffivity at w-points
@@ -658,15 +650,5 @@
          END DO
       END DO
-      CALL lbc_lnk( avm, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
-      !
-      DO jk = 2, jpkm1            !* vertical eddy viscosity at wu- and wv-points
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * wumask(ji,jj,jk)
-               avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * wvmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions
+      !
       !
       IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt
@@ -677,7 +659,4 @@
 # if defined key_c1d
                   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
@@ -685,10 +664,11 @@
          END DO
       ENDIF
-      CALL lbc_lnk( avt, 'W', 1. )                      ! Lateral boundary conditions on avt  (sign unchanged)
+!!gm I'm sure this is needed to compute the shear term at next time-step
+      CALL lbc_lnk( avm, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
+      CALL lbc_lnk( avt, 'W', 1. )      ! Lateral boundary conditions on avt  (sign unchanged)
 
       IF(ln_ctl) THEN
-         CALL prt_ctl( tab3d_1=en  , clinfo1=' tke  - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk)
-         CALL prt_ctl( tab3d_1=avmu, clinfo1=' tke  - u: ', mask1=umask,                   &
-            &          tab3d_2=avmv, clinfo2=       ' v: ', mask2=vmask, ovlap=1, kdim=jpk )
+         CALL prt_ctl( tab3d_1=en , clinfo1=' tke  - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk)
+         CALL prt_ctl( tab3d_1=avm, clinfo1=' tke  - m: ', ovlap=1, kdim=jpk )
       ENDIF
       !
@@ -759,10 +739,9 @@
       ENDIF
       !
-      IF( ln_zdftmx ) THEN          ! Internal wave driven mixing
-         !                          ! specific values of rn_emin & rmxl_min are used
-         rn_emin  = 1.e-10_wp
-         rmxl_min = 1.e-03_wp
+      IF( ln_zdfiwm ) THEN          ! Internal wave-driven mixing
+         rn_emin  = 1.e-10_wp             ! specific values of rn_emin & rmxl_min are used
+         rmxl_min = 1.e-03_wp             ! associated avt minimum = molecular salt diffusivity (10^-9 m2/s)
          IF(lwp) WRITE(numout,*) '      Internal wave-driven mixing case:   force   rn_emin = 1.e-10 and rmxl_min = 1.e-3 '
-      ELSE
+      ELSE                          ! standard case : associated avt minimum = molecular viscosity (10^-6 m2/s)
          rmxl_min = 1.e-6_wp / ( rn_ediff * SQRT( rn_emin ) )    ! resulting minimum length to recover molecular viscosity
          IF(lwp) WRITE(numout,*) '      minimum mixing length with your parameters rmxl_min = ', rmxl_min
@@ -798,6 +777,4 @@
          avt (:,:,jk) = avtb(jk) * wmask (:,:,jk)
          avm (:,:,jk) = avmb(jk) * wmask (:,:,jk)
-         avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk)
-         avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk)
       END DO
       dissl(:,:,:) = 1.e-12_wp
@@ -821,6 +798,6 @@
      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
      !
-     INTEGER ::   jit, jk   ! dummy loop indices
-     INTEGER ::   id1, id2, id3, id4, id5, id6   ! local integers
+     INTEGER ::   jit, jk              ! dummy loop indices
+     INTEGER ::   id1, id2, id3, id4   ! local integers
      !!----------------------------------------------------------------------
      !
@@ -829,40 +806,33 @@
         IF( ln_rstart ) THEN                   !* Read the restart file
            id1 = iom_varid( numror, 'en'   , ldstop = .FALSE. )
-           id2 = iom_varid( numror, 'avt'  , ldstop = .FALSE. )
-           id3 = iom_varid( numror, 'avm'  , ldstop = .FALSE. )
-           id4 = iom_varid( numror, 'avmu' , ldstop = .FALSE. )
-           id5 = iom_varid( numror, 'avmv' , ldstop = .FALSE. )
-           id6 = iom_varid( numror, 'dissl', ldstop = .FALSE. )
+           id2 = iom_varid( numror, 'avt_k', ldstop = .FALSE. )
+           id3 = iom_varid( numror, 'avm_k', ldstop = .FALSE. )
+           id4 = iom_varid( numror, 'dissl', ldstop = .FALSE. )
            !
            IF( id1 > 0 ) THEN                       ! 'en' exists
               CALL iom_get( numror, jpdom_autoglo, 'en', en )
-              IF( MIN( id2, id3, id4, id5, id6 ) > 0 ) THEN        ! all required arrays exist
-                 CALL iom_get( numror, jpdom_autoglo, 'avt'  , avt   )
-                 CALL iom_get( numror, jpdom_autoglo, 'avm'  , avm   )
-                 CALL iom_get( numror, jpdom_autoglo, 'avmu' , avmu  )
-                 CALL iom_get( numror, jpdom_autoglo, 'avmv' , avmv  )
+              IF( MIN( id2, id3, id4 ) > 0 ) THEN        ! all required arrays exist
+                 CALL iom_get( numror, jpdom_autoglo, 'avt_k', avt_k )
+                 CALL iom_get( numror, jpdom_autoglo, 'avm_k', avm_k )
                  CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl )
               ELSE                                                 ! one at least array is missing
-                 CALL tke_avn                                          ! compute avt, avm, avmu, avmv and dissl (approximation)
+                 CALL tke_avn                                          ! compute avt, avm, and dissl (approximation)
               ENDIF
            ELSE                                     ! No TKE array found: initialisation
               IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without tke scheme, en computed by iterative loop'
-              en (:,:,:) = rn_emin * tmask(:,:,:)
-              CALL tke_avn                               ! recompute avt, avm, avmu, avmv and dissl (approximation)
-              !
-              avt_k (:,:,:) = avt (:,:,:)
-              avm_k (:,:,:) = avm (:,:,:)
-              avmu_k(:,:,:) = avmu(:,:,:)
-              avmv_k(:,:,:) = avmv(:,:,:)
+              en (:,:,:) = rn_emin * wmask(:,:,:)
+              CALL tke_avn                               ! recompute avt, avm, and dissl (approximation)
+              avt_k(:,:,:) = avt(:,:,:)
+              avm_k(:,:,:) = avm(:,:,:)
               !
               DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_tke( jit )   ;   END DO
+              avt_k(:,:,:) = avt(:,:,:)
+              avm_k(:,:,:) = avm(:,:,:)
            ENDIF
         ELSE                                   !* Start from rest
            en(:,:,:) = rn_emin * tmask(:,:,:)
            DO jk = 1, jpk                           ! set the Kz to the background value
-              avt (:,:,jk) = avtb(jk) * wmask (:,:,jk)
-              avm (:,:,jk) = avmb(jk) * wmask (:,:,jk)
-              avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk)
-              avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk)
+              avt_k(:,:,jk) = avtb(jk) * wmask (:,:,jk)
+              avm_k(:,:,jk) = avmb(jk) * wmask (:,:,jk)
            END DO
         ENDIF
@@ -871,10 +841,8 @@
         !                                   ! -------------------
         IF(lwp) WRITE(numout,*) '---- tke-rst ----'
-        CALL iom_rstput( kt, nitrst, numrow, 'en'   , en     )
-        CALL iom_rstput( kt, nitrst, numrow, 'avt'  , avt_k  )
-        CALL iom_rstput( kt, nitrst, numrow, 'avm'  , avm_k  )
-        CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k )
-        CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv_k )
-        CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl  )
+        CALL iom_rstput( kt, nitrst, numrow, 'en'   , en    )
+        CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k )
+        CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k )
+        CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl )
         !
      ENDIF
Index: anches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90	(revision 7989)
+++ 	(revision )
@@ -1,495 +1,0 @@
-MODULE zdftmx
-   !!========================================================================
-   !!                       ***  MODULE  zdftmx  ***
-   !! Ocean physics: Internal gravity wave-driven vertical mixing
-   !!========================================================================
-   !! History :  1.0  !  2004-04  (L. Bessieres, G. Madec)  Original code
-   !!             -   !  2006-08  (A. Koch-Larrouy)  Indonesian strait
-   !!            3.3  !  2010-10  (C. Ethe, G. Madec)  reorganisation of initialisation phase
-   !!            3.6  !  2016-03  (C. de Lavergne)  New param: internal wave-driven mixing 
-   !!            4.0  !  2017-04  (G. Madec)  Remove the old tidal mixing param. and key zdftmx(_new)
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   zdf_tmx       : global     momentum & tracer Kz with wave induced Kz
-   !!   zdf_tmx_init  : global     momentum & tracer Kz with wave induced Kz
-   !!----------------------------------------------------------------------
-   USE oce            ! ocean dynamics and tracers variables
-   USE dom_oce        ! ocean space and time domain variables
-   USE zdf_oce        ! ocean vertical physics variables
-   USE zdfddm         ! ocean vertical physics: double diffusive mixing
-   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
-   USE eosbn2         ! ocean equation of state
-   USE phycst         ! physical constants
-   USE prtctl         ! Print control
-   USE in_out_manager ! I/O manager
-   USE iom            ! I/O Manager
-   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)  
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   zdf_tmx         ! called in step module 
-   PUBLIC   zdf_tmx_init    ! called in nemogcm module 
-   PUBLIC   zdf_tmx_alloc   ! called in nemogcm module
-
-   !                       !!* Namelist  namzdf_tmx : internal wave-driven mixing *
-   INTEGER  ::  nn_zpyc     ! pycnocline-intensified mixing energy proportional to N (=1) or N^2 (=2)
-   LOGICAL  ::  ln_mevar    ! variable (=T) or constant (=F) mixing efficiency
-   LOGICAL  ::  ln_tsdiff   ! account for differential T/S wave-driven mixing (=T) or not (=F)
-
-   REAL(wp) ::  r1_6 = 1._wp / 6._wp
-
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ebot_tmx     ! power available from high-mode wave breaking (W/m2)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   epyc_tmx     ! power available from low-mode, pycnocline-intensified wave breaking (W/m2)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ecri_tmx     ! power available from low-mode, critical slope wave breaking (W/m2)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hbot_tmx     ! WKB decay scale for high-mode energy dissipation (m)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hcri_tmx     ! decay scale for low-mode critical slope dissipation (m)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   emix_tmx     ! local energy density available for mixing (W/kg)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   bflx_tmx     ! buoyancy flux Kz * N^2 (W/kg)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   pcmap_tmx    ! vertically integrated buoyancy flux (W/m2)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zav_ratio    ! S/T diffusivity ratio (only for ln_tsdiff=T)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zav_wave     ! Internal wave-induced diffusivity
-
-   !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 4.0 , NEMO Consortium (2016)
-   !! $Id$
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   INTEGER FUNCTION zdf_tmx_alloc()
-      !!----------------------------------------------------------------------
-      !!                ***  FUNCTION zdf_tmx_alloc  ***
-      !!----------------------------------------------------------------------
-      ALLOCATE(     ebot_tmx(jpi,jpj),  epyc_tmx(jpi,jpj),  ecri_tmx(jpi,jpj)    ,   &
-      &             hbot_tmx(jpi,jpj),  hcri_tmx(jpi,jpj),  emix_tmx(jpi,jpj,jpk),   &
-      &         bflx_tmx(jpi,jpj,jpk), pcmap_tmx(jpi,jpj), zav_ratio(jpi,jpj,jpk),   & 
-      &         zav_wave(jpi,jpj,jpk), STAT=zdf_tmx_alloc     )
-      !
-      IF( lk_mpp             )   CALL mpp_sum ( zdf_tmx_alloc )
-      IF( zdf_tmx_alloc /= 0 )   CALL ctl_warn('zdf_tmx_alloc: failed to allocate arrays')
-   END FUNCTION zdf_tmx_alloc
-
-
-   SUBROUTINE zdf_tmx( kt )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE zdf_tmx  ***
-      !!                   
-      !! ** Purpose :   add to the vertical mixing coefficients the effect of
-      !!              breaking internal waves.
-      !!
-      !! ** Method  : - internal wave-driven vertical mixing is given by:
-      !!                  Kz_wave = min(  100 cm2/s, f(  Reb = emix_tmx /( Nu * N^2 )  )
-      !!              where emix_tmx is the 3D space distribution of the wave-breaking 
-      !!              energy and Nu the molecular kinematic viscosity.
-      !!              The function f(Reb) is linear (constant mixing efficiency)
-      !!              if the namelist parameter ln_mevar = F and nonlinear if ln_mevar = T.
-      !!
-      !!              - Compute emix_tmx, the 3D power density that allows to compute
-      !!              Reb and therefrom the wave-induced vertical diffusivity.
-      !!              This is divided into three components:
-      !!                 1. Bottom-intensified low-mode dissipation at critical slopes
-      !!                     emix_tmx(z) = ( ecri_tmx / rau0 ) * EXP( -(H-z)/hcri_tmx )
-      !!                                   / ( 1. - EXP( - H/hcri_tmx ) ) * hcri_tmx
-      !!              where hcri_tmx is the characteristic length scale of the bottom 
-      !!              intensification, ecri_tmx a map of available power, and H the ocean depth.
-      !!                 2. Pycnocline-intensified low-mode dissipation
-      !!                     emix_tmx(z) = ( epyc_tmx / rau0 ) * ( sqrt(rn2(z))^nn_zpyc )
-      !!                                   / SUM( sqrt(rn2(z))^nn_zpyc * e3w(z) )
-      !!              where epyc_tmx is a map of available power, and nn_zpyc
-      !!              is the chosen stratification-dependence of the internal wave
-      !!              energy dissipation.
-      !!                 3. WKB-height dependent high mode dissipation
-      !!                     emix_tmx(z) = ( ebot_tmx / rau0 ) * rn2(z) * EXP(-z_wkb(z)/hbot_tmx)
-      !!                                   / SUM( rn2(z) * EXP(-z_wkb(z)/hbot_tmx) * e3w(z) )
-      !!              where hbot_tmx is the characteristic length scale of the WKB bottom 
-      !!              intensification, ebot_tmx is a map of available power, and z_wkb is the
-      !!              WKB-stretched height above bottom defined as
-      !!                    z_wkb(z) = H * SUM( sqrt(rn2(z'>=z)) * e3w(z'>=z) )
-      !!                                 / SUM( sqrt(rn2(z'))    * e3w(z')    )
-      !!
-      !!              - update the model vertical eddy viscosity and diffusivity: 
-      !!                     avt  = avt  +    av_wave
-      !!                     avm  = avm  +    av_wave
-      !!                     avmu = avmu + mi(av_wave)
-      !!                     avmv = avmv + mj(av_wave)
-      !!
-      !!              - if namelist parameter ln_tsdiff = T, account for differential mixing:
-      !!                     avs  = avt  +    av_wave * diffusivity_ratio(Reb)
-      !!
-      !! ** Action  : - Define emix_tmx used to compute internal wave-induced mixing
-      !!              - avt, avs, avm, avmu, avmv increased by internal wave-driven mixing    
-      !!
-      !! References :  de Lavergne et al. 2015, JPO; 2016, in prep.
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step 
-      !
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   ztpc         ! scalar workspace
-      REAL(wp), DIMENSION(:,:)  , POINTER ::  zfact     ! Used for vertical structure
-      REAL(wp), DIMENSION(:,:)  , POINTER ::  zhdep     ! Ocean depth
-      REAL(wp), DIMENSION(:,:,:), POINTER ::  zwkb      ! WKB-stretched height above bottom
-      REAL(wp), DIMENSION(:,:,:), POINTER ::  zweight   ! Weight for high mode vertical distribution
-      REAL(wp), DIMENSION(:,:,:), POINTER ::  znu_t     ! Molecular kinematic viscosity (T grid)
-      REAL(wp), DIMENSION(:,:,:), POINTER ::  znu_w     ! Molecular kinematic viscosity (W grid)
-      REAL(wp), DIMENSION(:,:,:), POINTER ::  zReb      ! Turbulence intensity parameter
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )   CALL timing_start('zdf_tmx')
-      !
-      CALL wrk_alloc( jpi,jpj,       zfact, zhdep )
-      CALL wrk_alloc( jpi,jpj,jpk,   zwkb, zweight, znu_t, znu_w, zReb )
-
-      !                          ! ----------------------------- !
-      !                          !  Internal wave-driven mixing  !  (compute zav_wave)
-      !                          ! ----------------------------- !
-      !                             
-      !                        !* Critical slope mixing: distribute energy over the time-varying ocean depth,
-      !                                                 using an exponential decay from the seafloor.
-      DO jj = 1, jpj                ! part independent of the level
-         DO ji = 1, jpi
-            zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1)       ! depth of the ocean
-            zfact(ji,jj) = rau0 * (  1._wp - EXP( -zhdep(ji,jj) / hcri_tmx(ji,jj) )  )
-            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ecri_tmx(ji,jj) / zfact(ji,jj)
-         END DO
-      END DO
-
-      DO jk = 2, jpkm1              ! complete with the level-dependent part
-         emix_tmx(:,:,jk) = zfact(:,:) * (  EXP( ( gde3w_n(:,:,jk  ) - zhdep(:,:) ) / hcri_tmx(:,:) )                      &
-            &                             - EXP( ( gde3w_n(:,:,jk-1) - zhdep(:,:) ) / hcri_tmx(:,:) )  ) * wmask(:,:,jk)   &
-!!gm delta(gde3w_n) = e3t_n  !!  Please verify the grid-point position w versus t-point
-            &                          / ( gde3w_n(:,:,jk) - gde3w_n(:,:,jk-1) )
-      END DO
-
-      !                        !* Pycnocline-intensified mixing: distribute energy over the time-varying 
-      !                        !* ocean depth as proportional to sqrt(rn2)^nn_zpyc
-
-      SELECT CASE ( nn_zpyc )
-
-      CASE ( 1 )               ! Dissipation scales as N (recommended)
-
-         zfact(:,:) = 0._wp
-         DO jk = 2, jpkm1              ! part independent of the level
-            zfact(:,:) = zfact(:,:) + e3w_n(:,:,jk) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk)
-         END DO
-
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_tmx(ji,jj) / ( rau0 * zfact(ji,jj) )
-            END DO
-         END DO
-
-         DO jk = 2, jpkm1              ! complete with the level-dependent part
-            emix_tmx(:,:,jk) = emix_tmx(:,:,jk) + zfact(:,:) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk)
-         END DO
-
-      CASE ( 2 )               ! Dissipation scales as N^2
-
-         zfact(:,:) = 0._wp
-         DO jk = 2, jpkm1              ! part independent of the level
-            zfact(:,:) = zfact(:,:) + e3w_n(:,:,jk) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk)
-         END DO
-
-         DO jj= 1, jpj
-            DO ji = 1, jpi
-               IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_tmx(ji,jj) / ( rau0 * zfact(ji,jj) )
-            END DO
-         END DO
-
-         DO jk = 2, jpkm1              ! complete with the level-dependent part
-            emix_tmx(:,:,jk) = emix_tmx(:,:,jk) + zfact(:,:) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk)
-         END DO
-
-      END SELECT
-
-      !                        !* WKB-height dependent mixing: distribute energy over the time-varying 
-      !                        !* ocean depth as proportional to rn2 * exp(-z_wkb/rn_hbot)
-      
-      zwkb(:,:,:) = 0._wp
-      zfact(:,:) = 0._wp
-      DO jk = 2, jpkm1
-         zfact(:,:) = zfact(:,:) + e3w_n(:,:,jk) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk)
-         zwkb(:,:,jk) = zfact(:,:)
-      END DO
-
-      DO jk = 2, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF( zfact(ji,jj) /= 0 )   zwkb(ji,jj,jk) = zhdep(ji,jj) * ( zfact(ji,jj) - zwkb(ji,jj,jk) )   &
-                                            &           * tmask(ji,jj,jk) / zfact(ji,jj)
-            END DO
-         END DO
-      END DO
-      zwkb(:,:,1) = zhdep(:,:) * tmask(:,:,1)
-
-      zweight(:,:,:) = 0._wp
-      DO jk = 2, jpkm1
-         zweight(:,:,jk) = MAX( 0._wp, rn2(:,:,jk) ) * hbot_tmx(:,:) * wmask(:,:,jk)                    &
-            &   * (  EXP( -zwkb(:,:,jk) / hbot_tmx(:,:) ) - EXP( -zwkb(:,:,jk-1) / hbot_tmx(:,:) )  )
-      END DO
-
-      zfact(:,:) = 0._wp
-      DO jk = 2, jpkm1              ! part independent of the level
-         zfact(:,:) = zfact(:,:) + zweight(:,:,jk)
-      END DO
-
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ebot_tmx(ji,jj) / ( rau0 * zfact(ji,jj) )
-         END DO
-      END DO
-
-      DO jk = 2, jpkm1              ! complete with the level-dependent part
-         emix_tmx(:,:,jk) = emix_tmx(:,:,jk) + zweight(:,:,jk) * zfact(:,:) * wmask(:,:,jk)   &
-            &                                / ( gde3w_n(:,:,jk) - gde3w_n(:,:,jk-1) )
-      END DO
-
-
-      ! Calculate molecular kinematic viscosity
-      znu_t(:,:,:) = 1.e-4_wp * (  17.91_wp - 0.53810_wp * tsn(:,:,:,jp_tem) + 0.00694_wp * tsn(:,:,:,jp_tem) * tsn(:,:,:,jp_tem)  &
-         &                                  + 0.02305_wp * tsn(:,:,:,jp_sal)  ) * tmask(:,:,:) * r1_rau0
-      DO jk = 2, jpkm1
-         znu_w(:,:,jk) = 0.5_wp * ( znu_t(:,:,jk-1) + znu_t(:,:,jk) ) * wmask(:,:,jk)
-      END DO
-
-      ! Calculate turbulence intensity parameter Reb
-      DO jk = 2, jpkm1
-         zReb(:,:,jk) = emix_tmx(:,:,jk) / MAX( 1.e-20_wp, znu_w(:,:,jk) * rn2(:,:,jk) )
-      END DO
-
-      ! Define internal wave-induced diffusivity
-      DO jk = 2, jpkm1
-         zav_wave(:,:,jk) = znu_w(:,:,jk) * zReb(:,:,jk) * r1_6   ! This corresponds to a constant mixing efficiency of 1/6
-      END DO
-
-      IF( ln_mevar ) THEN              ! Variable mixing efficiency case : modify zav_wave in the
-         DO jk = 2, jpkm1              ! energetic (Reb > 480) and buoyancy-controlled (Reb <10.224 ) regimes
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  IF( zReb(ji,jj,jk) > 480.00_wp ) THEN
-                     zav_wave(ji,jj,jk) = 3.6515_wp * znu_w(ji,jj,jk) * SQRT( zReb(ji,jj,jk) )
-                  ELSEIF( zReb(ji,jj,jk) < 10.224_wp ) THEN
-                     zav_wave(ji,jj,jk) = 0.052125_wp * znu_w(ji,jj,jk) * zReb(ji,jj,jk) * SQRT( zReb(ji,jj,jk) )
-                  ENDIF
-               END DO
-            END DO
-         END DO
-      ENDIF
-
-      DO jk = 2, jpkm1                 ! Bound diffusivity by molecular value and 100 cm2/s
-         zav_wave(:,:,jk) = MIN(  MAX( 1.4e-7_wp, zav_wave(:,:,jk) ), 1.e-2_wp  ) * wmask(:,:,jk)
-      END DO
-
-      IF( kt == nit000 ) THEN        !* Control print at first time-step: diagnose the energy consumed by zav_wave
-         ztpc = 0._wp
-!!gm used of glosum 3D....
-         DO jk = 2, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ztpc = ztpc + e3w_n(ji,jj,jk) * e1e2t(ji,jj)   &
-                     &         * MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj)
-               END DO
-            END DO
-         END DO
-         IF( lk_mpp )   CALL mpp_sum( ztpc )
-         ztpc = rau0 * ztpc ! Global integral of rauo * Kz * N^2 = power contributing to mixing 
- 
-         IF(lwp) THEN
-            WRITE(numout,*)
-            WRITE(numout,*) 'zdf_tmx : Internal wave-driven mixing (tmx)'
-            WRITE(numout,*) '~~~~~~~ '
-            WRITE(numout,*)
-            WRITE(numout,*) '      Total power consumption by av_wave: ztpc =  ', ztpc * 1.e-12_wp, 'TW'
-         ENDIF
-      ENDIF
-
-      !                          ! ----------------------- !
-      !                          !   Update  mixing coefs  !                          
-      !                          ! ----------------------- !
-      !      
-      IF( ln_tsdiff ) THEN          !* Option for differential mixing of salinity and temperature
-         DO jk = 2, jpkm1              ! Calculate S/T diffusivity ratio as a function of Reb
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zav_ratio(ji,jj,jk) = ( 0.505_wp + 0.495_wp *                                                                  &
-                      &   TANH(    0.92_wp * (   LOG10(  MAX( 1.e-20_wp, zReb(ji,jj,jk) * 5._wp * r1_6 )  ) - 0.60_wp   )    )   &
-                      &                 ) * wmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-         CALL iom_put( "av_ratio", zav_ratio )
-         DO jk = 2, jpkm1           !* update momentum & tracer diffusivity with wave-driven mixing
-            avs(:,:,jk) = avs(:,:,jk) + zav_wave(:,:,jk) * zav_ratio(:,:,jk)
-            avt(:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk)
-            avm(:,:,jk) = avm(:,:,jk) + zav_wave(:,:,jk)
-         END DO
-         !
-      ELSE                          !* update momentum & tracer diffusivity with wave-driven mixing
-         DO jk = 2, jpkm1
-            avs(:,:,jk) = avs(:,:,jk) + zav_wave(:,:,jk)
-            avt(:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk)
-            avm(:,:,jk) = avm(:,:,jk) + zav_wave(:,:,jk)
-         END DO
-      ENDIF
-
-      DO jk = 2, jpkm1              !* update momentum diffusivity at wu and wv points
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1  ! vector opt.
-               avmu(ji,jj,jk) = avmu(ji,jj,jk) + 0.5_wp * ( zav_wave(ji,jj,jk) + zav_wave(ji+1,jj  ,jk) ) * wumask(ji,jj,jk)
-               avmv(ji,jj,jk) = avmv(ji,jj,jk) + 0.5_wp * ( zav_wave(ji,jj,jk) + zav_wave(ji  ,jj+1,jk) ) * wvmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! lateral boundary condition
-
-      !                             !* output internal wave-driven mixing coefficient
-      CALL iom_put( "av_wave", zav_wave )
-                                    !* output useful diagnostics: N^2, Kz * N^2 (bflx_tmx), 
-                                    !  vertical integral of rau0 * Kz * N^2 (pcmap_tmx), energy density (emix_tmx)
-      IF( iom_use("bflx_tmx") .OR. iom_use("pcmap_tmx") ) THEN
-         bflx_tmx(:,:,:) = MAX( 0._wp, rn2(:,:,:) ) * zav_wave(:,:,:)
-         pcmap_tmx(:,:) = 0._wp
-         DO jk = 2, jpkm1
-            pcmap_tmx(:,:) = pcmap_tmx(:,:) + e3w_n(:,:,jk) * bflx_tmx(:,:,jk) * wmask(:,:,jk)
-         END DO
-         pcmap_tmx(:,:) = rau0 * pcmap_tmx(:,:)
-         CALL iom_put( "bflx_tmx", bflx_tmx )
-         CALL iom_put( "pcmap_tmx", pcmap_tmx )
-      ENDIF
-      CALL iom_put( "bn2", rn2 )
-      CALL iom_put( "emix_tmx", emix_tmx )
-      
-      CALL wrk_dealloc( jpi,jpj,       zfact, zhdep )
-      CALL wrk_dealloc( jpi,jpj,jpk,   zwkb, zweight, znu_t, znu_w, zReb )
-
-      IF(ln_ctl)   CALL prt_ctl(tab3d_1=zav_wave , clinfo1=' tmx - av_wave: ', tab3d_2=avt, clinfo2=' avt: ', ovlap=1, kdim=jpk)
-      !
-      IF( nn_timing == 1 )   CALL timing_stop('zdf_tmx')
-      !
-   END SUBROUTINE zdf_tmx
-
-
-   SUBROUTINE zdf_tmx_init
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE zdf_tmx_init  ***
-      !!                     
-      !! ** Purpose :   Initialization of the wave-driven vertical mixing, reading
-      !!              of input power maps and decay length scales in netcdf files.
-      !!
-      !! ** Method  : - Read the namzdf_tmx namelist and check the parameters
-      !!
-      !!              - Read the input data in NetCDF files :
-      !!              power available from high-mode wave breaking (mixing_power_bot.nc)
-      !!              power available from pycnocline-intensified wave-breaking (mixing_power_pyc.nc)
-      !!              power available from critical slope wave-breaking (mixing_power_cri.nc)
-      !!              WKB decay scale for high-mode wave-breaking (decay_scale_bot.nc)
-      !!              decay scale for critical slope wave-breaking (decay_scale_cri.nc)
-      !!
-      !! ** input   : - Namlist namzdf_tmx
-      !!              - NetCDF files : mixing_power_bot.nc, mixing_power_pyc.nc, mixing_power_cri.nc,
-      !!              decay_scale_bot.nc decay_scale_cri.nc
-      !!
-      !! ** Action  : - Increase by 1 the nstop flag is setting problem encounter
-      !!              - Define ebot_tmx, epyc_tmx, ecri_tmx, hbot_tmx, hcri_tmx
-      !!
-      !! References : de Lavergne et al. 2015, JPO; 2016, in prep.
-      !!			  
-      !!----------------------------------------------------------------------
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      INTEGER  ::   inum         ! local integer
-      INTEGER  ::   ios
-      REAL(wp) ::   zbot, zpyc, zcri   ! local scalars
-      !!
-      NAMELIST/namzdf_tmx_new/ nn_zpyc, ln_mevar, ln_tsdiff
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('zdf_tmx_init')
-      !
-      REWIND( numnam_ref )              ! Namelist namzdf_tmx in reference namelist : Wave-driven mixing
-      READ  ( numnam_ref, namzdf_tmx_new, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tmx in reference namelist', lwp )
-      !
-      REWIND( numnam_cfg )              ! Namelist namzdf_tmx in configuration namelist : Wave-driven mixing
-      READ  ( numnam_cfg, namzdf_tmx_new, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tmx in configuration namelist', lwp )
-      IF(lwm) WRITE ( numond, namzdf_tmx_new )
-      !
-      IF(lwp) THEN                  ! Control print
-         WRITE(numout,*)
-         WRITE(numout,*) 'zdf_tmx_init : internal wave-driven mixing'
-         WRITE(numout,*) '~~~~~~~~~~~~'
-         WRITE(numout,*) '   Namelist namzdf_tmx_new : set wave-driven mixing parameters'
-         WRITE(numout,*) '      Pycnocline-intensified diss. scales as N (=1) or N^2 (=2) = ', nn_zpyc
-         WRITE(numout,*) '      Variable (T) or constant (F) mixing efficiency            = ', ln_mevar
-         WRITE(numout,*) '      Differential internal wave-driven mixing (T) or not (F)   = ', ln_tsdiff
-      ENDIF
-      
-      ! The new wave-driven mixing parameterization elevates avt and avm in the interior, and
-      ! ensures that avt remains larger than its molecular value (=1.4e-7). Therefore, avtb should 
-      ! be set here to a very small value, and avmb to its (uniform) molecular value (=1.4e-6).
-      avmb(:) = 1.4e-6_wp        ! viscous molecular value
-      avtb(:) = 1.e-10_wp        ! very small diffusive minimum (background avt is specified in zdf_tmx)    
-      avtb_2d(:,:) = 1.e0_wp     ! uniform 
-      IF(lwp) THEN                  ! Control print
-         WRITE(numout,*)
-         WRITE(numout,*) '   Force the background value applied to avm & avt in TKE to be everywhere ',   &
-            &               'the viscous molecular value & a very small diffusive value, resp.'
-      ENDIF
-            
-      !                             ! allocate tmx arrays
-      IF( zdf_tmx_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_tmx_init : unable to allocate tmx arrays' )
-      !
-      !                             ! read necessary fields
-      CALL iom_open('mixing_power_bot',inum)       ! energy flux for high-mode wave breaking [W/m2]
-      CALL iom_get  (inum, jpdom_data, 'field', ebot_tmx, 1 ) 
-      CALL iom_close(inum)
-      !
-      CALL iom_open('mixing_power_pyc',inum)       ! energy flux for pynocline-intensified wave breaking [W/m2]
-      CALL iom_get  (inum, jpdom_data, 'field', epyc_tmx, 1 )
-      CALL iom_close(inum)
-      !
-      CALL iom_open('mixing_power_cri',inum)       ! energy flux for critical slope wave breaking [W/m2]
-      CALL iom_get  (inum, jpdom_data, 'field', ecri_tmx, 1 )
-      CALL iom_close(inum)
-      !
-      CALL iom_open('decay_scale_bot',inum)        ! spatially variable decay scale for high-mode wave breaking [m]
-      CALL iom_get  (inum, jpdom_data, 'field', hbot_tmx, 1 )
-      CALL iom_close(inum)
-      !
-      CALL iom_open('decay_scale_cri',inum)        ! spatially variable decay scale for critical slope wave breaking [m]
-      CALL iom_get  (inum, jpdom_data, 'field', hcri_tmx, 1 )
-      CALL iom_close(inum)
-
-      ebot_tmx(:,:) = ebot_tmx(:,:) * ssmask(:,:)
-      epyc_tmx(:,:) = epyc_tmx(:,:) * ssmask(:,:)
-      ecri_tmx(:,:) = ecri_tmx(:,:) * ssmask(:,:)
-
-      ! Set once for all to zero the first and last vertical levels of appropriate variables
-      emix_tmx (:,:, 1 ) = 0._wp
-      emix_tmx (:,:,jpk) = 0._wp
-      zav_ratio(:,:, 1 ) = 0._wp
-      zav_ratio(:,:,jpk) = 0._wp
-      zav_wave (:,:, 1 ) = 0._wp
-      zav_wave (:,:,jpk) = 0._wp
-
-      zbot = glob_sum( e1e2t(:,:) * ebot_tmx(:,:) )
-      zpyc = glob_sum( e1e2t(:,:) * epyc_tmx(:,:) )
-      zcri = glob_sum( e1e2t(:,:) * ecri_tmx(:,:) )
-      IF(lwp) THEN
-         WRITE(numout,*) '      High-mode wave-breaking energy:             ', zbot * 1.e-12_wp, 'TW'
-         WRITE(numout,*) '      Pycnocline-intensifed wave-breaking energy: ', zpyc * 1.e-12_wp, 'TW'
-         WRITE(numout,*) '      Critical slope wave-breaking energy:        ', zcri * 1.e-12_wp, 'TW'
-      ENDIF
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('zdf_tmx_init')
-      !
-   END SUBROUTINE zdf_tmx_init
-
-   !!======================================================================
-END MODULE zdftmx
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/module_example
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/module_example	(revision 7989)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/module_example	(revision 7990)
@@ -86,9 +86,4 @@
       !!                Give references if exist otherwise suppress these lines
       !!----------------------------------------------------------------------
-      USE toto_module      ! description of the module
-      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
-      USE wrk_nemo, ONLY:   zztab => wrk_2d_5                     ! 2D workspace
-      USE wrk_nemo, ONLY:   zwx => wrk_3d_12 , zwy => wrk_3d_13   ! 3D workspace
-      !!
       INTEGER , INTENT(in   )                     ::   kt      ! short description 
       INTEGER , INTENT(inout)                     ::   pvar1   !   -         -
@@ -100,10 +95,7 @@
       REAL(wp) ::   zmlmin, zbbrau   ! temporary scalars     (DOCTOR : start with z)
       REAL(wp) ::   zfact1, zfact2   ! do not use continuation lines in declaration
+      REAL(wp), DIMENSION(jpi,jpj) ::   zwrk_2d   ! 2D workspace
       !!--------------------------------------------------------------------
-
-      IF( wrk_in_use(3, 12,13) .OR. wrk_in_use(2, 5 ) THEN
-         CALL ctl_stop('exa_mpl: requested workspace arrays unavailable')   ;   RETURN
-      ENDIF
-
+      !
       IF( kt == nit000  )   CALL exa_mpl_init    ! Initialization (first time-step only)
 
@@ -119,5 +111,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  avmv(ji,jj,jk) = ....
+                  avm(ji,jj,jk) = ....
                END DO
             END DO
@@ -128,5 +120,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  avmv(ji,jj,jk) = ...
+                  avm(ji,jj,jk) = ...
                END DO
             END DO
@@ -135,9 +127,5 @@
       END SELECT
       !
-      CALL mpplnk2( avmu, 'U', 1. )              ! Lateral boundary conditions (unchanged sign)
-      !
-      IF( wrk_not_released(3, 12,13) .OR. wrk_not_released(2, 5 ) THEN
-         CALL ctl_stop('exa_mpl: failed to release workspace arrays')   ;   RETURN
-      ENDIF
+      CALL lbc_lnk( avm, 'T', 1. )              ! Lateral boundary conditions (unchanged sign)
       !
    END SUBROUTINE exa_mpl
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 7989)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 7990)
@@ -64,15 +64,4 @@
 
    USE zdfphy         ! vertical physics manager      (zdf_phy_init routine)
-!!gm to be suppressed
-   USE zdftmx           ! tide-induced vertical mixing     (zdf_tmx routine)
-   USE zdfbfr           ! bottom friction                  (zdf_bfr routine)
-   USE zdftke           ! TKE vertical mixing              (zdf_tke routine)
-   USE zdfgls           ! GLS vertical mixing              (zdf_gls routine)
-   USE zdfddm           ! double diffusion mixing          (zdf_ddm routine)
-   USE zdfevd           ! enhanced vertical diffusion      (zdf_evd routine)
-   USE zdfric           ! Richardson vertical mixing       (zdf_ric routine)
-   USE zdfmxl           ! Mixed-layer depth                (zdf_mxl routine)
-   USE zdfqiao          !Qiao module wave induced mixing   (zdf_qiao routine)
-!!gm end
 
    USE step_diu        ! Time stepping for diurnal sst
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zopt.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zopt.F90	(revision 7989)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zopt.F90	(revision 7990)
@@ -119,5 +119,5 @@
       !                                          ! --------------
       neln(:,:) = 1                                   ! euphotic layer level
-      DO jk = 1, jpk                                  ! (i.e. 1rst T-level strictly below EL bottom)
+      DO jk = 1, jpkm1                                ! (i.e. 1rst T-level strictly below EL bottom)
          DO jj = 1, jpj
            DO ji = 1, jpi
@@ -147,4 +147,5 @@
    END SUBROUTINE p2z_opt
 
+
    SUBROUTINE p2z_opt_init
       !!----------------------------------------------------------------------
