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 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/AMM12/EXP00/namelist_cfg	(revision 7953)
@@ -29,9 +29,4 @@
 !-----------------------------------------------------------------------
    rn_rdt      =   600.    !  time step for the dynamics (and tracer if nn_acc=0)
-/
-!-----------------------------------------------------------------------
-&namcrs        !   Grid coarsening for dynamics output and/or
-!              !   passive tracer coarsened online simulations
-!-----------------------------------------------------------------------
 /
 !-----------------------------------------------------------------------
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 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/GYRE_PISCES/EXP00/namelist_cfg	(revision 7953)
@@ -12,7 +12,4 @@
 !!                   ***  Run management namelists  ***
 !!======================================================================
-!!   namrun        parameters of the run
-!!======================================================================
-!
 !-----------------------------------------------------------------------
 &namrun        !   parameters of the run
@@ -25,4 +22,7 @@
    nn_write    =      60   !  frequency of write in the output file   (modulo referenced to nn_it000)
 /
+!!======================================================================
+!!                      ***  Domain namelists  ***
+!!======================================================================
 !-----------------------------------------------------------------------
 &namcfg     !   parameters of the configuration   
@@ -37,13 +37,6 @@
 !-----------------------------------------------------------------------
    ln_linssh   = .true.    !  =T  linear free surface  ==>>  model level are fixed in time
-   !
    nn_msh      =    0      !  create (>0) a mesh file or not (=0)
-   !
    rn_rdt      = 7200.     !  time step for the dynamics (and tracer if nn_acc=0)
-/
-!-----------------------------------------------------------------------
-&namcrs        !   Grid coarsening for dynamics output and/or
-!              !   passive tracer coarsened online simulations
-!-----------------------------------------------------------------------
 /
 !-----------------------------------------------------------------------
@@ -56,4 +49,8 @@
    ln_tsd_tradmp = .false.   !  damping of ocean T & S toward T &S input data (T) or not (F)
 /
+
+!!======================================================================
+!!            ***  Surface Boundary Condition namelists  ***
+!!======================================================================
 !-----------------------------------------------------------------------
 &namsbc        !   Surface Boundary Condition (surface module)
@@ -76,25 +73,4 @@
 /
 !-----------------------------------------------------------------------
-&namsbc_rnf    !   runoffs namelist surface boundary condition
-!-----------------------------------------------------------------------
-   ln_rnf_mouth = .false.   !  specific treatment at rivers mouths
-/
-!-----------------------------------------------------------------------
-&namsbc_apr    !   Atmospheric pressure used as ocean forcing or in bulk
-!-----------------------------------------------------------------------
-/
-!-----------------------------------------------------------------------
-&namsbc_ssr    !   surface boundary condition : sea surface restoring
-!-----------------------------------------------------------------------
-/
-!-----------------------------------------------------------------------
-&namsbc_alb    !   albedo parameters
-!-----------------------------------------------------------------------
-/
-!-----------------------------------------------------------------------
-&namberg       !   iceberg parameters
-!-----------------------------------------------------------------------
-/
-!-----------------------------------------------------------------------
 &namlbc        !   lateral momentum boundary condition
 !-----------------------------------------------------------------------
@@ -102,31 +78,7 @@
 /
 !-----------------------------------------------------------------------
-&namagrif      !  AGRIF zoom                                            ("key_agrif")
-!-----------------------------------------------------------------------
-/
-!-----------------------------------------------------------------------
-&nam_tide      !    tide parameters
-!-----------------------------------------------------------------------
-/
-!-----------------------------------------------------------------------
-&nambdy        !  unstructured open boundaries                          
-!-----------------------------------------------------------------------
-/
-!-----------------------------------------------------------------------
-&nambdy_dta      !  open boundaries - external data           
-!-----------------------------------------------------------------------
-/
-!-----------------------------------------------------------------------
-&nambdy_tide     ! tidal forcing at open boundaries
-!-----------------------------------------------------------------------
-/
-!-----------------------------------------------------------------------
 &nambfr        !   bottom friction
 !-----------------------------------------------------------------------
    nn_bfr      =    2      !  type of bottom friction :   = 0 : free slip,  = 1 : linear friction
-/
-!-----------------------------------------------------------------------
-&nambbc        !   bottom temperature boundary condition                (default: NO)
-!-----------------------------------------------------------------------
 /
 !-----------------------------------------------------------------------
@@ -282,22 +234,22 @@
 /
 !-----------------------------------------------------------------------
-&namzdf_ric    !   richardson number dependent vertical diffusion       ("key_zdfric" )
-!-----------------------------------------------------------------------
-/
-!-----------------------------------------------------------------------
-&namzdf_tke    !   turbulent eddy kinetic dependent vertical diffusion  ("key_zdftke")
+&namzdf_ric    !   richardson number dependent vertical diffusion       (ln_zdfric=T)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namzdf_tke    !   turbulent eddy kinetic dependent vertical diffusion  (ln_zdftke=T)
 !-----------------------------------------------------------------------
    nn_etau     =   0       !  penetration of tke below the mixed layer (ML) due to internal & intertial waves
 /
 !-----------------------------------------------------------------------
-&namzdf_gls                !   GLS vertical diffusion                   ("key_zdfgls")
-!-----------------------------------------------------------------------
-/
-!-----------------------------------------------------------------------
-&namzdf_ddm    !   double diffusive mixing parameterization             ("key_zdfddm")
-!-----------------------------------------------------------------------
-/
-!-----------------------------------------------------------------------
-&namzdf_tmx    !   tidal mixing parameterization                        ("key_zdftmx")
+&namzdf_gls                !   GLS vertical diffusion                   (ln_zdfgls=T)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namzdf_ddm    !   double diffusive mixing parameterization             (ln_zdfddm=T)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namzdf_tmx    !   tidal mixing parameterization                        (ln_zdftmx=T)
 !-----------------------------------------------------------------------
    ln_tmx_itf  = .false.   !  ITF specific parameterisation
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/GYRE_PISCES/EXP00/namelist_top_cfg
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/GYRE_PISCES/EXP00/namelist_top_cfg	(revision 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/GYRE_PISCES/EXP00/namelist_top_cfg	(revision 7953)
@@ -43,8 +43,4 @@
 /
 !-----------------------------------------------------------------------
-&namtrc_zdf        !   vertical physics
-!-----------------------------------------------------------------------
-/
-!-----------------------------------------------------------------------
 &namtrc_rad        !  treatment of negative concentrations
 !-----------------------------------------------------------------------
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/GYRE_PISCES/cpp_GYRE_PISCES.fcm
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/GYRE_PISCES/cpp_GYRE_PISCES.fcm	(revision 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/GYRE_PISCES/cpp_GYRE_PISCES.fcm	(revision 7953)
@@ -1,1 +1,1 @@
-bld::tool::fppkeys key_zdftke key_top key_mpp_mpi
+bld::tool::fppkeys key_zdftke key_top key_mpp_mpi key_nosignedzero
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 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/ORCA2_LIM3_PISCES/EXP00/namelist_cfg	(revision 7953)
@@ -235,13 +235,37 @@
 /
 !-----------------------------------------------------------------------
-&namzdf        !   vertical physics
-!-----------------------------------------------------------------------
-/
-!-----------------------------------------------------------------------
-&namzdf_tke    !   turbulent eddy kinetic dependent vertical diffusion  ("key_zdftke")
-!-----------------------------------------------------------------------
-/
-!-----------------------------------------------------------------------
-&namzdf_ddm    !   double diffusive mixing parameterization             ("key_zdfddm")
+&namzdf        !   vertical physics                                     (default: NO selection)
+!-----------------------------------------------------------------------
+   !                       ! type of vertical closure
+   ln_zdfcst   = .false.      !  constant mixing
+   ln_zdfric   = .false.      !  local Richardson dependent formulation (T =>   fill namzdf_ric)
+   ln_zdftke   = .true.       !  Turbulent Kinetic Energy closure       (T =>   fill namzdf_tke)
+   ln_zdfgls   = .false.      !  Generic Length Scale closure           (T =>   fill namzdf_gls)
+   !
+   !                       ! convection
+   ln_zdfevd   = .true.       !  Enhanced Vertical Diffusion scheme
+      nn_evdm  =    0            !  evd apply on tracer (=0) or on tracer and momentum (=1)
+      rn_evd   =  100.           !  evd mixing coefficient [m2/s]
+   !
+   ln_zdfddm   = .true.    ! double diffusive mixing
+      rn_avts  =    1.e-4     !  maximum avs (vertical mixing on salinity)
+      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)
+   !
+   !                       ! 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
+   !
+   !                       !  Coefficients
+   rn_avm0     =   1.2e-4     !  vertical eddy viscosity   [m2/s]       (background Kz if ln_zdfcst=F)
+   rn_avt0     =   1.2e-5     !  vertical eddy diffusivity [m2/s]       (background Kz if ln_zdfcst=F)
+   nn_avb      =    0         !  profile for background avt & avm (=1) or not (=0)
+   nn_havtb    =    1         !  horizontal shape for avtb (=1) or not (=0)
+/
+!-----------------------------------------------------------------------
+&namzdf_tke    !   turbulent eddy kinetic dependent vertical diffusion 
 !-----------------------------------------------------------------------
 /
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/ORCA2_LIM3_PISCES/EXP00/namelist_top_cfg
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/ORCA2_LIM3_PISCES/EXP00/namelist_top_cfg	(revision 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/ORCA2_LIM3_PISCES/EXP00/namelist_top_cfg	(revision 7953)
@@ -81,8 +81,4 @@
 /
 !-----------------------------------------------------------------------
-&namtrc_zdf        !   vertical physics
-!-----------------------------------------------------------------------
-/
-!-----------------------------------------------------------------------
 &namtrc_rad        !  treatment of negative concentrations 
 !-----------------------------------------------------------------------
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/ORCA2_LIM3_PISCES/cpp_ORCA2_LIM3_PISCES.fcm
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/ORCA2_LIM3_PISCES/cpp_ORCA2_LIM3_PISCES.fcm	(revision 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/ORCA2_LIM3_PISCES/cpp_ORCA2_LIM3_PISCES.fcm	(revision 7953)
@@ -1,1 +1,1 @@
-bld::tool::fppkeys key_trabbl key_lim3 key_zdftke key_zdfddm key_zdftmx_new key_iomput key_mpp_mpi key_top key_nosignedzero
+bld::tool::fppkeys key_trabbl key_lim3 key_zdftke key_zdftmx_new key_iomput key_mpp_mpi key_top key_nosignedzero
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 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/SHARED/field_def_nemo-opa.xml	(revision 7953)
@@ -390,5 +390,5 @@
          <field id="avm"          long_name="vertical eddy viscosity"     standard_name="ocean_vertical_momentum_diffusivity"   unit="m2/s" />
 
-         <!-- avs: available with key_zdfddm -->
+         <!-- avs: if ln_zdfddm=F avs=avt -->
          <field id="avs"          long_name="salt vertical eddy diffusivity"   standard_name="ocean_vertical_salt_diffusivity"   unit="m2/s" />
          <field id="logavs"       long_name="logarithm of salt vertical eddy diffusivity"   standard_name="ocean_vertical_heat_diffusivity"       unit="m2/s" />
@@ -398,8 +398,5 @@
          <field id="avm_evd"      long_name="convective enhancement of vertical viscosity"     standard_name="ocean_vertical_momentum_diffusivity_due_to_convection"   unit="m2/s" />
 
-         <!-- avt_tide: available with key_zdftmx -->
-         <field id="av_tide"      long_name="tidal vertical diffusivity"   standard_name="ocean_vertical_tracer_diffusivity_due_to_tides"   unit="m2/s" />
-
-         <!-- variables available with key_zdftmx_new -->
+         <!-- variables available with ln_zdftmx =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" />
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 7953)
@@ -236,5 +236,5 @@
 /
 !-----------------------------------------------------------------------
-&namsbc_blk   !   namsbc_blk  generic Bulk formula                      (ln_blk = T)
+&namsbc_blk   !   namsbc_blk  generic Bulk formula                      (ln_blk =T)
 !-----------------------------------------------------------------------
 !              !  file name                   ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights                              ! rotation ! land/sea mask !
@@ -326,5 +326,5 @@
 /
 !-----------------------------------------------------------------------
-&namtra_qsr    !   penetrative solar radiation                          (ln_traqsr=T)
+&namtra_qsr    !   penetrative solar radiation                          (ln_traqsr =T)
 !-----------------------------------------------------------------------
 !              !  file name  ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
@@ -343,5 +343,5 @@
 /
 !-----------------------------------------------------------------------
-&namsbc_rnf    !   runoffs namelist surface boundary condition          (ln_rnf=T)
+&namsbc_rnf    !   runoffs namelist surface boundary condition          (ln_rnf =T)
 !-----------------------------------------------------------------------
 !              !  file name           ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
@@ -419,5 +419,5 @@
 /
 !-----------------------------------------------------------------------
-&namsbc_ssr    !   surface boundary condition : sea surface restoring   (ln_ssr=T)
+&namsbc_ssr    !   surface boundary condition : sea surface restoring   (ln_ssr =T)
 !-----------------------------------------------------------------------
 !              ! file name ! frequency (hours) ! variable ! time interp.!  clim  ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
@@ -448,5 +448,5 @@
 /
 !-----------------------------------------------------------------------
-&namsbc_wave   ! External fields from wave model                        (ln_wave=T)
+&namsbc_wave   ! External fields from wave model                        (ln_wave =T)
 !-----------------------------------------------------------------------
 !              !  file name  ! frequency (hours) ! variable     ! time interp. !  clim   ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
@@ -667,5 +667,5 @@
 !
 !-----------------------------------------------------------------------
-&nameos        !   ocean Equation Of State                              (default: NO)
+&nameos        !   ocean Equation Of Seawater                           (default: NO)
 !-----------------------------------------------------------------------
    ln_teos10   = .false.         !  = Use TEOS-10 equation of state
@@ -894,5 +894,4 @@
       nn_npc      =    1         ! frequency of application of npc
       nn_npcp     =  365         ! npc control print frequency
-!!gm new
    !
    ln_zdfddm   = .false.   ! double diffusive mixing
@@ -902,6 +901,4 @@
    ln_zdftmx   = .false.   ! tidal mixing parameterization              (T =>   fill namzdf_tmx)
    !
-!!gm
-   !
    ln_zdfqiao  = .false.   ! surface wave-induced mixing (Qiao et al. 2010) (T =>   ln_wave=ln_sdw=T. & fill namsbc_wave)
    !
@@ -909,4 +906,5 @@
    ln_zdfexp   = .false.   ! split-explicit (T) or implicit (F) scheme
       nn_zdfexp=    3         !  number of sub-timestep for ln_zdfexp=T
+   !
    !                       ! coefficients
    rn_avm0     =   1.2e-4     !  vertical eddy viscosity   [m2/s]       (background Kz if ln_zdfcst=F)
@@ -977,15 +975,5 @@
 /
 !-----------------------------------------------------------------------
-&namzdf_tmx    !   tidal mixing parameterization                        ("key_zdftmx")
-!-----------------------------------------------------------------------
-   rn_htmx     = 500.      !  vertical decay scale for turbulence (meters)
-   rn_n2min    = 1.e-8     !  threshold of the Brunt-Vaisala frequency (s-1)
-   rn_tfe      = 0.333     !  tidal dissipation efficiency
-   rn_me       = 0.2       !  mixing efficiency
-   ln_tmx_itf  = .true.    !  ITF specific parameterisation
-   rn_tfe_itf  = 1.        !  ITF tidal dissipation efficiency
-/
-!-----------------------------------------------------------------------
-&namzdf_tmx_new !   internal wave-driven mixing parameterization        ("key_zdftmx_new" & "key_zdfddm")
+&namzdf_tmx    !    internal wave-driven mixing parameterization        (ln_zdftmx =T)
 !-----------------------------------------------------------------------
    nn_zpyc     = 1         !  pycnocline-intensified dissipation scales as N (=1) or N^2 (=2)
@@ -996,5 +984,5 @@
 !!                  ***  Miscellaneous namelists  ***
 !!======================================================================
-!!   nammpp            Massively Parallel Processing                    ("key_mpp_mpi)
+!!   nammpp            Massively Parallel Processing                    ("key_mpp_mpi")
 !!   namctl            Control prints 
 !!   namsto            Stochastic parametrization of EOS
@@ -1002,5 +990,5 @@
 !
 !-----------------------------------------------------------------------
-&nammpp        !   Massively Parallel Processing                        ("key_mpp_mpi)
+&nammpp        !   Massively Parallel Processing                        ("key_mpp_mpi")
 !-----------------------------------------------------------------------
    cn_mpi_send =  'I'      !  mpi send/recieve type   ='S', 'B', or 'I' for standard send,
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/SHARED/namelist_top_ref
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/SHARED/namelist_top_ref	(revision 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/SHARED/namelist_top_ref	(revision 7953)
@@ -95,10 +95,4 @@
 /
 !-----------------------------------------------------------------------
-&namtrc_zdf      !   vertical physics
-!-----------------------------------------------------------------------
-   ln_trczdf_exp =  .false.  !  split explicit (T) or implicit (F) time stepping
-   nn_trczdf_exp =   3       !  number of sub-timestep for ln_trczdfexp=T
-/
-!-----------------------------------------------------------------------
 &namtrc_rad      !  treatment of negative concentrations 
 !-----------------------------------------------------------------------
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/cfg.txt
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/cfg.txt	(revision 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/cfg.txt	(revision 7953)
@@ -1,8 +1,9 @@
 AMM12 OPA_SRC
 C1D_PAPA OPA_SRC
-GYRE_PISCES OPA_SRC TOP_SRC
 GYRE_BFM OPA_SRC TOP_SRC
 ORCA2_SAS_LIM3 OPA_SRC SAS_SRC LIM_SRC_3 NST_SRC
 ORCA2_OFF_PISCES OPA_SRC OFF_SRC TOP_SRC
 ORCA2_OFF_TRC OPA_SRC OFF_SRC TOP_SRC
+GYRE_PISCES OPA_SRC TOP_SRC
 ORCA2_LIM3_PISCES OPA_SRC LIM_SRC_3 TOP_SRC NST_SRC
+GYRE_PISCES_XIOS OPA_SRC TOP_SRC
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/NST_SRC/agrif_ice.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/NST_SRC/agrif_ice.F90	(revision 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/NST_SRC/agrif_ice.F90	(revision 7953)
@@ -18,20 +18,19 @@
    PUBLIC agrif_ice_alloc ! routine called by nemo_init in nemogcm.F90
 
-   INTEGER, PUBLIC :: u_ice_id, v_ice_id, adv_ice_id
-   REAL(wp), PUBLIC :: lim_nbstep = 0.    ! child time position in sea-ice model
+   INTEGER , PUBLIC ::   u_ice_id, v_ice_id, adv_ice_id
+   REAL(wp), PUBLIC ::   lim_nbstep = 0.    ! child time position in sea-ice model
 #if defined key_lim2_vp
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)     :: u_ice_nst, v_ice_nst   
 #else
    REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)   :: u_ice_oe, u_ice_sn     !: boundaries arrays
-   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)   :: v_ice_oe, v_ice_sn     !:  "          " 
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)   :: v_ice_oe, v_ice_sn     !:     -        - 
 #endif
-   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:,:) :: adv_ice_oe, adv_ice_sn !:  "          "
+   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:,:) :: adv_ice_oe, adv_ice_sn !:     -        -
 
    !!----------------------------------------------------------------------
-   !! NEMO/NST 3.3.4 , NEMO Consortium (2012)
+   !! NEMO/NST 4.0 , NEMO Consortium (2017)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS 
 
@@ -50,11 +49,11 @@
 
 #if ! defined key_lim2_vp
-      u_ice_oe(:,:,:) =  0.e0
-      v_ice_oe(:,:,:) =  0.e0
-      u_ice_sn(:,:,:) =  0.e0
-      v_ice_sn(:,:,:) =  0.e0
+      u_ice_oe(:,:,:) =  0._wp
+      v_ice_oe(:,:,:) =  0._wp
+      u_ice_sn(:,:,:) =  0._wp
+      v_ice_sn(:,:,:) =  0._wp
 #endif
-      adv_ice_oe (:,:,:,:) = 0.e0 
-      adv_ice_sn (:,:,:,:) = 0.e0 
+      adv_ice_oe (:,:,:,:) = 0._wp 
+      adv_ice_sn (:,:,:,:) = 0._wp
       !
    END FUNCTION agrif_ice_alloc
@@ -71,5 +70,5 @@
 
    !!----------------------------------------------------------------------
-   !! NEMO/NST 3.6 , NEMO Consortium (2016)
+   !! NEMO/NST 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/NST_SRC/agrif_lim3_interp.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/NST_SRC/agrif_lim3_interp.F90	(revision 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/NST_SRC/agrif_lim3_interp.F90	(revision 7953)
@@ -28,5 +28,5 @@
    PRIVATE
 
-   PUBLIC agrif_interp_lim3
+   PUBLIC   agrif_interp_lim3   ! called by ???
 
    !!----------------------------------------------------------------------
@@ -46,13 +46,13 @@
       !!  computing factor for time interpolation
       !!-----------------------------------------------------------------------
-      CHARACTER(len=1), INTENT( in )           :: cd_type
-      INTEGER         , INTENT( in ), OPTIONAL :: kiter, kitermax
-      !!
-      REAL(wp) :: zbeta
+      CHARACTER(len=1), INTENT(in   )           ::   cd_type
+      INTEGER         , INTENT(in   ), OPTIONAL ::   kiter, kitermax
+      !!
+      REAL(wp) ::   zbeta   ! local scalar
       !!-----------------------------------------------------------------------
       !
       IF( Agrif_Root() )  RETURN
       !
-      SELECT CASE(cd_type)
+      SELECT CASE( cd_type )
       CASE('U','V')
          IF( PRESENT( kiter ) ) THEN  ! interpolation at the child sub-time step (only for ice rheology)
@@ -66,22 +66,17 @@
       END SELECT
       !
-      Agrif_SpecialValue=-9999.
+      Agrif_SpecialValue    = -9999.
       Agrif_UseSpecialValue = .TRUE.
-      SELECT CASE(cd_type)
-      CASE('U')
-         CALL Agrif_Bc_variable( u_ice_id  , procname=interp_u_ice  , calledweight=zbeta )
-      CASE('V')
-         CALL Agrif_Bc_variable( v_ice_id  , procname=interp_v_ice  , calledweight=zbeta )
-      CASE('T')
-         CALL Agrif_Bc_variable( tra_ice_id, procname=interp_tra_ice, calledweight=zbeta )
+      SELECT CASE( cd_type )
+      CASE('U')   ;   CALL Agrif_Bc_variable( u_ice_id  , procname=interp_u_ice  , calledweight=zbeta )
+      CASE('V')   ;   CALL Agrif_Bc_variable( v_ice_id  , procname=interp_v_ice  , calledweight=zbeta )
+      CASE('T')   ;   CALL Agrif_Bc_variable( tra_ice_id, procname=interp_tra_ice, calledweight=zbeta )
       END SELECT
-      Agrif_SpecialValue=0.
+      Agrif_SpecialValue    = 0._wp
       Agrif_UseSpecialValue = .FALSE.
       !
    END SUBROUTINE agrif_interp_lim3
 
-   !!------------------
-   !! Local subroutines
-   !!------------------
+
    SUBROUTINE interp_u_ice( ptab, i1, i2, j1, j2, before )
       !!-----------------------------------------------------------------------
@@ -92,9 +87,9 @@
       !! put -9999 WHERE the parent grid is masked. The child solution will be found in the 9(?) points around
       !!-----------------------------------------------------------------------
-      INTEGER , INTENT(in) :: i1, i2, j1, j2
-      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
-      LOGICAL , INTENT(in) :: before
-      !!
-      REAL(wp) :: zrhoy
+      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
+      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
+      LOGICAL                         , INTENT(in   ) ::   before
+      !!
+      REAL(wp) ::   zrhoy   ! local scalar
       !!-----------------------------------------------------------------------
       !
@@ -118,9 +113,9 @@
       !! put -9999 WHERE the parent grid is masked. The child solution will be found in the 9(?) points around
       !!-----------------------------------------------------------------------      
-      INTEGER , INTENT(in) :: i1, i2, j1, j2
-      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
-      LOGICAL , INTENT(in) :: before
-      !!
-      REAL(wp) :: zrhox
+      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
+      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
+      LOGICAL                         , INTENT(in   ) ::   before
+      !!
+      REAL(wp) ::   zrhox   ! local scalar
       !!-----------------------------------------------------------------------
       !
@@ -144,17 +139,16 @@
       !! put -9999 WHERE the parent grid is masked. The child solution will be found in the 9(?) points around
       !!-----------------------------------------------------------------------
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
-      INTEGER , INTENT(in) :: i1, i2, j1, j2, k1, k2
-      LOGICAL , INTENT(in) :: before
-      INTEGER , INTENT(in) :: nb, ndir
-      !!
-      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztab
+      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
+      INTEGER                               , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2
+      LOGICAL                               , INTENT(in   ) ::   before
+      INTEGER                               , INTENT(in   ) ::   nb, ndir
+      !!
       INTEGER  ::   ji, jj, jk, jl, jm
       INTEGER  ::   imin, imax, jmin, jmax
+      LOGICAL  ::   western_side, eastern_side, northern_side, southern_side
       REAL(wp) ::   zrhox, z1, z2, z3, z4, z5, z6, z7
-      LOGICAL  ::   western_side, eastern_side, northern_side, southern_side
-
-      !!-----------------------------------------------------------------------
-      ! tracers are not multiplied by grid cell here => before: * e12t ; after: * r1_e12t / rhox / rhoy
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztab
+      !!-----------------------------------------------------------------------
+      ! tracers are not multiplied by grid cell here => before: * e1e2t ; after: * r1_e1e2t / rhox / rhoy
       ! and it is ok since we conserve tracers (same as in the ocean).
       ALLOCATE( ztab(SIZE(a_i_b,1),SIZE(a_i_b,2),SIZE(ptab,3)) )
@@ -163,20 +157,20 @@
          jm = 1
          DO jl = 1, jpl
-            ptab(i1:i2,j1:j2,jm) = a_i_b  (i1:i2,j1:j2,jl) ; jm = jm + 1
-            ptab(i1:i2,j1:j2,jm) = v_i_b  (i1:i2,j1:j2,jl) ; jm = jm + 1
-            ptab(i1:i2,j1:j2,jm) = v_s_b  (i1:i2,j1:j2,jl) ; jm = jm + 1
-            ptab(i1:i2,j1:j2,jm) = smv_i_b(i1:i2,j1:j2,jl) ; jm = jm + 1
-            ptab(i1:i2,j1:j2,jm) = oa_i_b (i1:i2,j1:j2,jl) ; jm = jm + 1
+            ptab(i1:i2,j1:j2,jm) = a_i_b  (i1:i2,j1:j2,jl)   ;   jm = jm + 1
+            ptab(i1:i2,j1:j2,jm) = v_i_b  (i1:i2,j1:j2,jl)   ;   jm = jm + 1
+            ptab(i1:i2,j1:j2,jm) = v_s_b  (i1:i2,j1:j2,jl)   ;   jm = jm + 1
+            ptab(i1:i2,j1:j2,jm) = smv_i_b(i1:i2,j1:j2,jl)   ;   jm = jm + 1
+            ptab(i1:i2,j1:j2,jm) = oa_i_b (i1:i2,j1:j2,jl)   ;   jm = jm + 1
             DO jk = 1, nlay_s
-               ptab(i1:i2,j1:j2,jm) = e_s_b(i1:i2,j1:j2,jk,jl) ; jm = jm + 1
-            ENDDO
+               ptab(i1:i2,j1:j2,jm) = e_s_b(i1:i2,j1:j2,jk,jl)   ;   jm = jm + 1
+            END DO
             DO jk = 1, nlay_i
-               ptab(i1:i2,j1:j2,jm) = e_i_b(i1:i2,j1:j2,jk,jl) ; jm = jm + 1
-            ENDDO
-         ENDDO
+               ptab(i1:i2,j1:j2,jm) = e_i_b(i1:i2,j1:j2,jk,jl)   ;   jm = jm + 1
+            END DO
+         END DO
          
          DO jk = k1, k2
-            WHERE( tmask(i1:i2,j1:j2,1) == 0. )  ptab(i1:i2,j1:j2,jk) = -9999.
-         ENDDO
+            WHERE( tmask(i1:i2,j1:j2,1) == 0._wp )   ptab(i1:i2,j1:j2,jk) = -9999.
+         END DO
          
       ELSE               ! child grid
@@ -184,16 +178,16 @@
          jm = 1
          DO jl = 1, jpl
-            a_i  (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
-            v_i  (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
-            v_s  (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
-            smv_i(i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
-            oa_i (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
+            a_i  (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1)   ;   jm = jm + 1
+            v_i  (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1)   ;   jm = jm + 1
+            v_s  (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1)   ;   jm = jm + 1
+            smv_i(i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1)   ;   jm = jm + 1
+            oa_i (i1:i2,j1:j2,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1)   ;   jm = jm + 1
             DO jk = 1, nlay_s
-               e_s(i1:i2,j1:j2,jk,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
-            ENDDO
+               e_s(i1:i2,j1:j2,jk,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1)   ;   jm = jm + 1
+            END DO
             DO jk = 1, nlay_i
-               e_i(i1:i2,j1:j2,jk,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
-            ENDDO
-         ENDDO
+               e_i(i1:i2,j1:j2,jk,jl) = ptab(i1:i2,j1:j2,jm) * tmask(i1:i2,j1:j2,1)   ;   jm = jm + 1
+            END DO
+         END DO
 
 !! ==> this is a more complex interpolation since we mix solutions over a couple of grid points
@@ -319,5 +313,5 @@
          et_s(i1:i2,j1:j2)  = SUM( SUM( e_s(i1:i2,j1:j2,:,:), dim=4 ), dim=3 )
          et_i(i1:i2,j1:j2)  = SUM( SUM( e_i(i1:i2,j1:j2,:,:), dim=4 ), dim=3 )
-
+         !
       ENDIF
       
@@ -327,11 +321,13 @@
 
 #else
+   !!----------------------------------------------------------------------
+   !!   Empty module                                             no sea-ice
+   !!----------------------------------------------------------------------
 CONTAINS
    SUBROUTINE agrif_lim3_interp_empty
-      !!---------------------------------------------
-      !!   *** ROUTINE agrif_lim3_interp_empty ***
-      !!---------------------------------------------
       WRITE(*,*)  'agrif_lim3_interp : You should not have seen this print! error?'
    END SUBROUTINE agrif_lim3_interp_empty
 #endif
+
+   !!======================================================================
 END MODULE agrif_lim3_interp
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/NST_SRC/agrif_lim3_update.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/NST_SRC/agrif_lim3_update.F90	(revision 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/NST_SRC/agrif_lim3_update.F90	(revision 7953)
@@ -31,12 +31,11 @@
    PRIVATE
 
-   PUBLIC agrif_update_lim3
+   PUBLIC   agrif_update_lim3   ! called by ????
 
    !!----------------------------------------------------------------------
-   !! NEMO/NST 3.6 , LOCEAN-IPSL (2016)
+   !! NEMO/NST 4.0 , LOCEAN-IPSL (2017)
    !! $Id: agrif_lim3_update.F90 6204 2016-01-04 13:47:06Z cetlod $
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -49,5 +48,4 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) :: kt
-      !!
       !!----------------------------------------------------------------------
       !
@@ -57,5 +55,5 @@
                                                                                                                            ! i.e. update only at the parent time step
       Agrif_UseSpecialValueInUpdate = .TRUE.
-      Agrif_SpecialValueFineGrid = -9999.
+      Agrif_SpecialValueFineGrid    = -9999.
 # if defined TWO_WAY
       IF( MOD(nbcline,nbclineupdate) == 0) THEN ! update the whole basin at each nbclineupdate (=nn_cln_update) baroclinic parent time steps
@@ -75,7 +73,4 @@
 
 
-   !!------------------
-   !! Local subroutines
-   !!------------------
    SUBROUTINE update_tra_ice( ptab, i1, i2, j1, j2, k1, k2, before )
       !!-----------------------------------------------------------------------
@@ -84,7 +79,7 @@
       !!              the properties per mass on the coarse grid
       !!-----------------------------------------------------------------------
-      INTEGER , INTENT(in) :: i1, i2, j1, j2, k1, k2
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
-      LOGICAL , INTENT(in) :: before
+      INTEGER                               , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2
+      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
+      LOGICAL                               , INTENT(in   ) ::   before
       !!
       INTEGER  :: jk, jl, jm
@@ -94,36 +89,36 @@
          jm = 1
          DO jl = 1, jpl
-            ptab(:,:,jm) = a_i  (i1:i2,j1:j2,jl) ; jm = jm + 1
-            ptab(:,:,jm) = v_i  (i1:i2,j1:j2,jl) ; jm = jm + 1
-            ptab(:,:,jm) = v_s  (i1:i2,j1:j2,jl) ; jm = jm + 1
-            ptab(:,:,jm) = smv_i(i1:i2,j1:j2,jl) ; jm = jm + 1
-            ptab(:,:,jm) = oa_i (i1:i2,j1:j2,jl) ; jm = jm + 1
+            ptab(:,:,jm) = a_i  (i1:i2,j1:j2,jl)   ;   jm = jm + 1
+            ptab(:,:,jm) = v_i  (i1:i2,j1:j2,jl)   ;   jm = jm + 1
+            ptab(:,:,jm) = v_s  (i1:i2,j1:j2,jl)   ;   jm = jm + 1
+            ptab(:,:,jm) = smv_i(i1:i2,j1:j2,jl)   ;   jm = jm + 1
+            ptab(:,:,jm) = oa_i (i1:i2,j1:j2,jl)   ;   jm = jm + 1
             DO jk = 1, nlay_s
-               ptab(:,:,jm) = e_s(i1:i2,j1:j2,jk,jl) ; jm = jm + 1
-            ENDDO
+               ptab(:,:,jm) = e_s(i1:i2,j1:j2,jk,jl)   ;   jm = jm + 1
+            END DO
             DO jk = 1, nlay_i
-               ptab(:,:,jm) = e_i(i1:i2,j1:j2,jk,jl) ; jm = jm + 1
-            ENDDO
-         ENDDO
+               ptab(:,:,jm) = e_i(i1:i2,j1:j2,jk,jl)   ;   jm = jm + 1
+            END DO
+         END DO
 
          DO jk = k1, k2
             WHERE( tmask(i1:i2,j1:j2,1) == 0. )  ptab(:,:,jk) = -9999.
-         ENDDO
-                  
+         END DO
+         !       
       ELSE
          jm = 1
          DO jl = 1, jpl
-            a_i  (i1:i2,j1:j2,jl) = ptab(:,:,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
-            v_i  (i1:i2,j1:j2,jl) = ptab(:,:,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
-            v_s  (i1:i2,j1:j2,jl) = ptab(:,:,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
-            smv_i(i1:i2,j1:j2,jl) = ptab(:,:,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
-            oa_i (i1:i2,j1:j2,jl) = ptab(:,:,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
+            a_i  (i1:i2,j1:j2,jl) = ptab(:,:,jm) * tmask(i1:i2,j1:j2,1)   ;   jm = jm + 1
+            v_i  (i1:i2,j1:j2,jl) = ptab(:,:,jm) * tmask(i1:i2,j1:j2,1)   ;   jm = jm + 1
+            v_s  (i1:i2,j1:j2,jl) = ptab(:,:,jm) * tmask(i1:i2,j1:j2,1)   ;   jm = jm + 1
+            smv_i(i1:i2,j1:j2,jl) = ptab(:,:,jm) * tmask(i1:i2,j1:j2,1)   ;   jm = jm + 1
+            oa_i (i1:i2,j1:j2,jl) = ptab(:,:,jm) * tmask(i1:i2,j1:j2,1)   ;   jm = jm + 1
             DO jk = 1, nlay_s
-               e_s(i1:i2,j1:j2,jk,jl) = ptab(:,:,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
+               e_s(i1:i2,j1:j2,jk,jl) = ptab(:,:,jm) * tmask(i1:i2,j1:j2,1)   ;   jm = jm + 1
             ENDDO
             DO jk = 1, nlay_i
-               e_i(i1:i2,j1:j2,jk,jl) = ptab(:,:,jm) * tmask(i1:i2,j1:j2,1) ; jm = jm + 1
-            ENDDO
-         ENDDO
+               e_i(i1:i2,j1:j2,jk,jl) = ptab(:,:,jm) * tmask(i1:i2,j1:j2,1)   ;   jm = jm + 1
+            END DO
+         END DO
 
          ! integrated values
@@ -144,9 +139,9 @@
       !! ** Method  : Update the fluxes and recover the properties (C-grid)
       !!-----------------------------------------------------------------------
-      INTEGER , INTENT(in) :: i1, i2, j1, j2
-      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
-      LOGICAL , INTENT(in) :: before
+      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
+      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
+      LOGICAL                         , INTENT(in   ) ::   before
       !!
-      REAL(wp) :: zrhoy
+      REAL(wp) ::   zrhoy   ! local scalar
       !!-----------------------------------------------------------------------
       !
@@ -154,5 +149,5 @@
          zrhoy = Agrif_Rhoy()
          ptab(:,:) = e2u(i1:i2,j1:j2) * u_ice(i1:i2,j1:j2) * zrhoy
-         WHERE( umask(i1:i2,j1:j2,1) == 0. )  ptab(:,:) = -9999.
+         WHERE( umask(i1:i2,j1:j2,1) == 0. )   ptab(:,:) = -9999.
       ELSE
          u_ice(i1:i2,j1:j2) = ptab(:,:) / e2u(i1:i2,j1:j2) * umask(i1:i2,j1:j2,1)
@@ -167,9 +162,9 @@
       !! ** Method  : Update the fluxes and recover the properties (C-grid)
       !!-----------------------------------------------------------------------
-      INTEGER , INTENT(in) :: i1,i2,j1,j2
-      REAL(wp), DIMENSION(i1:i2,j1:j2),  INTENT(inout) :: ptab
-      LOGICAL , INTENT(in) :: before
+      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
+      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
+      LOGICAL                         , INTENT(in   ) ::   before
       !!
-      REAL(wp) :: zrhox
+      REAL(wp) ::   zrhox   ! local scalar
       !!-----------------------------------------------------------------------
       !
@@ -177,5 +172,5 @@
          zrhox = Agrif_Rhox()
          ptab(:,:) = e1v(i1:i2,j1:j2) * v_ice(i1:i2,j1:j2) * zrhox
-         WHERE( vmask(i1:i2,j1:j2,1) == 0. )  ptab(:,:) = -9999.
+         WHERE( vmask(i1:i2,j1:j2,1) == 0. )   ptab(:,:) = -9999.
       ELSE
          v_ice(i1:i2,j1:j2) = ptab(:,:) / e1v(i1:i2,j1:j2) * vmask(i1:i2,j1:j2,1)
@@ -185,11 +180,13 @@
 
 #else
+   !!----------------------------------------------------------------------
+   !!   Empty module                                             no sea-ice
+   !!----------------------------------------------------------------------
 CONTAINS
    SUBROUTINE agrif_lim3_update_empty
-      !!---------------------------------------------
-      !!   *** ROUTINE agrif_lim3_update_empty ***
-      !!---------------------------------------------
       WRITE(*,*)  'agrif_lim3_update : You should not have seen this print! error?'
    END SUBROUTINE agrif_lim3_update_empty
 #endif
+
+   !!======================================================================
 END MODULE agrif_lim3_update
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/NST_SRC/agrif_oce.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/NST_SRC/agrif_oce.F90	(revision 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/NST_SRC/agrif_oce.F90	(revision 7953)
@@ -47,6 +47,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE,  DIMENSION(:,:) :: fsahm_spt, fsahm_spf !: sponge viscosities
 
-   ! Barotropic arrays used to store open boundary data during
-   ! time-splitting loop:
+   ! Barotropic arrays used to store open boundary data during time-splitting loop:
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) ::  ubdy_w, vbdy_w, hbdy_w
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) ::  ubdy_e, vbdy_e, hbdy_e
@@ -54,20 +53,22 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) ::  ubdy_s, vbdy_s, hbdy_s
 
-   INTEGER :: tsn_id                                                  ! AGRIF profile for tracers interpolation and update
-   INTEGER :: un_interp_id, vn_interp_id                              ! AGRIF profiles for interpolations
-   INTEGER :: un_update_id, vn_update_id                              ! AGRIF profiles for udpates
-   INTEGER :: tsn_sponge_id, un_sponge_id, vn_sponge_id               ! AGRIF profiles for sponge layers
+!!gm   add PUBLIC in all variable below:
+
+   INTEGER, PUBLIC :: tsn_id                                                  ! AGRIF profile for tracers interpolation and update
+   INTEGER, PUBLIC :: un_interp_id, vn_interp_id                              ! AGRIF profiles for interpolations
+   INTEGER, PUBLIC :: un_update_id, vn_update_id                              ! AGRIF profiles for udpates
+   INTEGER, PUBLIC :: tsn_sponge_id, un_sponge_id, vn_sponge_id               ! AGRIF profiles for sponge layers
 # if defined key_top
-   INTEGER :: trn_id, trn_sponge_id
+   INTEGER, PUBLIC :: trn_id, trn_sponge_id
 # endif  
-   INTEGER :: unb_id, vnb_id, ub2b_interp_id, vb2b_interp_id
-   INTEGER :: ub2b_update_id, vb2b_update_id
-   INTEGER :: e3t_id, e1u_id, e2v_id, sshn_id
-   INTEGER :: scales_t_id
-# if defined key_zdftke
-   INTEGER :: avt_id, avm_id, en_id
-# endif  
-   INTEGER :: umsk_id, vmsk_id
-   INTEGER :: kindic_agr
+   INTEGER, PUBLIC :: unb_id, vnb_id, ub2b_interp_id, vb2b_interp_id
+   INTEGER, PUBLIC :: ub2b_update_id, vb2b_update_id
+   INTEGER, PUBLIC :: e3t_id, e1u_id, e2v_id, sshn_id
+   INTEGER, PUBLIC :: scales_t_id
+   INTEGER, PUBLIC :: avt_id, avm_id, en_id                ! TKE related identificators
+   INTEGER, PUBLIC :: umsk_id, vmsk_id
+   INTEGER, PUBLIC :: kindic_agr
+   
+!!gm end public addition
 
    !!----------------------------------------------------------------------
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90	(revision 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90	(revision 7953)
@@ -21,5 +21,5 @@
    USE oce
    USE dom_oce      
-   USE zdf_oce
+   USE zdf_oce          ! vertical physics
    USE agrif_oce
    USE phycst
@@ -34,11 +34,9 @@
 
    PUBLIC   Agrif_tra, Agrif_dyn, Agrif_ssh, Agrif_dyn_ts, Agrif_ssh_ts, Agrif_dta_ts
-   PUBLIC   interpun, interpvn
-   PUBLIC   interptsn,  interpsshn
-   PUBLIC   interpunb, interpvnb, interpub2b, interpvb2b
+   PUBLIC   interpun , interpvn
+   PUBLIC   interptsn, interpsshn
+   PUBLIC   interpunb, interpvnb , interpub2b, interpvb2b
    PUBLIC   interpe3t, interpumsk, interpvmsk
-# if defined key_zdftke
    PUBLIC   Agrif_tke, interpavm
-# endif
 
    INTEGER ::   bdy_tinterp = 0
@@ -46,5 +44,5 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/NST 3.7 , NEMO Consortium (2015)
+   !! NEMO/NST 4.0 , NEMO Consortium (2017)
    !! $Id$
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
@@ -385,5 +383,4 @@
       !!                  ***  ROUTINE Agrif_dyn_ts  ***
       !!----------------------------------------------------------------------  
-      !! 
       INTEGER, INTENT(in) ::   jn
       !!
@@ -444,5 +441,4 @@
       !!                  ***  ROUTINE Agrif_dta_ts  ***
       !!----------------------------------------------------------------------  
-      !! 
       INTEGER, INTENT(in) ::   kt
       !!
@@ -504,5 +500,4 @@
       !!----------------------------------------------------------------------  
       INTEGER, INTENT(in) ::   kt
-      !!
       !!----------------------------------------------------------------------  
       !
@@ -541,5 +536,5 @@
       !!----------------------------------------------------------------------  
       !
-      IF((nbondi == -1).OR.(nbondi == 2)) THEN
+      IF( nbondi == -1 .OR. nbondi == 2 ) THEN
          DO jj = 1, jpj
             ssha_e(2,jj) = hbdy_w(jj)
@@ -547,5 +542,5 @@
       ENDIF
       !
-      IF((nbondi == 1).OR.(nbondi == 2)) THEN
+      IF( nbondi == 1 .OR. nbondi == 2 ) THEN
          DO jj = 1, jpj
             ssha_e(nlci-1,jj) = hbdy_e(jj)
@@ -553,5 +548,5 @@
       ENDIF
       !
-      IF((nbondj == -1).OR.(nbondj == 2)) THEN
+      IF( nbondj == -1 .OR.(nbondj == 2 ) THEN
          DO ji = 1, jpi
             ssha_e(ji,2) = hbdy_s(ji)
@@ -559,5 +554,5 @@
       ENDIF
       !
-      IF((nbondj == 1).OR.(nbondj == 2)) THEN
+      IF( nbondj == 1 .OR. nbondj == 2 ) THEN
          DO ji = 1, jpi
             ssha_e(ji,nlcj-1) = hbdy_n(ji)
@@ -567,5 +562,4 @@
    END SUBROUTINE Agrif_ssh_ts
 
-# if defined key_zdftke
 
    SUBROUTINE Agrif_tke
@@ -579,18 +573,17 @@
       IF( zalpha > 1. )   zalpha = 1.
       !
-      Agrif_SpecialValue    = 0.e0
+      Agrif_SpecialValue    = 0._wp
       Agrif_UseSpecialValue = .TRUE.
       !
-      CALL Agrif_Bc_variable(avm_id ,calledweight=zalpha, procname=interpavm)       
+      CALL Agrif_Bc_variable( avm_id , calledweight=zalpha, procname=interpavm )       
       !
       Agrif_UseSpecialValue = .FALSE.
       !
    END SUBROUTINE Agrif_tke
-   
-# endif
+
 
    SUBROUTINE interptsn( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir )
       !!----------------------------------------------------------------------
-      !!   *** ROUTINE interptsn ***
+      !!                  *** ROUTINE interptsn ***
       !!----------------------------------------------------------------------
       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) ::   ptab
@@ -599,9 +592,9 @@
       INTEGER                                     , INTENT(in   ) ::   nb , ndir
       !
-      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
-      INTEGER  ::   imin, imax, jmin, jmax
-      REAL(wp) ::   zrhox , zalpha1, zalpha2, zalpha3
-      REAL(wp) ::   zalpha4, zalpha5, zalpha6, zalpha7
-      LOGICAL  ::   western_side, eastern_side,northern_side,southern_side
+      INTEGER ::   ji, jj, jk, jn   ! dummy loop indices
+      INTEGER ::   imin, imax, jmin, jmax
+      REAL(wp)::   zrhox , zalpha1, zalpha2, zalpha3
+      REAL(wp)::   zalpha4, zalpha5, zalpha6, zalpha7
+      LOGICAL ::   western_side, eastern_side,northern_side,southern_side
       !!----------------------------------------------------------------------
       !
@@ -770,5 +763,5 @@
    SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, before )
       !!----------------------------------------------------------------------
-      !!   *** ROUTINE interpun ***
+      !!                         *** ROUTINE interpun ***
       !!----------------------------------------------------------------------
       INTEGER                               , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2
@@ -776,6 +769,6 @@
       LOGICAL                               , INTENT(in   ) ::   before
       !
-      INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   zrhoy	
+      INTEGER ::   ji, jj, jk
+      REAL(wp)::   zrhoy	
       !!----------------------------------------------------------------------
       !
@@ -798,5 +791,5 @@
    SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, before )
       !!----------------------------------------------------------------------
-      !!   *** ROUTINE interpvn ***
+      !!                      *** ROUTINE interpvn ***
       !!----------------------------------------------------------------------
       INTEGER                               , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2
@@ -804,6 +797,6 @@
       LOGICAL                               , INTENT(in   ) ::   before
       !
-      INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   zrhox	
+      INTEGER ::   ji, jj, jk
+      REAL(wp)::   zrhox	
       !!----------------------------------------------------------------------
       !      
@@ -831,7 +824,7 @@
       INTEGER                         , INTENT(in   ) ::   nb , ndir
       !
-      INTEGER  ::   ji, jj
-      REAL(wp) ::   zrhoy, zrhot, zt0, zt1, ztcoeff
-      LOGICAL  ::   western_side, eastern_side,northern_side,southern_side
+      INTEGER ::   ji, jj
+      REAL(wp)::   zrhoy, zrhot, zt0, zt1, ztcoeff
+      LOGICAL ::   western_side, eastern_side,northern_side,southern_side
       !!----------------------------------------------------------------------  
       !
@@ -901,7 +894,7 @@
       INTEGER                         , INTENT(in   ) ::   nb , ndir
       !
-      INTEGER  ::   ji,jj
-      REAL(wp) ::   zrhox, zrhot, zt0, zt1, ztcoeff   
-      LOGICAL  ::   western_side, eastern_side,northern_side,southern_side
+      INTEGER ::   ji,jj
+      REAL(wp)::   zrhox, zrhot, zt0, zt1, ztcoeff   
+      LOGICAL ::   western_side, eastern_side,northern_side,southern_side
       !!----------------------------------------------------------------------  
       ! 
@@ -919,9 +912,9 @@
          zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot      
          IF( bdy_tinterp == 1 ) THEN
-            ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
-               &               - zt0**2._wp * (       zt0 - 1._wp)        )
+            ztcoeff = zrhot * (  zt1**2._wp * ( zt1 - 1._wp)        &
+               &               - zt0**2._wp * ( zt0 - 1._wp)        )
          ELSEIF( bdy_tinterp == 2 ) THEN
-            ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
-               &               - zt0        * (       zt0 - 1._wp)**2._wp ) 
+            ztcoeff = zrhot * (  zt1        * ( zt1 - 1._wp)**2._wp &
+               &               - zt0        * ( zt0 - 1._wp)**2._wp ) 
          ELSE
             ztcoeff = 1
@@ -942,19 +935,8 @@
          !            
          IF( bdy_tinterp == 0 .OR. bdy_tinterp == 2) THEN
-            IF(western_side) THEN
-               vbdy_w(j1:j2) = vbdy_w(j1:j2) / (zrhox*e1v(i1,j1:j2))   &
-                     &                                  * vmask(i1,j1:j2,1)
-            ENDIF
-            IF(eastern_side) THEN
-               vbdy_e(j1:j2) = vbdy_e(j1:j2) / (zrhox*e1v(i1,j1:j2))   &
-                     &                                  * vmask(i1,j1:j2,1)
-            ENDIF
-            IF(southern_side) THEN
-               vbdy_s(i1:i2) = vbdy_s(i1:i2) / (zrhox*e1v(i1:i2,j1))   &
-                     &                                  * vmask(i1:i2,j1,1)
-            ENDIF
-            IF(northern_side) THEN
-               vbdy_n(i1:i2) = vbdy_n(i1:i2) / (zrhox*e1v(i1:i2,j1))   &
-                     &                                  * vmask(i1:i2,j1,1)
+            IF( western_side  )   vbdy_w(j1:j2) = vbdy_w(j1:j2) / (zrhox*e1v(i1,j1:j2)) * vmask(i1,j1:j2,1)
+            IF( eastern_side  )   vbdy_e(j1:j2) = vbdy_e(j1:j2) / (zrhox*e1v(i1,j1:j2)) * vmask(i1,j1:j2,1)
+            IF( southern_side )   vbdy_s(i1:i2) = vbdy_s(i1:i2) / (zrhox*e1v(i1:i2,j1)) * vmask(i1:i2,j1,1)
+            IF( northern_side )   vbdy_n(i1:i2) = vbdy_n(i1:i2) / (zrhox*e1v(i1:i2,j1)) * vmask(i1:i2,j1,1)
             ENDIF
          ENDIF
@@ -973,7 +955,7 @@
       INTEGER                         , INTENT(in   ) ::   nb , ndir
       !
-      INTEGER  ::   ji,jj
-      REAL(wp) ::   zrhot, zt0, zt1,zat
-      LOGICAL  ::   western_side, eastern_side,northern_side,southern_side
+      INTEGER ::   ji,jj
+      REAL(wp)::   zrhot, zt0, zt1,zat
+      LOGICAL ::   western_side, eastern_side,northern_side,southern_side
       !!----------------------------------------------------------------------  
       IF( before ) THEN
@@ -1030,8 +1012,8 @@
             &           - zt0**2._wp * (-2._wp*zt0 + 3._wp)    ) 
          !
-         IF(western_side )   vbdy_w(j1:j2) = zat * ptab(i1,j1:j2)  
-         IF(eastern_side )   vbdy_e(j1:j2) = zat * ptab(i1,j1:j2)  
-         IF(southern_side)   vbdy_s(i1:i2) = zat * ptab(i1:i2,j1) 
-         IF(northern_side)   vbdy_n(i1:i2) = zat * ptab(i1:i2,j1) 
+         IF( western_side  )   vbdy_w(j1:j2) = zat * ptab(i1,j1:j2)  
+         IF( eastern_side  )   vbdy_e(j1:j2) = zat * ptab(i1,j1:j2)  
+         IF( southern_side )   vbdy_s(i1:i2) = zat * ptab(i1:i2,j1) 
+         IF( northern_side )   vbdy_n(i1:i2) = zat * ptab(i1:i2,j1) 
       ENDIF
       !      
@@ -1048,7 +1030,7 @@
       INTEGER                              , INTENT(in   ) :: nb , ndir
       !
-      INTEGER :: ji, jj, jk
-      LOGICAL :: western_side, eastern_side, northern_side, southern_side
-      REAL(wp) :: ztmpmsk      
+      INTEGER ::   ji, jj, jk
+      LOGICAL ::   western_side, eastern_side, northern_side, southern_side
+      REAL(wp)::   ztmpmsk      
       !!----------------------------------------------------------------------  
       !    
@@ -1065,8 +1047,8 @@
                DO ji = i1, i2
                   ! Get velocity mask at boundary edge points:
-                  IF( western_side )   ztmpmsk = umask(ji    ,jj    ,1)
-                  IF( eastern_side )   ztmpmsk = umask(nlci-2,jj    ,1)
-                  IF( northern_side)   ztmpmsk = vmask(ji    ,nlcj-2,1)
-                  IF( southern_side)   ztmpmsk = vmask(ji    ,2     ,1)
+                  IF( western_side  )   ztmpmsk = umask(ji    ,jj    ,1)
+                  IF( eastern_side  )   ztmpmsk = umask(nlci-2,jj    ,1)
+                  IF( northern_side )   ztmpmsk = vmask(ji    ,nlcj-2,1)
+                  IF( southern_side )   ztmpmsk = vmask(ji    ,2     ,1)
                   !
                   IF( ABS( ptab(ji,jj,jk) - tmask(ji,jj,jk) * e3t_0(ji,jj,jk) )*ztmpmsk > 1.D-2) THEN
@@ -1141,5 +1123,5 @@
       REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
       LOGICAL                              , INTENT(in   ) ::   before
-      INTEGER                              , INTENT(in   ) :: nb , ndir
+      INTEGER                              , INTENT(in   ) ::   nb , ndir
       !
       INTEGER ::   ji, jj, jk
@@ -1175,5 +1157,4 @@
    END SUBROUTINE interpvmsk
 
-# if defined key_zdftke
 
    SUBROUTINE interpavm( ptab, i1, i2, j1, j2, k1, k2, before )
@@ -1186,13 +1167,9 @@
       !!----------------------------------------------------------------------  
       !      
-      IF( before ) THEN
-         ptab (i1:i2,j1:j2,k1:k2) = avm_k(i1:i2,j1:j2,k1:k2)
-      ELSE
-         avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2)
+      IF( before ) THEN   ;   ptab (i1:i2,j1:j2,k1:k2) = avm_k(i1:i2,j1:j2,k1:k2)
+      ELSE                ;   avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2)
       ENDIF
       !
    END SUBROUTINE interpavm
-
-# endif /* key_zdftke */
 
 #else
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90	(revision 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90	(revision 7953)
@@ -3,13 +3,20 @@
 MODULE agrif_opa_sponge
    !!======================================================================
-   !!                ***  MODULE agrif_opa_update  ***
-   !! AGRIF :   
+   !!                   ***  MODULE  agrif_opa_interp  ***
+   !! AGRIF: interpolation package
    !!======================================================================
-   !! History :  
+   !! History :  2.0  !  2002-06  (XXX)  Original cade
+   !!             -   !  2005-11  (XXX) 
+   !!            3.2  !  2009-04  (R. Benshila) 
+   !!            3.6  !  2014-09  (R. Benshila) 
    !!----------------------------------------------------------------------
 #if defined key_agrif
+   !!----------------------------------------------------------------------
+   !!   'key_agrif'                                              AGRIF zoom
+   !!----------------------------------------------------------------------
    USE par_oce
    USE oce
    USE dom_oce
+   !
    USE in_out_manager
    USE agrif_oce
@@ -24,5 +31,5 @@
 
    !!----------------------------------------------------------------------
-   !! NEMO/NST 3.7 , NEMO Consortium (2015)
+   !! NEMO/NST 4.0 , NEMO Consortium (2017)
    !! $Id$
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
@@ -31,20 +38,20 @@
 
    SUBROUTINE Agrif_Sponge_Tra
-      !!---------------------------------------------
-      !!   *** ROUTINE Agrif_Sponge_Tra ***
-      !!---------------------------------------------
-      REAL(wp) :: timecoeff
-      !!---------------------------------------------
+      !!----------------------------------------------------------------------
+      !!                 *** ROUTINE Agrif_Sponge_Tra ***
+      !!----------------------------------------------------------------------
+      REAL(wp) ::   timecoeff   ! local scalar
+      !!----------------------------------------------------------------------
       !
 #if defined SPONGE
       timecoeff = REAL(Agrif_NbStepint(),wp)/Agrif_rhot()
-
+      !
       CALL Agrif_Sponge
-      Agrif_SpecialValue=0.
+      Agrif_SpecialValue    = 0._wp
       Agrif_UseSpecialValue = .TRUE.
-      tabspongedone_tsn = .FALSE.
-
+      tabspongedone_tsn     = .FALSE.
+      !
       CALL Agrif_Bc_Variable(tsn_sponge_id,calledweight=timecoeff,procname=interptsn_sponge)
-
+      !
       Agrif_UseSpecialValue = .FALSE.
 #endif
@@ -54,24 +61,24 @@
 
    SUBROUTINE Agrif_Sponge_dyn
-      !!---------------------------------------------
-      !!   *** ROUTINE Agrif_Sponge_dyn ***
-      !!---------------------------------------------
-      REAL(wp) :: timecoeff
-      !!---------------------------------------------
-
+      !!----------------------------------------------------------------------
+      !!                 *** ROUTINE Agrif_Sponge_dyn ***
+      !!----------------------------------------------------------------------
+      REAL(wp) ::   timecoeff   ! local scalar
+      !!----------------------------------------------------------------------
+      !
 #if defined SPONGE
       timecoeff = REAL(Agrif_NbStepint(),wp)/Agrif_rhot()
-
-      Agrif_SpecialValue=0.
+      !
+      Agrif_SpecialValue    = 0._wp
       Agrif_UseSpecialValue = ln_spc_dyn
-
+      !
       tabspongedone_u = .FALSE.
       tabspongedone_v = .FALSE.         
       CALL Agrif_Bc_Variable(un_sponge_id,calledweight=timecoeff,procname=interpun_sponge)
-
+      !
       tabspongedone_u = .FALSE.
       tabspongedone_v = .FALSE.
       CALL Agrif_Bc_Variable(vn_sponge_id,calledweight=timecoeff,procname=interpvn_sponge)
-
+      !
       Agrif_UseSpecialValue = .FALSE.
 #endif
@@ -81,7 +88,7 @@
 
    SUBROUTINE Agrif_Sponge
-      !!---------------------------------------------
-      !!   *** ROUTINE  Agrif_Sponge ***
-      !!---------------------------------------------
+      !!----------------------------------------------------------------------
+      !!                 *** ROUTINE  Agrif_Sponge ***
+      !!----------------------------------------------------------------------
       INTEGER  :: ji,jj,jk
       INTEGER  :: ispongearea, ilci, ilcj
@@ -89,5 +96,6 @@
       REAL(wp) :: z1spongearea, zramp
       REAL(wp), POINTER, DIMENSION(:,:) :: ztabramp
-
+      !!----------------------------------------------------------------------
+      !
 #if defined SPONGE || defined SPONGE_TOP
       ll_spdone=.TRUE.
@@ -176,8 +184,8 @@
                fsahm_spt(ji,jj) = visc_dyn * ztabramp(ji,jj)
                fsahm_spf(ji,jj) = 0.25_wp * visc_dyn * ( ztabramp(ji,jj) + ztabramp(ji  ,jj+1) &
-                                                     &  +ztabramp(ji,jj) + ztabramp(ji+1,jj  ) )
-            END DO
-         END DO
-
+                  &                                     +ztabramp(ji,jj) + ztabramp(ji+1,jj  ) )
+            END DO
+         END DO
+         !
          CALL lbc_lnk( fsahm_spt, 'T', 1. )   ! Lateral boundary conditions
          CALL lbc_lnk( fsahm_spf, 'F', 1. )
@@ -192,11 +200,11 @@
 
 
-   SUBROUTINE interptsn_sponge(tabres,i1,i2,j1,j2,k1,k2,n1,n2,before)
-      !!---------------------------------------------
-      !!   *** ROUTINE interptsn_sponge ***
-      !!---------------------------------------------
-      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
-      LOGICAL, INTENT(in) :: before
+   SUBROUTINE interptsn_sponge( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
+      !!----------------------------------------------------------------------
+      !!                 *** ROUTINE interptsn_sponge ***
+      !!----------------------------------------------------------------------
+      INTEGER                                     , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2
+      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) ::   tabres
+      LOGICAL                                     , INTENT(in   ) ::   before
       !
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
@@ -205,4 +213,5 @@
       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ztu, ztv
       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2) ::tsbdiff
+      !!----------------------------------------------------------------------
       !
       IF( before ) THEN
@@ -258,20 +267,18 @@
 
 
-   SUBROUTINE interpun_sponge(tabres,i1,i2,j1,j2,k1,k2, before)
-      !!---------------------------------------------
-      !!   *** ROUTINE interpun_sponge ***
-      !!---------------------------------------------    
-      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
-      LOGICAL, INTENT(in) :: before
-
+   SUBROUTINE interpun_sponge( tabres, i1, i2, j1, j2, k1, k2, before )
+      !!----------------------------------------------------------------------
+      !!                 *** ROUTINE interpun_sponge ***
+      !!----------------------------------------------------------------------
+      INTEGER                               , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2
+      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   tabres
+      LOGICAL                               , INTENT(in   ) ::   before
+      !!
       INTEGER :: ji,jj,jk
-
-      ! sponge parameters 
+      INTEGER :: jmax
       REAL(wp) :: ze2u, ze1v, zua, zva, zbtr
       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ubdiff
       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: rotdiff, hdivdiff
-      INTEGER :: jmax
-      !!---------------------------------------------    
+      !!----------------------------------------------------------------------
       !
       IF( before ) THEN
@@ -356,19 +363,18 @@
 
 
-   SUBROUTINE interpvn_sponge(tabres,i1,i2,j1,j2,k1,k2, before,nb,ndir)
-      !!---------------------------------------------
-      !!   *** ROUTINE interpvn_sponge ***
-      !!--------------------------------------------- 
-      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
-      LOGICAL, INTENT(in) :: before
-      INTEGER, INTENT(in) :: nb , ndir
-      !
-      INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   ze2u, ze1v, zua, zva, zbtr
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: vbdiff
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: rotdiff, hdivdiff
-      INTEGER :: imax
-      !!--------------------------------------------- 
+   SUBROUTINE interpvn_sponge( tabres, i1, i2, j1, j2, k1, k2, before, nb, ndir )
+      !!----------------------------------------------------------------------
+      !!                 *** ROUTINE interpvn_sponge ***
+      !!----------------------------------------------------------------------
+      INTEGER                               , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2
+      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   tabres
+      LOGICAL                               , INTENT(in   ) ::   before
+      INTEGER                               , INTENT(in   ) ::   nb , ndir
+      !
+      INTEGER ::   ji, jj, jk
+      INTEGER ::   imax
+      REAL(wp)::   ze2u, ze1v, zua, zva, zbtr
+      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) ::   vbdiff, rotdiff, hdivdiff
+      !!----------------------------------------------------------------------
 
       IF( before ) THEN 
@@ -403,5 +409,5 @@
          !                                                
 
-         imax = i2-1
+         imax = i2 - 1
          IF ((nbondi == 1).OR.(nbondi == 2))   imax = MIN(imax,nlci-3)
 
@@ -437,9 +443,12 @@
 
 #else
+   !!----------------------------------------------------------------------
+   !!   Empty module                                          no AGRIF zoom
+   !!----------------------------------------------------------------------
 CONTAINS
    SUBROUTINE agrif_opa_sponge_empty
-      !!---------------------------------------------
-      !!   *** ROUTINE agrif_OPA_sponge_empty ***
-      !!---------------------------------------------
+      !!----------------------------------------------------------------------
+      !!                 *** ROUTINE agrif_OPA_sponge_empty ***
+      !!----------------------------------------------------------------------
       WRITE(*,*)  'agrif_opa_sponge : You should not have seen this print! error?'
    END SUBROUTINE agrif_opa_sponge_empty
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90	(revision 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90	(revision 7953)
@@ -3,23 +3,35 @@
  
 MODULE agrif_opa_update
+   !!======================================================================
+   !!                   ***  MODULE  agrif_opa_interp  ***
+   !! AGRIF: interpolation package
+   !!======================================================================
+   !! History :  2.0  !  2002-06  (XXX)  Original cade
+   !!             -   !  2005-11  (XXX) 
+   !!            3.2  !  2009-04  (R. Benshila) 
+   !!            3.6  !  2014-09  (R. Benshila) 
+   !!----------------------------------------------------------------------
 #if defined key_agrif 
+   !!----------------------------------------------------------------------
+   !!   'key_agrif'                                              AGRIF zoom
+   !!----------------------------------------------------------------------
    USE par_oce
    USE oce
    USE dom_oce
+   USE zdf_oce        ! vertical physics: ocean variables 
    USE agrif_oce
-   USE in_out_manager  ! I/O manager
+   !
+   USE in_out_manager ! I/O manager
    USE lib_mpp
    USE wrk_nemo  
-   USE zdf_oce        ! vertical physics: ocean variables 
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC Agrif_Update_Tra, Agrif_Update_Dyn,Update_Scales
-# if defined key_zdftke
-   PUBLIC Agrif_Update_Tke
-# endif
+   PUBLIC   Agrif_Update_Tra, Agrif_Update_Dyn, Update_Scales
+   PUBLIC   Agrif_Update_Tke
+
    !!----------------------------------------------------------------------
-   !! NEMO/NST 3.6 , NEMO Consortium (2010)
+   !! NEMO/NST 4.0 , NEMO Consortium (2017)
    !! $Id$
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
@@ -28,7 +40,7 @@
 
    RECURSIVE SUBROUTINE Agrif_Update_Tra( )
-      !!---------------------------------------------
-      !!   *** ROUTINE Agrif_Update_Tra ***
-      !!---------------------------------------------
+      !!----------------------------------------------------------------------
+      !!                   *** ROUTINE Agrif_Update_Tra ***
+      !!----------------------------------------------------------------------
       ! 
       IF (Agrif_Root()) RETURN
@@ -38,5 +50,5 @@
 
       Agrif_UseSpecialValueInUpdate = .TRUE.
-      Agrif_SpecialValueFineGrid = 0.
+      Agrif_SpecialValueFineGrid    = 0._wp
       ! 
       IF (MOD(nbcline,nbclineupdate) == 0) THEN
@@ -68,7 +80,7 @@
 
    RECURSIVE SUBROUTINE Agrif_Update_Dyn( )
-      !!---------------------------------------------
-      !!   *** ROUTINE Agrif_Update_Dyn ***
-      !!---------------------------------------------
+      !!----------------------------------------------------------------------
+      !!                   *** ROUTINE Agrif_Update_Dyn ***
+      !!----------------------------------------------------------------------
       ! 
       IF (Agrif_Root()) RETURN
@@ -106,5 +118,5 @@
 # endif
 
-      IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN
+      IF ( ln_dynspg_ts .AND. ln_bt_fw ) THEN
          ! Update time integrated transports
          IF (mod(nbcline,nbclineupdate) == 0) THEN
@@ -149,47 +161,47 @@
    END SUBROUTINE Agrif_Update_Dyn
 
-# if defined key_zdftke
+!!gm Missing GLS case !!!!!
 
    SUBROUTINE Agrif_Update_Tke( kt )
-      !!---------------------------------------------
-      !!   *** ROUTINE Agrif_Update_Tke ***
-      !!---------------------------------------------
-      !!
+      !!----------------------------------------------------------------------
+      !!                   *** ROUTINE Agrif_Update_Tke ***
+      !!----------------------------------------------------------------------
       INTEGER, INTENT(in) :: kt
-      !       
-      IF( ( Agrif_NbStepint() .NE. 0 ) .AND. (kt /= 0) ) RETURN
+      !!----------------------------------------------------------------------
+      !
+!!gm test on kt/=0  ????  why not nit000-1  ?  doesn't seem logic
+      IF( ( Agrif_NbStepint() /= 0 ) .AND. kt /= 0 )   RETURN
 #  if defined TWO_WAY
-
+      !
       Agrif_UseSpecialValueInUpdate = .TRUE.
-      Agrif_SpecialValueFineGrid = 0.
-
-      CALL Agrif_Update_Variable( en_id, locupdate=(/0,0/), procname=updateEN  )
-      CALL Agrif_Update_Variable(avt_id, locupdate=(/0,0/), procname=updateAVT )
-      CALL Agrif_Update_Variable(avm_id, locupdate=(/0,0/), procname=updateAVM )
-
+      Agrif_SpecialValueFineGrid    = 0._wp
+      !
+      CALL Agrif_Update_Variable(  en_id, locupdate=(/0,0/), procname=updateEN  )
+      CALL Agrif_Update_Variable( avt_id, locupdate=(/0,0/), procname=updateAVT )
+      CALL Agrif_Update_Variable( avm_id, locupdate=(/0,0/), procname=updateAVM )
+      !
       Agrif_UseSpecialValueInUpdate = .FALSE.
-
+      !
 #  endif
-      
+      !
    END SUBROUTINE Agrif_Update_Tke
    
-# endif /* key_zdftke */
 
    SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
-      !!---------------------------------------------
+      !!----------------------------------------------------------------------
       !!           *** ROUTINE updateT ***
-      !!---------------------------------------------
-      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
-      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
-      LOGICAL, INTENT(in) :: before
-      !!
-      INTEGER :: ji,jj,jk,jn
-      !!---------------------------------------------
-      !
-      IF (before) THEN
-         DO jn = n1,n2
-            DO jk=k1,k2
-               DO jj=j1,j2
-                  DO ji=i1,i2
+      !!----------------------------------------------------------------------
+      INTEGER                                    , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2
+      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) ::   tabres
+      LOGICAL                                    , INTENT(in   ) ::   before
+      !
+      INTEGER :: ji, jj, jk, jn
+      !!----------------------------------------------------------------------
+      !
+      IF( before ) THEN
+         DO jn = n1, n2
+            DO jk = k1, k2
+               DO jj = j1, j2
+                  DO ji = i1, i2
                      tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)
                   END DO
@@ -209,8 +221,8 @@
                                  &             - tsn(ji,jj,jk,jn) ) * tmask(ji,jj,jk)
                         ENDIF
-                     ENDDO
-                  ENDDO
-               ENDDO
-            ENDDO
+                     END DO
+                  END DO
+               END DO
+            END DO
          ENDIF
          DO jn = n1,n2
@@ -238,6 +250,6 @@
       LOGICAL                               , INTENT(in   ) :: before
       !
-      INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   zrhoy
+      INTEGER ::   ji, jj, jk
+      REAL(wp)::   zrhoy
       !!---------------------------------------------
       ! 
@@ -268,16 +280,16 @@
 
    SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, before )
-      !!---------------------------------------------
-      !!           *** ROUTINE updatev ***
-      !!---------------------------------------------
-      INTEGER :: i1,i2,j1,j2,k1,k2
-      INTEGER :: ji,jj,jk
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: tabres
-      LOGICAL :: before
+      !!----------------------------------------------------------------------
+      !!                      *** ROUTINE updatev ***
+      !!----------------------------------------------------------------------
+      INTEGER                               , INTENT(in   ) :: i1, i2, j1, j2, k1, k2
+      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
+      LOGICAL                               , INTENT(in   ) :: before
       !!
-      REAL(wp) :: zrhox
-      !!---------------------------------------------      
-      !
-      IF (before) THEN
+      INTEGER  ::   ji, jj, jk
+      REAL(wp) ::   zrhox
+      !!----------------------------------------------------------------------
+      !
+      IF( before ) THEN
          zrhox = Agrif_Rhox()
          DO jk=k1,k2
@@ -309,17 +321,16 @@
 
    SUBROUTINE updateu2d( tabres, i1, i2, j1, j2, before )
+      !!----------------------------------------------------------------------
+      !!                      *** ROUTINE updateu2d ***
+      !!----------------------------------------------------------------------
+      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
+      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   tabres
+      LOGICAL                         , INTENT(in   ) ::   before
+      !! 
+      INTEGER ::   ji, jj, jk
+      REAL(wp)::   zrhoy, zcorr
       !!---------------------------------------------
-      !!          *** ROUTINE updateu2d ***
-      !!---------------------------------------------
-      INTEGER, INTENT(in) :: i1, i2, j1, j2
-      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
-      LOGICAL, INTENT(in) :: before
-      !! 
-      INTEGER :: ji, jj, jk
-      REAL(wp) :: zrhoy
-      REAL(wp) :: zcorr
-      !!---------------------------------------------
-      !
-      IF (before) THEN
+      !
+      IF( before ) THEN
          zrhoy = Agrif_Rhoy()
          DO jj=j1,j2
@@ -374,17 +385,16 @@
 
    SUBROUTINE updatev2d( tabres, i1, i2, j1, j2, before )
-      !!---------------------------------------------
-      !!          *** ROUTINE updatev2d ***
-      !!---------------------------------------------
-      INTEGER, INTENT(in) :: i1, i2, j1, j2
-      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
-      LOGICAL, INTENT(in) :: before
-      !! 
+      !!----------------------------------------------------------------------
+      !!                   *** ROUTINE updatev2d ***
+      !!----------------------------------------------------------------------
+      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
+      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   tabres
+      LOGICAL                         , INTENT(in   ) ::   before
+      ! 
       INTEGER :: ji, jj, jk
-      REAL(wp) :: zrhox
-      REAL(wp) :: zcorr
-      !!---------------------------------------------
-      !
-      IF (before) THEN
+      REAL(wp) :: zrhox, zcorr
+      !!----------------------------------------------------------------------
+      !
+      IF( before ) THEN
          zrhox = Agrif_Rhox()
          DO jj=j1,j2
@@ -439,15 +449,15 @@
 
    SUBROUTINE updateSSH( tabres, i1, i2, j1, j2, before )
-      !!---------------------------------------------
-      !!          *** ROUTINE updateSSH ***
-      !!---------------------------------------------
-      INTEGER, INTENT(in) :: i1, i2, j1, j2
-      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
-      LOGICAL, INTENT(in) :: before
+      !!----------------------------------------------------------------------
+      !!                   *** ROUTINE updateSSH ***
+      !!----------------------------------------------------------------------
+      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
+      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   tabres
+      LOGICAL                         , INTENT(in   ) ::   before
       !!
       INTEGER :: ji, jj
-      !!---------------------------------------------
-      ! 
-      IF (before) THEN
+      !!----------------------------------------------------------------------
+      ! 
+      IF( before ) THEN
          DO jj=j1,j2
             DO ji=i1,i2
@@ -478,14 +488,14 @@
 
    SUBROUTINE updateub2b( tabres, i1, i2, j1, j2, before )
-      !!---------------------------------------------
-      !!          *** ROUTINE updateub2b ***
-      !!---------------------------------------------
-      INTEGER, INTENT(in) :: i1, i2, j1, j2
-      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
-      LOGICAL, INTENT(in) :: before
+      !!----------------------------------------------------------------------
+      !!                      *** ROUTINE updateub2b ***
+      !!----------------------------------------------------------------------
+      INTEGER                            , INTENT(in) ::   i1, i2, j1, j2
+      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   tabres
+      LOGICAL                            , INTENT(in) ::   before
       !!
-      INTEGER :: ji, jj
-      REAL(wp) :: zrhoy
-      !!---------------------------------------------
+      INTEGER ::   ji, jj
+      REAL(wp)::   zrhoy
+      !!----------------------------------------------------------------------
       !
       IF (before) THEN
@@ -509,16 +519,16 @@
 
    SUBROUTINE updatevb2b( tabres, i1, i2, j1, j2, before )
-      !!---------------------------------------------
-      !!          *** ROUTINE updatevb2b ***
-      !!---------------------------------------------
-      INTEGER, INTENT(in) :: i1, i2, j1, j2
-      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
-      LOGICAL, INTENT(in) :: before
+      !!----------------------------------------------------------------------
+      !!                      *** ROUTINE updatevb2b ***
+      !!----------------------------------------------------------------------
+      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
+      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   tabres
+      LOGICAL                         , INTENT(in   ) ::   before
       !!
-      INTEGER :: ji, jj
-      REAL(wp) :: zrhox
-      !!---------------------------------------------
-      !
-      IF (before) THEN
+      INTEGER ::   ji, jj
+      REAL(wp)::   zrhox
+      !!----------------------------------------------------------------------
+      !
+      IF( before ) THEN
          zrhox = Agrif_Rhox()
          DO jj=j1,j2
@@ -540,15 +550,17 @@
 
    SUBROUTINE update_scales( tabres, i1, i2, j1, j2, k1, k2, n1,n2, before )
-      ! currently not used
-      !!---------------------------------------------
-      !!           *** ROUTINE updateT ***
-      !!---------------------------------------------
-      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
-      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
-      LOGICAL, iNTENT(in) :: before
-      !
+      !
+      ! ====>>>>>>>>>>    currently not used
+      !
+      !!----------------------------------------------------------------------
+      !!                      *** ROUTINE updateT ***
+      !!----------------------------------------------------------------------
+      INTEGER                                    , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2
+      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) ::   tabres
+      LOGICAL                                    , INTENT(in   ) ::   before
+      !!
       INTEGER :: ji,jj,jk
       REAL(wp) :: ztemp
-      !!---------------------------------------------
+      !!----------------------------------------------------------------------
 
       IF (before) THEN
@@ -587,16 +599,15 @@
    END SUBROUTINE update_scales
 
-# if defined key_zdftke
 
    SUBROUTINE updateEN( ptab, i1, i2, j1, j2, k1, k2, before )
-      !!---------------------------------------------
-      !!           *** ROUTINE updateen ***
-      !!---------------------------------------------
-      INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
-      LOGICAL, INTENT(in) :: before
-      !!---------------------------------------------
-      !
-      IF (before) THEN
+      !!----------------------------------------------------------------------
+      !!                      *** ROUTINE updateen ***
+      !!----------------------------------------------------------------------
+      INTEGER                               , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2
+      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
+      LOGICAL                               , INTENT(in   ) ::   before
+      !!----------------------------------------------------------------------
+      !
+      IF( before ) THEN
          ptab (i1:i2,j1:j2,k1:k2) = en(i1:i2,j1:j2,k1:k2)
       ELSE
@@ -608,16 +619,14 @@
 
    SUBROUTINE updateAVT( ptab, i1, i2, j1, j2, k1, k2, before )
-      !!---------------------------------------------
-      !!           *** ROUTINE updateavt ***
-      !!---------------------------------------------
-      INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
-      LOGICAL, INTENT(in) :: before
-      !!---------------------------------------------
-      !
-      IF (before) THEN
-         ptab (i1:i2,j1:j2,k1:k2) = avt_k(i1:i2,j1:j2,k1:k2)
-      ELSE
-         avt_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2) 
+      !!----------------------------------------------------------------------
+      !!                      *** ROUTINE updateavt ***
+      !!----------------------------------------------------------------------
+      INTEGER                               , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2
+      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
+      LOGICAL                               , INTENT(in   ) ::   before
+      !!----------------------------------------------------------------------
+      !
+      IF( before ) THEN   ;   ptab (i1:i2,j1:j2,k1:k2) = avt_k(i1:i2,j1:j2,k1:k2)
+      ELSE                ;   avt_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2) 
       ENDIF
       !
@@ -628,29 +637,27 @@
       !!---------------------------------------------
       !!           *** ROUTINE updateavm ***
-      !!---------------------------------------------
-      INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
-      LOGICAL, INTENT(in) :: before
-      !!---------------------------------------------
-      !
-      IF (before) THEN
-         ptab (i1:i2,j1:j2,k1:k2) = avm_k(i1:i2,j1:j2,k1:k2)
-      ELSE
-         avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2) 
+      !!----------------------------------------------------------------------
+      INTEGER                               , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2
+      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
+      LOGICAL                               , INTENT(in   ) ::   before
+      !!----------------------------------------------------------------------
+      !
+      IF( before ) THEN   ;   ptab (i1:i2,j1:j2,k1:k2) = avm_k(i1:i2,j1:j2,k1:k2)
+      ELSE                ;   avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2) 
       ENDIF
       !
    END SUBROUTINE updateAVM
 
-# endif /* key_zdftke */ 
-
 #else
+   !!----------------------------------------------------------------------
+   !!   Empty module                                          no AGRIF zoom
+   !!----------------------------------------------------------------------
 CONTAINS
    SUBROUTINE agrif_opa_update_empty
-      !!---------------------------------------------
-      !!   *** ROUTINE agrif_opa_update_empty ***
-      !!---------------------------------------------
       WRITE(*,*)  'agrif_opa_update : You should not have seen this print! error?'
    END SUBROUTINE agrif_opa_update_empty
 #endif
+
+   !!======================================================================
 END MODULE agrif_opa_update
 
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/NST_SRC/agrif_top_interp.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/NST_SRC/agrif_top_interp.F90	(revision 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/NST_SRC/agrif_top_interp.F90	(revision 7953)
@@ -1,4 +1,14 @@
 MODULE agrif_top_interp
+   !!======================================================================
+   !!                   ***  MODULE  agrif_top_interp  ***
+   !! AGRIF: interpolation package
+   !!======================================================================
+   !! History :  2.0  !  ??? 
+   !!----------------------------------------------------------------------
 #if defined key_agrif && defined key_top
+   !!----------------------------------------------------------------------
+   !!   'key_agrif'                                              AGRIF zoom
+   !!   'key_top'                                           on-line tracers
+   !!----------------------------------------------------------------------
    USE par_oce
    USE oce
@@ -8,4 +18,5 @@
    USE par_trc
    USE trc
+   !
    USE lib_mpp
    USE wrk_nemo  
@@ -16,7 +27,6 @@
    PUBLIC Agrif_trc, interptrn
 
-#  include "vectopt_loop_substitute.h90"
   !!----------------------------------------------------------------------
-   !! NEMO/NST 3.6 , NEMO Consortium (2010)
+   !! NEMO/NST 4.0 , NEMO Consortium (2017)
    !! $Id$
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
@@ -26,12 +36,12 @@
    SUBROUTINE Agrif_trc
       !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE Agrif_trc  ***
+      !!                   ***  ROUTINE Agrif_trc  ***
       !!----------------------------------------------------------------------
       !
       IF( Agrif_Root() )   RETURN
-
-      Agrif_SpecialValue    = 0.e0
+      !
+      Agrif_SpecialValue    = 0._wp
       Agrif_UseSpecialValue = .TRUE.
-
+      !
       CALL Agrif_Bc_variable( trn_id, procname=interptrn )
       Agrif_UseSpecialValue = .FALSE.
@@ -40,20 +50,21 @@
 
 
-   SUBROUTINE interptrn(ptab,i1,i2,j1,j2,k1,k2,n1,n2,before,nb,ndir)
-      !!---------------------------------------------
-      !!   *** ROUTINE interptrn ***
-      !!---------------------------------------------
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: ptab
-      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
-      LOGICAL, INTENT(in) :: before
-      INTEGER, INTENT(in) :: nb , ndir
-      !
-      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
-      INTEGER :: imin, imax, jmin, jmax
-      REAL(wp) ::   zrhox , zalpha1, zalpha2, zalpha3
-      REAL(wp) ::   zalpha4, zalpha5, zalpha6, zalpha7
-      LOGICAL :: western_side, eastern_side,northern_side,southern_side
-
-      IF (before) THEN         
+   SUBROUTINE interptrn( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir )
+      !!----------------------------------------------------------------------
+      !!                   *** ROUTINE interptrn ***
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) ::   ptab
+      INTEGER                                     , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2
+      LOGICAL                                     , INTENT(in   ) ::   before
+      INTEGER                                     , INTENT(in   ) ::   nb , ndir
+      !!
+      INTEGER ::   ji, jj, jk, jn   ! dummy loop indices
+      INTEGER ::   imin, imax, jmin, jmax
+      LOGICAL ::   western_side, eastern_side,northern_side,southern_side
+      REAL(wp)::   zrhox , zalpha1, zalpha2, zalpha3
+      REAL(wp)::   zalpha4, zalpha5, zalpha6, zalpha7
+      !!----------------------------------------------------------------------
+      !
+      IF( before ) THEN         
          ptab(i1:i2,j1:j2,k1:k2,n1:n2) = trn(i1:i2,j1:j2,k1:k2,n1:n2)
       ELSE
@@ -185,4 +196,7 @@
 
 #else
+   !!----------------------------------------------------------------------
+   !!   Empty module                                           no TOP AGRIF
+   !!----------------------------------------------------------------------
 CONTAINS
    SUBROUTINE Agrif_TOP_Interp_empty
@@ -193,3 +207,5 @@
    END SUBROUTINE Agrif_TOP_Interp_empty
 #endif
+
+   !!======================================================================
 END MODULE agrif_top_interp
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/NST_SRC/agrif_top_sponge.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/NST_SRC/agrif_top_sponge.F90	(revision 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/NST_SRC/agrif_top_sponge.F90	(revision 7953)
@@ -4,14 +4,13 @@
    !!======================================================================
    !!                ***  MODULE agrif_top_sponge  ***
-   !! AGRIF :   define in memory AGRIF variables for sea-ice
+   !! AGRIF :   TOP sponge layer
    !!======================================================================
    !! History :  2.0  ! 2006-08  (R. Benshila, L. Debreu)  Original code
    !!----------------------------------------------------------------------
-
+#if defined key_agrif && defined key_top
    !!----------------------------------------------------------------------
    !!   Agrif_Sponge_trc : 
    !!   interptrn_sponge :  
    !!----------------------------------------------------------------------
-#if defined key_agrif && defined key_top
    USE par_oce
    USE par_trc
@@ -32,5 +31,5 @@
 
    !!----------------------------------------------------------------------
-   !! NEMO/NST 3.7 , NEMO Consortium (2015)
+   !! NEMO/NST 4.0 , NEMO Consortium (2017)
    !! $Id$
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
@@ -42,5 +41,5 @@
       !!                   *** ROUTINE Agrif_Sponge_Trc ***
       !!----------------------------------------------------------------------
-      REAL(wp) ::   timecoeff
+      REAL(wp) ::   timecoeff   ! local scalar
       !!----------------------------------------------------------------------
       !
@@ -107,5 +106,7 @@
 
 #else
-
+   !!----------------------------------------------------------------------
+   !!   Empty module                                           no TOP AGRIF
+   !!----------------------------------------------------------------------
 CONTAINS
    SUBROUTINE agrif_top_sponge_empty
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/NST_SRC/agrif_top_update.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/NST_SRC/agrif_top_update.F90	(revision 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/NST_SRC/agrif_top_update.F90	(revision 7953)
@@ -6,15 +6,19 @@
    !!                ***  MODULE agrif_top_update  ***
    !! AGRIF :   
-   !!----------------------------------------------------------------------
+   !!======================================================================
    !! History :  
    !!----------------------------------------------------------------------
-
 #if defined key_agrif && defined key_top
+   !!----------------------------------------------------------------------
+   !!   'key_agrif'                                              AGRIF zoom
+   !!   'key_TOP'                                           on-line tracers
+   !!----------------------------------------------------------------------
    USE par_oce
    USE oce
+   USE dom_oce
+   USE agrif_oce
    USE par_trc
    USE trc
-   USE dom_oce
-   USE agrif_oce
+   !
    USE wrk_nemo  
 
@@ -27,5 +31,5 @@
 
    !!----------------------------------------------------------------------
-   !! NEMO/NST 3.7 , NEMO Consortium (2015)
+   !! NEMO/NST 4.0 , NEMO Consortium (2017)
    !! $Id$
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
@@ -112,9 +116,9 @@
 
 #else
+   !!----------------------------------------------------------------------
+   !!   Empty module                                           no TOP AGRIF
+   !!----------------------------------------------------------------------
 CONTAINS
    SUBROUTINE agrif_top_update_empty
-      !!---------------------------------------------
-      !!   *** ROUTINE agrif_Top_update_empty ***
-      !!---------------------------------------------
       WRITE(*,*)  'agrif_top_update : You should not have seen this print! error?'
    END SUBROUTINE agrif_top_update_empty
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/NST_SRC/agrif_user.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/NST_SRC/agrif_user.F90	(revision 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/NST_SRC/agrif_user.F90	(revision 7953)
@@ -1,5 +1,5 @@
 #if defined key_agrif
 !!----------------------------------------------------------------------
-!! NEMO/NST 3.7 , NEMO Consortium (2016)
+!! NEMO/NST 4.0 , NEMO Consortium (2017)
 !! $Id$
 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
@@ -107,4 +107,5 @@
    !!
    IMPLICIT NONE
+   !
    !!----------------------------------------------------------------------
    !
@@ -125,5 +126,5 @@
    USE par_oce       
    USE oce
-   !!
+   !
    IMPLICIT NONE
    !!----------------------------------------------------------------------
@@ -136,20 +137,20 @@
    ! 2. Type of interpolation
    !-------------------------
-   CALL Agrif_Set_bcinterp(e1u_id,interp1=Agrif_linear,interp2=AGRIF_ppm)
-   CALL Agrif_Set_bcinterp(e2v_id,interp1=AGRIF_ppm,interp2=Agrif_linear)
+   CALL Agrif_Set_bcinterp( e1u_id, interp1=Agrif_linear, interp2=AGRIF_ppm    )
+   CALL Agrif_Set_bcinterp( e2v_id, interp1=AGRIF_ppm   , interp2=Agrif_linear )
 
    ! 3. Location of interpolation
    !-----------------------------
-   CALL Agrif_Set_bc(e1u_id,(/0,0/))
-   CALL Agrif_Set_bc(e2v_id,(/0,0/))
+   CALL Agrif_Set_bc( e1u_id, (/0,0/) )
+   CALL Agrif_Set_bc( e2v_id, (/0,0/) )
 
    ! 5. Update type
    !--------------- 
-   CALL Agrif_Set_Updatetype(e1u_id,update1 = Agrif_Update_Copy, update2=Agrif_Update_Average)
-   CALL Agrif_Set_Updatetype(e2v_id,update1 = Agrif_Update_Average, update2=Agrif_Update_Copy)
+   CALL Agrif_Set_Updatetype( e1u_id, update1=Agrif_Update_Copy   , update2=Agrif_Update_Average )
+   CALL Agrif_Set_Updatetype( e2v_id, update1=Agrif_Update_Average, update2=Agrif_Update_Copy    )
 
 ! High order updates
-!   CALL Agrif_Set_Updatetype(e1u_id,update1 = Agrif_Update_Average,            update2=Agrif_Update_Full_Weighting)
-!   CALL Agrif_Set_Updatetype(e2v_id,update1 = Agrif_Update_Full_Weighting,     update2=Agrif_Update_Average)
+!   CALL Agrif_Set_Updatetype( e1u_id, update1=Agrif_Update_Average       , update2=Agrif_Update_Full_Weighting )
+!   CALL Agrif_Set_Updatetype( e2v_id, update1=Agrif_Update_Full_Weighting, update2=Agrif_Update_Average        )
     !
 END SUBROUTINE agrif_declare_var_dom
@@ -165,5 +166,7 @@
    USE oce 
    USE dom_oce
+   USE zdf_oce
    USE nemogcm
+   !
    USE lib_mpp
    USE in_out_manager
@@ -171,5 +174,5 @@
    USE agrif_opa_interp
    USE agrif_opa_sponge
-   !!
+   !
    IMPLICIT NONE
    !
@@ -184,5 +187,5 @@
    ! 2. First interpolations of potentially non zero fields
    !-------------------------------------------------------
-   Agrif_SpecialValue=0.
+   Agrif_SpecialValue    = 0._wp
    Agrif_UseSpecialValue = .TRUE.
    CALL Agrif_Bc_variable(tsn_id,calledweight=1.,procname=interptsn)
@@ -319,7 +322,5 @@
    ENDIF
    !
-# if defined key_zdftke
-   CALL Agrif_Update_tke(0)
-# endif
+   IF( ln_zdftke )   CALL Agrif_Update_tke( 0 )
    !
    Agrif_UseSpecialValueInUpdate = .FALSE.
@@ -337,7 +338,8 @@
    !!----------------------------------------------------------------------
    USE agrif_util
-   USE par_oce       !   ONLY : jpts
+   USE agrif_oce
+   USE par_oce       ! ocean parameters
+   USE zdf_oce       ! vertical physics
    USE oce
-   USE agrif_oce
    !!
    IMPLICIT NONE
@@ -371,9 +373,9 @@
    CALL agrif_declare_variable((/2,2/),(/3,3/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),sshn_id)
 
-# if defined key_zdftke
-   CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/jpi,jpj,jpk/), en_id)
-   CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/jpi,jpj,jpk/),avt_id)
-   CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/jpi,jpj,jpk/),avm_id)
-# endif
+   IF( ln_zdftke ) THEN
+      CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/jpi,jpj,jpk/), en_id)
+      CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/jpi,jpj,jpk/),avt_id)
+      CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/jpi,jpj,jpk/),avm_id)
+   ENDIF
 
    ! 2. Type of interpolation
@@ -400,8 +402,5 @@
    CALL Agrif_Set_bcinterp(vmsk_id,interp=AGRIF_constant)
 
-# if defined key_zdftke
-   CALL Agrif_Set_bcinterp(avm_id ,interp=AGRIF_linear)
-# endif
-
+   IF( ln_zdftke )   CALL Agrif_Set_bcinterp( avm_id, interp=AGRIF_linear )
 
    ! 3. Location of interpolation
@@ -418,9 +417,9 @@
    CALL Agrif_Set_bc(vn_sponge_id ,(/-nn_sponge_len*Agrif_irhox()-1,0/))
 
-   CALL Agrif_Set_bc(sshn_id,(/0,0/))
-   CALL Agrif_Set_bc(unb_id ,(/0,0/))
-   CALL Agrif_Set_bc(vnb_id ,(/0,0/))
-   CALL Agrif_Set_bc(ub2b_interp_id,(/0,0/))
-   CALL Agrif_Set_bc(vb2b_interp_id,(/0,0/))
+   CALL Agrif_Set_bc( sshn_id       , (/0,0/) )
+   CALL Agrif_Set_bc( unb_id        , (/0,0/) )
+   CALL Agrif_Set_bc( vnb_id        , (/0,0/) )
+   CALL Agrif_Set_bc( ub2b_interp_id, (/0,0/) )
+   CALL Agrif_Set_bc( vb2b_interp_id, (/0,0/) )
 
    CALL Agrif_Set_bc(e3t_id,(/-2*Agrif_irhox()-1,0/))   ! if west and rhox=3: column 2 to 9
@@ -428,7 +427,5 @@
    CALL Agrif_Set_bc(vmsk_id,(/0,0/))
 
-# if defined key_zdftke
-   CALL Agrif_Set_bc(avm_id ,(/0,1/))
-# endif
+   IF( ln_zdftke )   CALL Agrif_Set_bc( avm_id, (/0,1/) )
 
    ! 5. Update type
@@ -446,9 +443,9 @@
    CALL Agrif_Set_Updatetype(vb2b_update_id,update1 = Agrif_Update_Average, update2 = Agrif_Update_Copy)
 
-# if defined key_zdftke
-   CALL Agrif_Set_Updatetype( en_id, update = AGRIF_Update_Average)
-   CALL Agrif_Set_Updatetype(avt_id, update = AGRIF_Update_Average)
-   CALL Agrif_Set_Updatetype(avm_id, update = AGRIF_Update_Average)
-# endif
+   IF( ln_zdftke) THEN
+      CALL Agrif_Set_Updatetype( en_id, update = AGRIF_Update_Average)
+      CALL Agrif_Set_Updatetype(avt_id, update = AGRIF_Update_Average)
+      CALL Agrif_Set_Updatetype(avm_id, update = AGRIF_Update_Average)
+   ENDIF
 
 ! High order updates
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OFF_SRC/nemogcm.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OFF_SRC/nemogcm.F90	(revision 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OFF_SRC/nemogcm.F90	(revision 7953)
@@ -28,11 +28,10 @@
    USE trabbl          ! bottom boundary layer          (tra_bbl_init routine)
    USE traldf          ! lateral physics                (tra_ldf_init routine)
-   USE zdfini          ! vertical physics: initialization
-   USE sbcmod          ! surface boundary condition       (sbc_init     routine)
-   USE phycst          ! physical constant                  (par_cst routine)
+   USE sbcmod          ! surface boundary condition     (sbc_init     routine)
+   USE phycst          ! physical constant                   (par_cst routine)
    USE dtadyn          ! Lecture and Interpolation of the dynamical fields
    USE trcini          ! Initilization of the passive tracers
-   USE daymod          ! calendar                         (day     routine)
-   USE trcstp          ! passive tracer time-stepping      (trc_stp routine)
+   USE daymod          ! calendar                            (day     routine)
+   USE trcstp          ! passive tracer time-stepping        (trc_stp routine)
    USE dtadyn          ! Lecture and interpolation of the dynamical fields
    !              ! Passive tracers needs
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90	(revision 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90	(revision 7953)
@@ -28,7 +28,5 @@
    USE ldfslp             ! Lateral diffusion: slopes of neutral surfaces
    USE tradmp             ! Tracer damping
-#if defined key_zdftke
    USE zdftke             ! TKE vertical physics
-#endif
    USE eosbn2             ! Equation of state (eos_bn2 routine)
    USE zdfmxl             ! Mixed layer depth
@@ -94,10 +92,8 @@
             IF( nitbkg_r == nit000 - 1 ) THEN      ! Treat special case when nitbkg = 0
                zdate = REAL( ndastp )
-#if defined key_zdftke
-               ! lk_zdftke=T :   Read turbulent kinetic energy ( en )
-               IF(lwp) WRITE(numout,*) ' Reading TKE (en) from restart...'
-               CALL tke_rst( nit000, 'READ' )               ! lk_zdftke=T :   Read turbulent kinetic energy ( en )
-
-#endif
+               IF( ln_zdftke ) THEN                   ! read turbulent kinetic energy ( en )
+                  IF(lwp) WRITE(numout,*) ' Reading TKE (en) from restart...'
+                  CALL tke_rst( nit000, 'READ' )
+               ENDIF
             ELSE
                zdate = REAL( ndastp )
@@ -111,7 +107,5 @@
             CALL iom_rstput( kt, nitbkg_r, inum, 'sn'     , tsn(:,:,:,jp_sal) )
             CALL iom_rstput( kt, nitbkg_r, inum, 'sshn'   , sshn              )
-#if defined key_zdftke
-            CALL iom_rstput( kt, nitbkg_r, inum, 'en'     , en                )
-#endif
+            IF( ln_zdftke )   CALL iom_rstput( kt, nitbkg_r, inum, 'en'     , en                )
             !
             CALL iom_close( inum )
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/C1D/step_c1d.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/C1D/step_c1d.F90	(revision 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/C1D/step_c1d.F90	(revision 7953)
@@ -74,28 +74,7 @@
                          CALL bn2( tsb, rab_b, rn2b ) ! before Brunt-Vaisala frequency
                          CALL bn2( tsn, rab_n, rn2  ) ! now    Brunt-Vaisala frequency
-      !  VERTICAL PHYSICS   
-                         CALL zdf_bfr( kstp )         ! bottom friction
-      !                                               ! Vertical eddy viscosity and diffusivity coefficients
-      IF( lk_zdfric  )   CALL zdf_ric( kstp )            ! Richardson number dependent Kz
-      IF( lk_zdftke  )   CALL zdf_tke( kstp )            ! TKE closure scheme for Kz
-      IF( lk_zdfgls  )   CALL zdf_gls( kstp )            ! GLS closure scheme for Kz
-      IF( lk_zdfcst  )   THEN                            ! Constant Kz (reset avt, avm[uv] to the background value)
-         avt (:,:,:) = rn_avt0 * tmask(:,:,:)
-         avmu(:,:,:) = rn_avm0 * umask(:,:,:)
-         avmv(:,:,:) = rn_avm0 * vmask(:,:,:)
-      ENDIF
-
-      IF( ln_rnf_mouth ) THEN                         ! increase diffusivity at rivers mouths
-         DO jk = 2, nkrnf   ;   avt(:,:,jk) = avt(:,:,jk) + 2.e0 * rn_avt_rnf * rnfmsk(:,:)   ;   END DO
-      ENDIF
-      IF( ln_zdfevd  )   CALL zdf_evd( kstp )         ! enhanced vertical eddy diffusivity
-      IF( lk_zdftmx  )   CALL zdf_tmx( kstp )         ! tidal vertical mixing
-      IF( ln_zdfddm  )   CALL zdf_ddm( kstp )         ! double diffusive mixing
-                         CALL zdf_mxl( kstp )         ! mixed layer depth
-
-                                                      ! write tke information in the restart file
-      IF( lrst_oce .AND. lk_zdftke )   CALL tke_rst( kstp, 'WRITE' )
-                                                      ! write gls information in the restart file
-      IF( lrst_oce .AND. lk_zdfgls )   CALL gls_rst( kstp, 'WRITE' )
+      
+      !  VERTICAL PHYSICS
+                         CALL zdf_phy( kstp )         ! vertical physics update (bfr, avt, avs, avm + MLD)
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/CRS/crs.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/CRS/crs.F90	(revision 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/CRS/crs.F90	(revision 7953)
@@ -151,8 +151,5 @@
 
       ! Vertical diffusion
-      REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)  ::  avt_crs           !: vert. diffusivity coef. [m2/s] at w-point for temp  
-# if defined key_zdfddm
       REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)  ::  avs_crs           !: salinity vertical diffusivity coeff. [m2/s] at w-point
-# endif
 
       ! Mixing and Mixed Layer Depth
@@ -239,21 +236,16 @@
          &     fr_i_crs(jpi_crs,jpj_crs), sfx_crs(jpi_crs ,jpj_crs),  STAT=ierr(12)  )
 
-     ALLOCATE( tsn_crs(jpi_crs,jpj_crs,jpk,jpts), avt_crs(jpi_crs,jpj_crs,jpk),    &
-# if defined key_zdfddm
-         &      avs_crs(jpi_crs,jpj_crs,jpk),    &
-# endif
-         &      STAT=ierr(13) )
+     ALLOCATE( tsn_crs(jpi_crs,jpj_crs,jpk,jpts), avs_crs(jpi_crs,jpj_crs,jpk), STAT=ierr(13) )
 
       ALLOCATE( nmln_crs(jpi_crs,jpj_crs) , hmld_crs(jpi_crs,jpj_crs) , &
          &      hmlp_crs(jpi_crs,jpj_crs) , hmlpt_crs(jpi_crs,jpj_crs) , STAT=ierr(14) )
          
-      ALLOCATE( nimppt_crs(jpnij) , nlcit_crs(jpnij) , nldit_crs(jpnij) , nleit_crs(jpnij), &
-       &  nimppt_full(jpnij) , nlcit_full(jpnij) , nldit_full(jpnij) , nleit_full(jpnij),   &
-                njmppt_crs(jpnij) , nlcjt_crs(jpnij) , nldjt_crs(jpnij) , nlejt_crs(jpnij), &
-       &  njmppt_full(jpnij) , nlcjt_full(jpnij) , nldjt_full(jpnij) , nlejt_full(jpnij)  , STAT=ierr(15) )
-
-         
+      ALLOCATE( nimppt_crs (jpnij) , nlcit_crs (jpnij) , nldit_crs (jpnij) , nleit_crs (jpnij),   &
+         &      nimppt_full(jpnij) , nlcit_full(jpnij) , nldit_full(jpnij) , nleit_full(jpnij),   &
+                njmppt_crs (jpnij) , nlcjt_crs (jpnij) , nldjt_crs (jpnij) , nlejt_crs (jpnij),   &
+         &      njmppt_full(jpnij) , nlcjt_full(jpnij) , nldjt_full(jpnij) , nlejt_full(jpnij)  , STAT=ierr(15) )
+   
       crs_dom_alloc = MAXVAL(ierr)
-
+      !
    END FUNCTION crs_dom_alloc
 
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/CRS/crsfld.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/CRS/crsfld.F90	(revision 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/CRS/crsfld.F90	(revision 7953)
@@ -84,5 +84,5 @@
          vn_crs   (:,:,:  ) = 0._wp    ! v-velocity
          wn_crs   (:,:,:  ) = 0._wp    ! w
-         avt_crs  (:,:,:  ) = 0._wp    ! avt
+         avs_crs  (:,:,:  ) = 0._wp    ! avt
          hdivn_crs(:,:,:  ) = 0._wp    ! hdiv
          rke_crs  (:,:,:  ) = 0._wp    ! rke
@@ -200,12 +200,12 @@
       SELECT CASE ( nn_crs_kz )
          CASE ( 0 )
-            CALL crs_dom_ope( avt, 'VOL', 'W', tmask, avt_crs, p_e12=e1e2t, p_e3=ze3w, psgn=1.0 )
+            CALL crs_dom_ope( avt, 'VOL', 'W', tmask, avs_crs, p_e12=e1e2t, p_e3=ze3w, psgn=1.0 )
          CASE ( 1 )
-            CALL crs_dom_ope( avt, 'MAX', 'W', tmask, avt_crs, p_e12=e1e2t, p_e3=ze3w, psgn=1.0 )
+            CALL crs_dom_ope( avt, 'MAX', 'W', tmask, avs_crs, p_e12=e1e2t, p_e3=ze3w, psgn=1.0 )
          CASE ( 2 )
-            CALL crs_dom_ope( avt, 'MIN', 'W', tmask, avt_crs, p_e12=e1e2t, p_e3=ze3w, psgn=1.0 )
+            CALL crs_dom_ope( avt, 'MIN', 'W', tmask, avs_crs, p_e12=e1e2t, p_e3=ze3w, psgn=1.0 )
       END SELECT
       !
-      CALL iom_put( "avt", avt_crs )   !  Kz
+      CALL iom_put( "avt", avs_crs )   !  Kz
       
       !  sbc fields  
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 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90	(revision 7953)
@@ -8,37 +8,25 @@
    USE oce             ! ocean dynamics and tracers variables
    USE dom_oce         ! ocean space and time domain
+   USE zdf_oce         ! ocean vertical physics    
+   USE zdfgls, ONLY: mxln
    USE in_out_manager  ! I/O units
    USE iom             ! I/0 library
-   USE wrk_nemo        ! working arrays
-#if defined key_zdftke 
-   USE zdf_oce, ONLY: en
-#endif
-   USE zdf_oce, ONLY: avt, avm
-#if defined key_zdfgls
-   USE zdf_oce, ONLY: en
-   USE zdfgls, ONLY: mxln
-#endif
+   USE wrk_nemo        ! work arrays
 
    IMPLICIT NONE
    PRIVATE
 
-   LOGICAL , PUBLIC ::   ln_dia25h     !:  25h mean output
    PUBLIC   dia_25h_init               ! routine called by nemogcm.F90
    PUBLIC   dia_25h                    ! routine called by diawri.F90
 
-  !! * variables for calculating 25-hourly means
-   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   tn_25h  , sn_25h
-   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:)   ::   sshn_25h 
-   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   un_25h  , vn_25h  , wn_25h
-   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   avt_25h , avm_25h
-#if defined key_zdfgls || key_zdftke
-   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   en_25h
-#endif
-#if defined key_zdfgls 
-   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   rmxln_25h
-#endif
-   INTEGER, SAVE :: cnt_25h     ! Counter for 25 hour means
-
-
+   LOGICAL, PUBLIC ::   ln_dia25h      !:  25h mean output
+
+   ! variables for calculating 25-hourly means
+   INTEGER , SAVE ::   cnt_25h     ! Counter for 25 hour means
+   REAL(wp), SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   tn_25h  , sn_25h
+   REAL(wp), SAVE, ALLOCATABLE,   DIMENSION(:,:)   ::   sshn_25h 
+   REAL(wp), SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   un_25h  , vn_25h  , wn_25h
+   REAL(wp), SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   avt_25h , avm_25h
+   REAL(wp), SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   en_25h  , rmxln_25h
 
    !!----------------------------------------------------------------------
@@ -56,8 +44,5 @@
       !!        
       !! ** Method : Read namelist
-      !!   History
-      !!   3.6  !  08-14  (E. O'Dea) Routine to initialize dia_25h
       !!---------------------------------------------------------------------------
-      !!
       INTEGER ::   ios                 ! Local integer output status for namelist read
       INTEGER ::   ierror              ! Local integer for memory allocation
@@ -79,6 +64,6 @@
          WRITE(numout,*) 'dia_25h_init : Output 25 hour mean diagnostics'
          WRITE(numout,*) '~~~~~~~~~~~~'
-         WRITE(numout,*) 'Namelist nam_dia25h : set 25h outputs '
-         WRITE(numout,*) 'Switch for 25h diagnostics (T) or not (F)  ln_dia25h  = ', ln_dia25h
+         WRITE(numout,*) '   Namelist nam_dia25h : set 25h outputs '
+         WRITE(numout,*) '      Switch for 25h diagnostics (T) or not (F)  ln_dia25h  = ', ln_dia25h
       ENDIF
       IF( .NOT. ln_dia25h )   RETURN
@@ -86,47 +71,22 @@
       ! 1 - Allocate memory !
       ! ------------------- !
-      ALLOCATE( tn_25h(jpi,jpj,jpk), STAT=ierror )
+      !                                ! ocean arrays
+      ALLOCATE( tn_25h (jpi,jpj,jpk), sn_25h (jpi,jpj,jpk), sshn_25h(jpi,jpj)  ,     &
+         &      un_25h (jpi,jpj,jpk), vn_25h (jpi,jpj,jpk), wn_25h(jpi,jpj,jpk),     &
+         &      avt_25h(jpi,jpj,jpk), avm_25h(jpi,jpj,jpk),                      STAT=ierror )
       IF( ierror > 0 ) THEN
-         CALL ctl_stop( 'dia_25h: unable to allocate tn_25h' )   ;   RETURN
-      ENDIF
-      ALLOCATE( sn_25h(jpi,jpj,jpk), STAT=ierror )
-      IF( ierror > 0 ) THEN
-         CALL ctl_stop( 'dia_25h: unable to allocate sn_25h' )   ;   RETURN
-      ENDIF
-      ALLOCATE( un_25h(jpi,jpj,jpk), STAT=ierror )
-      IF( ierror > 0 ) THEN
-         CALL ctl_stop( 'dia_25h: unable to allocate un_25h' )   ;   RETURN
-      ENDIF
-      ALLOCATE( vn_25h(jpi,jpj,jpk), STAT=ierror )
-      IF( ierror > 0 ) THEN
-         CALL ctl_stop( 'dia_25h: unable to allocate vn_25h' )   ;   RETURN
-      ENDIF
-      ALLOCATE( wn_25h(jpi,jpj,jpk), STAT=ierror )
-      IF( ierror > 0 ) THEN
-         CALL ctl_stop( 'dia_25h: unable to allocate wn_25h' )   ;   RETURN
-      ENDIF
-      ALLOCATE( avt_25h(jpi,jpj,jpk), STAT=ierror )
-      IF( ierror > 0 ) THEN
-         CALL ctl_stop( 'dia_25h: unable to allocate avt_25h' )   ;   RETURN
-      ENDIF
-      ALLOCATE( avm_25h(jpi,jpj,jpk), STAT=ierror )
-      IF( ierror > 0 ) THEN
-         CALL ctl_stop( 'dia_25h: unable to allocate avm_25h' )   ;   RETURN
-      ENDIF
-# if defined key_zdfgls || defined key_zdftke
-      ALLOCATE( en_25h(jpi,jpj,jpk), STAT=ierror )
-      IF( ierror > 0 ) THEN
-         CALL ctl_stop( 'dia_25h: unable to allocate en_25h' )   ;   RETURN
-      ENDIF
-#endif
-# if defined key_zdfgls 
-      ALLOCATE( rmxln_25h(jpi,jpj,jpk), STAT=ierror )
-      IF( ierror > 0 ) THEN
-         CALL ctl_stop( 'dia_25h: unable to allocate rmxln_25h' )   ;   RETURN
-      ENDIF
-#endif
-      ALLOCATE( sshn_25h(jpi,jpj), STAT=ierror )
-      IF( ierror > 0 ) THEN
-         CALL ctl_stop( 'dia_25h: unable to allocate sshn_25h' )   ;   RETURN
+         CALL ctl_stop( 'dia_25h: unable to allocate ocean arrays' )   ;   RETURN
+      ENDIF
+      IF( ln_zdftke ) THEN             ! TKE physics
+         ALLOCATE( en_25h(jpi,jpj,jpk), STAT=ierror )
+         IF( ierror > 0 ) THEN
+            CALL ctl_stop( 'dia_25h: unable to allocate en_25h' )   ;   RETURN
+         ENDIF
+      ENDIF
+      IF( ln_zdfgls ) THEN             ! GLS physics
+         ALLOCATE( en_25h(jpi,jpj,jpk), rmxln_25h(jpi,jpj,jpk), STAT=ierror )
+         IF( ierror > 0 ) THEN
+            CALL ctl_stop( 'dia_25h: unable to allocate en_25h and rmxln_25h' )   ;   RETURN
+         ENDIF
       ENDIF
       ! ------------------------- !
@@ -142,19 +102,15 @@
       avt_25h(:,:,:) = avt(:,:,:)
       avm_25h(:,:,:) = avm(:,:,:)
-# if defined key_zdfgls || defined key_zdftke
+      IF( ln_zdftke ) THEN
          en_25h(:,:,:) = en(:,:,:)
-#endif
-# if defined key_zdfgls
+      ENDIF
+      IF( ln_zdfgls ) THEN
+         en_25h(:,:,:) = en(:,:,:)
          rmxln_25h(:,:,:) = mxln(:,:,:)
-#endif
+      ENDIF
 #if defined key_lim3 || defined key_lim2
          CALL ctl_stop('STOP', 'dia_25h not setup yet to do tidemean ice')
 #endif 
-
-      ! -------------------------- !
-      ! 3 - Return to dia_wri      !
-      ! -------------------------- !
-
-
+      !
    END SUBROUTINE dia_25h_init
 
@@ -164,36 +120,18 @@
       !!                 ***  ROUTINE dia_25h  ***
       !!         
-      !!
-      !!--------------------------------------------------------------------
-      !!                   
       !! ** Purpose :   Write diagnostics with M2/S2 tide removed
       !!
-      !! ** Method  :   
-      !!      25hr mean outputs for shelf seas
+      !! ** Method  :   25hr mean outputs for shelf seas
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
       !!
-      !! History :
-      !!   ?.0  !  07-04  (A. Hines) New routine, developed from dia_wri_foam
-      !!   3.4  !  02-13  (J. Siddorn) Routine taken from old dia_wri_foam
-      !!   3.6  !  08-14  (E. O'Dea) adapted for VN3.6
-      !!----------------------------------------------------------------------
-      !! * Modules used
-
-      IMPLICIT NONE
-
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
-
-
-      !! * Local declarations
       INTEGER ::   ji, jj, jk
-
+      INTEGER                          ::   iyear0, nimonth0,iday0            ! start year,imonth,day
       LOGICAL ::   ll_print = .FALSE.    ! =T print and flush numout
-      REAL(wp)                         ::   zsto, zout, zmax, zjulian, zmdi       ! temporary reals
-      INTEGER                          ::   i_steps                               ! no of timesteps per hour
-      REAL(wp), DIMENSION(jpi,jpj    ) ::   zw2d, un_dm, vn_dm                    ! temporary workspace
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zw3d                                  ! temporary workspace
-      REAL(wp), DIMENSION(jpi,jpj,3)   ::   zwtmb                                 ! temporary workspace
-      INTEGER                          ::   iyear0, nimonth0,iday0                ! start year,imonth,day
-
+      REAL(wp)                         ::   zsto, zout, zmax, zjulian, zmdi   ! local scalars
+      INTEGER                          ::   i_steps                           ! no of timesteps per hour
+      REAL(wp), DIMENSION(jpi,jpj    ) ::   zw2d, un_dm, vn_dm                ! workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zw3d                              ! workspace
+      REAL(wp), DIMENSION(jpi,jpj,3)   ::   zwtmb                             ! workspace
       !!----------------------------------------------------------------------
 
@@ -207,14 +145,9 @@
       ENDIF
 
-#if defined key_lim3 || defined key_lim2
-      CALL ctl_stop('STOP', 'dia_wri_tide not setup yet to do tidemean ice')
-#endif
-
       ! local variable for debugging
       ll_print = ll_print .AND. lwp
 
-      ! Sum of 25 hourly instantaneous values to give a 25h mean from 24hours
-      ! every day
-      IF( MOD( kt, i_steps ) == 0  .and. kt .ne. nn_it000 ) THEN
+      ! Sum of 25 hourly instantaneous values to give a 25h mean from 24hours every day
+      IF( MOD( kt, i_steps ) == 0  .AND. kt /= nn_it000 ) THEN
 
          IF (lwp) THEN
@@ -231,94 +164,96 @@
          avt_25h(:,:,:)       = avt_25h(:,:,:) + avt(:,:,:)
          avm_25h(:,:,:)       = avm_25h(:,:,:) + avm(:,:,:)
-# if defined key_zdfgls || defined key_zdftke
-         en_25h(:,:,:)        = en_25h(:,:,:) + en(:,:,:)
-#endif
-# if defined key_zdfgls
-         rmxln_25h(:,:,:)      = rmxln_25h(:,:,:) + mxln(:,:,:)
-#endif
+         IF( ln_zdftke ) THEN
+            en_25h(:,:,:)        = en_25h(:,:,:) + en(:,:,:)
+         ENDIF
+         IF( ln_zdfgls ) THEN
+            en_25h(:,:,:)        = en_25h(:,:,:) + en(:,:,:)
+            rmxln_25h(:,:,:)      = rmxln_25h(:,:,:) + mxln(:,:,:)
+         ENDIF
          cnt_25h = cnt_25h + 1
-
+         !
          IF (lwp) THEN
             WRITE(numout,*) 'dia_tide : Summed the following number of hourly values so far',cnt_25h
          ENDIF
-
+         !
       ENDIF ! MOD( kt, i_steps ) == 0
 
-         ! Write data for 25 hour mean output streams
-      IF( cnt_25h .EQ. 25 .AND.  MOD( kt, i_steps*24) == 0 .AND. kt .NE. nn_it000 ) THEN
-
-            IF(lwp) THEN
-               WRITE(numout,*) 'dia_wri_tide : Writing 25 hour mean tide diagnostics at timestep', kt
-               WRITE(numout,*) '~~~~~~~~~~~~ '
-            ENDIF
-
-            tn_25h(:,:,:)        = tn_25h(:,:,:) / 25.0_wp
-            sn_25h(:,:,:)        = sn_25h(:,:,:) / 25.0_wp
-            sshn_25h(:,:)        = sshn_25h(:,:) / 25.0_wp
-            un_25h(:,:,:)        = un_25h(:,:,:) / 25.0_wp
-            vn_25h(:,:,:)        = vn_25h(:,:,:) / 25.0_wp
-            wn_25h(:,:,:)        = wn_25h(:,:,:) / 25.0_wp
-            avt_25h(:,:,:)       = avt_25h(:,:,:) / 25.0_wp
-            avm_25h(:,:,:)       = avm_25h(:,:,:) / 25.0_wp
-# if defined key_zdfgls || defined key_zdftke
+      ! Write data for 25 hour mean output streams
+      IF( cnt_25h == 25 .AND.  MOD( kt, i_steps*24) == 0 .AND. kt /= nn_it000 ) THEN
+         !
+         IF(lwp) THEN
+            WRITE(numout,*) 'dia_wri_tide : Writing 25 hour mean tide diagnostics at timestep', kt
+            WRITE(numout,*) '~~~~~~~~~~~~ '
+         ENDIF
+         !
+         tn_25h(:,:,:)        = tn_25h(:,:,:) / 25.0_wp
+         sn_25h(:,:,:)        = sn_25h(:,:,:) / 25.0_wp
+         sshn_25h(:,:)        = sshn_25h(:,:) / 25.0_wp
+         un_25h(:,:,:)        = un_25h(:,:,:) / 25.0_wp
+         vn_25h(:,:,:)        = vn_25h(:,:,:) / 25.0_wp
+         wn_25h(:,:,:)        = wn_25h(:,:,:) / 25.0_wp
+         avt_25h(:,:,:)       = avt_25h(:,:,:) / 25.0_wp
+         avm_25h(:,:,:)       = avm_25h(:,:,:) / 25.0_wp
+         IF( ln_zdftke ) THEN
             en_25h(:,:,:)        = en_25h(:,:,:) / 25.0_wp
-#endif
-# if defined key_zdfgls
+         ENDIF
+         IF( ln_zdfgls ) THEN
+            en_25h(:,:,:)        = en_25h(:,:,:) / 25.0_wp
             rmxln_25h(:,:,:)       = rmxln_25h(:,:,:) / 25.0_wp
-#endif
-
-            IF (lwp)  WRITE(numout,*) 'dia_wri_tide : Mean calculated by dividing 25 hour sums and writing output'
-            zmdi=1.e+20 !missing data indicator for masking
-            ! write tracers (instantaneous)
-            zw3d(:,:,:) = tn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
-            CALL iom_put("temper25h", zw3d)   ! potential temperature
-            zw3d(:,:,:) = sn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
-            CALL iom_put( "salin25h", zw3d  )   ! salinity
-            zw2d(:,:) = sshn_25h(:,:)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1))
-            CALL iom_put( "ssh25h", zw2d )   ! sea surface 
-
-
-            ! Write velocities (instantaneous)
-            zw3d(:,:,:) = un_25h(:,:,:)*umask(:,:,:) + zmdi*(1.0-umask(:,:,:))
-            CALL iom_put("vozocrtx25h", zw3d)    ! i-current
-            zw3d(:,:,:) = vn_25h(:,:,:)*vmask(:,:,:) + zmdi*(1.0-vmask(:,:,:))
-            CALL iom_put("vomecrty25h", zw3d  )   ! j-current
-
-            zw3d(:,:,:) = wn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
-            CALL iom_put("vomecrtz25h", zw3d )   ! k-current
-            zw3d(:,:,:) = avt_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
-            CALL iom_put("avt25h", zw3d )   ! diffusivity
-            zw3d(:,:,:) = avm_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
-            CALL iom_put("avm25h", zw3d)   ! viscosity
-#if defined key_zdftke || defined key_zdfgls 
+         ENDIF
+         !
+         IF(lwp)  WRITE(numout,*) 'dia_wri_tide : Mean calculated by dividing 25 hour sums and writing output'
+         zmdi=1.e+20 !missing data indicator for masking
+         ! write tracers (instantaneous)
+         zw3d(:,:,:) = tn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+         CALL iom_put("temper25h", zw3d)   ! potential temperature
+         zw3d(:,:,:) = sn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+         CALL iom_put( "salin25h", zw3d  )   ! salinity
+         zw2d(:,:) = sshn_25h(:,:)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1))
+         CALL iom_put( "ssh25h", zw2d )   ! sea surface 
+         ! Write velocities (instantaneous)
+         zw3d(:,:,:) = un_25h(:,:,:)*umask(:,:,:) + zmdi*(1.0-umask(:,:,:))
+         CALL iom_put("vozocrtx25h", zw3d)    ! i-current
+         zw3d(:,:,:) = vn_25h(:,:,:)*vmask(:,:,:) + zmdi*(1.0-vmask(:,:,:))
+         CALL iom_put("vomecrty25h", zw3d  )   ! j-current
+         zw3d(:,:,:) = wn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+         CALL iom_put("vomecrtz25h", zw3d )   ! k-current
+         ! Write vertical physics
+         zw3d(:,:,:) = avt_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+         CALL iom_put("avt25h", zw3d )   ! diffusivity
+         zw3d(:,:,:) = avm_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+         CALL iom_put("avm25h", zw3d)   ! viscosity
+         IF( ln_zdftke ) THEN
             zw3d(:,:,:) = en_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
             CALL iom_put("tke25h", zw3d)   ! tke
-#endif
-#if defined key_zdfgls 
+         ENDIF
+         IF( ln_zdfgls ) THEN
+            zw3d(:,:,:) = en_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+            CALL iom_put("tke25h", zw3d)   ! tke
             zw3d(:,:,:) = rmxln_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
             CALL iom_put( "mxln25h",zw3d)
-#endif
-
-            ! After the write reset the values to cnt=1 and sum values equal current value 
-            tn_25h(:,:,:) = tsn(:,:,:,jp_tem)
-            sn_25h(:,:,:) = tsn(:,:,:,jp_sal)
-            sshn_25h(:,:) = sshn (:,:)
-            un_25h(:,:,:) = un(:,:,:)
-            vn_25h(:,:,:) = vn(:,:,:)
-            wn_25h(:,:,:) = wn(:,:,:)
-            avt_25h(:,:,:) = avt(:,:,:)
-            avm_25h(:,:,:) = avm(:,:,:)
-# if defined key_zdfgls || defined key_zdftke
+         ENDIF
+         !
+         ! After the write reset the values to cnt=1 and sum values equal current value 
+         tn_25h(:,:,:) = tsn(:,:,:,jp_tem)
+         sn_25h(:,:,:) = tsn(:,:,:,jp_sal)
+         sshn_25h(:,:) = sshn (:,:)
+         un_25h(:,:,:) = un(:,:,:)
+         vn_25h(:,:,:) = vn(:,:,:)
+         wn_25h(:,:,:) = wn(:,:,:)
+         avt_25h(:,:,:) = avt(:,:,:)
+         avm_25h(:,:,:) = avm(:,:,:)
+         IF( ln_zdftke ) THEN
             en_25h(:,:,:) = en(:,:,:)
-#endif
-# if defined key_zdfgls
+         ENDIF
+         IF( ln_zdfgls ) THEN
+            en_25h(:,:,:) = en(:,:,:)
             rmxln_25h(:,:,:) = mxln(:,:,:)
-#endif
-            cnt_25h = 1
-            IF (lwp)  WRITE(numout,*) 'dia_wri_tide : After 25hr mean write, reset sum to current value and cnt_25h to one for overlapping average',cnt_25h
-
+         ENDIF
+         cnt_25h = 1
+         IF (lwp)  WRITE(numout,*) 'dia_wri_tide : After 25hr mean write, reset sum to current value and cnt_25h to one for overlapping average',cnt_25h
+         !
       ENDIF !  cnt_25h .EQ. 25 .AND.  MOD( kt, i_steps * 24) == 0 .AND. kt .NE. nn_it000
-
-
+      !
    END SUBROUTINE dia_25h 
 
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 7953)
@@ -39,5 +39,5 @@
    USE zdfmxl          ! mixed layer
    USE dianam          ! build name of file (routine)
-   USE zdfddm          ! vertical  physics: double diffusion
+!   USE zdfddm          ! vertical  physics: double diffusion
    USE diahth          ! thermocline diagnostics
    USE wet_dry         ! wetting and drying
@@ -233,5 +233,5 @@
 
       CALL iom_put( "avt" , avt                        )    ! T vert. eddy diff. coef.
-      CALL iom_put( "avs" , avs                        )    ! S vert. eddy diff. coef. (useful only with key_zdfddm)
+      CALL iom_put( "avs" , avs                        )    ! S vert. eddy diff. coef.
       CALL iom_put( "avm" , avmu                       )    ! T vert. eddy visc. coef.
 
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 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90	(revision 7953)
@@ -96,5 +96,5 @@
       !!                 ***  ROUTINE dyn_zdf_init  ***
       !!
-      !! ** Purpose :   initializations of the vertical diffusion scheme
+      !! ** Purpose :   initialization of the vertical diffusion scheme
       !!
       !! ** Method  :   implicit (euler backward) scheme (default)
@@ -105,5 +105,5 @@
       !!----------------------------------------------------------------------
       !
-      ! Choice from ln_zdfexp read in namelist in zdfini
+      ! Choice from ln_zdfexp (namzdf namelist variable read in zdfphy module)
       IF( ln_zdfexp ) THEN   ;   nzdf = 0           ! use explicit scheme
       ELSE                   ;   nzdf = 1           ! use implicit scheme
@@ -111,5 +111,5 @@
       !
       ! Force implicit schemes
-      IF( lk_zdftke .OR. lk_zdfgls   )   nzdf = 1   ! TKE or GLS physics
+      IF( ln_zdftke .OR. ln_zdfgls   )   nzdf = 1   ! TKE or GLS physics
       IF( ln_dynldf_iso              )   nzdf = 1   ! iso-neutral lateral physics
       IF( ln_dynldf_hor .AND. ln_sco )   nzdf = 1   ! horizontal lateral physics in s-coordinate
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90	(revision 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90	(revision 7953)
@@ -124,5 +124,5 @@
       !!----------------------------------------------------------------------
       !
-      ! Choice from ln_zdfexp already read in namelist in zdfini module
+      ! Choice from ln_zdfexp (namzdf namelist variable read in zdfphy module)
       IF( ln_zdfexp ) THEN   ;   nzdf = 0           ! use explicit scheme
       ELSE                   ;   nzdf = 1           ! use implicit scheme
@@ -130,5 +130,5 @@
       !
       ! Force implicit schemes
-      IF( lk_zdftke .OR. lk_zdfgls   )   nzdf = 1   ! TKE, or GLS physics
+      IF( ln_zdftke .OR. ln_zdfgls   )   nzdf = 1   ! TKE, or GLS physics
       IF( ln_traldf_iso              )   nzdf = 1   ! iso-neutral lateral physics
       IF( ln_traldf_hor .AND. ln_sco )   nzdf = 1   ! horizontal lateral physics in s-coordinate
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 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90	(revision 7953)
@@ -15,10 +15,4 @@
 
    PUBLIC  zdf_oce_alloc    ! Called in nemogcm.F90
-
-#if defined key_zdfcst
-   LOGICAL, PARAMETER, PUBLIC ::   lk_zdfcst        = .TRUE.         !: constant vertical mixing flag
-#else
-   LOGICAL, PARAMETER, PUBLIC ::   lk_zdfcst        = .FALSE.        !: constant vertical mixing flag
-#endif
 
    !                                 !!* namelist namzdf: vertical diffusion *
Index: anches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm_substitute.h90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm_substitute.h90	(revision 7952)
+++ 	(revision )
@@ -1,18 +1,0 @@
-   !!----------------------------------------------------------------------
-   !!                    *** zdfddm_substitute.h90  ***
-   !!----------------------------------------------------------------------
-   !! ** purpose :   substitute fsaht. the eddy diffusivity coeff.
-   !!      with a constant or 1D or 2D or 3D array, using CPP macro.
-   !!----------------------------------------------------------------------
-#if defined key_zdfddm
-!   'key_zdfddm' :                      avs: 3D array defined in zdfddm module
-#   define   fsavs(i,j,k)       avs(i,j,k)
-#else
-!   Defautl option :                     avs = avt
-#   define   fsavs(i,j,k)       avt(i,j,k)
-#endif
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
-   !! $Id$ 
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
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 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 7953)
@@ -8,7 +8,5 @@
    !!             3.3  !  2010-10  (C. Bricaud)  Add in the reference
    !!----------------------------------------------------------------------
-#if defined key_zdfgls
-   !!----------------------------------------------------------------------
-   !!   'key_zdfgls'                 Generic Length Scale vertical physics
+
    !!----------------------------------------------------------------------
    !!   zdf_gls       : update momentum and tracer Kz from a gls scheme
@@ -39,8 +37,7 @@
 
    PUBLIC   zdf_gls        ! routine called in step module
-   PUBLIC   zdf_gls_init   ! routine called in opa module
-   PUBLIC   gls_rst        ! routine called in step module
-
-   LOGICAL , PUBLIC, PARAMETER ::   lk_zdfgls = .TRUE.   !: TKE vertical mixing flag
+   PUBLIC   zdf_gls_init   ! routine called in zdfphy module
+   PUBLIC   gls_rst        ! routine called in zdfphy module
+
    !
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   mxln    !: now mixing length
@@ -1211,23 +1208,4 @@
    END SUBROUTINE gls_rst
 
-#else
-   !!----------------------------------------------------------------------
-   !!   Dummy module :                                        NO TKE scheme
-   !!----------------------------------------------------------------------
-   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfgls = .FALSE.   !: TKE flag
-CONTAINS
-   SUBROUTINE zdf_gls_init           ! Empty routine
-      WRITE(*,*) 'zdf_gls_init: You should not have seen this print! error?'
-   END SUBROUTINE zdf_gls_init
-   SUBROUTINE zdf_gls( kt )          ! Empty routine
-      WRITE(*,*) 'zdf_gls: You should not have seen this print! error?', kt
-   END SUBROUTINE zdf_gls
-   SUBROUTINE gls_rst( kt, cdrw )          ! Empty routine
-      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
-      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
-      WRITE(*,*) 'gls_rst: You should not have seen this print! error?', kt, cdrw
-   END SUBROUTINE gls_rst
-#endif
-
    !!======================================================================
 END MODULE zdfgls
Index: anches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90	(revision 7952)
+++ 	(revision )
@@ -1,207 +1,0 @@
-MODULE zdfini
-   !!======================================================================
-   !!                      ***  MODULE  zdfini  ***
-   !! Ocean physics :   read vertical mixing namelist and check consistancy
-   !!======================================================================
-   !! History :  8.0  ! 1997-06  (G. Madec)  Original code from inimix
-   !!            1.0  ! 2002-08  (G. Madec)  F90 : free form
-   !!             -   ! 2005-06  (C. Ethe) KPP scheme
-   !!             -   ! 2009-07  (G. Madec) add avmb, avtb in restart for cen2 advection
-   !!            3.7  ! 2014-12  (G. Madec) remove KPP scheme
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   zdf_init    : initialization, namelist read, and parameters control
-   !!----------------------------------------------------------------------
-   USE par_oce         ! mesh and scale factors
-   USE zdf_oce         ! TKE vertical mixing          
-   USE sbc_oce         ! surface module (only for nn_isf in the option compatibility test)
-   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 tranpc          ! convection: non penetrative adjustment
-   USE ldfslp          ! iso-neutral slopes
-   !
-   USE in_out_manager  ! I/O manager
-   USE iom             ! IOM library
-   USE lib_mpp         ! distribued memory computing
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   zdf_init   ! routine called by opa.F90
-   
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
-   !! $Id$
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE zdf_init
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE zdf_init  ***
-      !! 
-      !! ** Purpose :   initializations of the vertical ocean physics
-      !!
-      !! ** Method  :   Read namelist namzdf, control logicals 
-      !!----------------------------------------------------------------------
-      INTEGER ::   ioptio, ios       ! local integers
-      !!
-      NAMELIST/namzdf/ ln_zdfcst, ln_zdfric, ln_zdftke, ln_zdfgls,   &     ! type of closure scheme
-         &             ln_zdfevd, nn_evdm, rn_evd ,                  &     ! convection : evd
-         &             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_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
-      !!----------------------------------------------------------------------
-
-      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 )
-
-      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
-         WRITE(numout,*)
-         WRITE(numout,*) 'zdf_init : vertical physics'
-         WRITE(numout,*) '~~~~~~~~'
-         WRITE(numout,*) '   Namelist namzdf : set vertical mixing mixing parameters'
-         WRITE(numout,*) '      vertical closure scheme'
-         WRITE(numout,*) '         constant vertical mixing coefficient    ln_zdfcst = ', ln_zdfcst
-         WRITE(numout,*) '         constant vertical mixing coefficient    ln_zdfric = ', ln_zdfric
-         WRITE(numout,*) '         constant vertical mixing coefficient    ln_zdftke = ', ln_zdftke
-         WRITE(numout,*) '         constant vertical mixing coefficient    ln_zdfgls = ', ln_zdfgls
-         WRITE(numout,*) '      convection: '
-         WRITE(numout,*) '         enhanced vertical diffusion             ln_zdfevd = ', ln_zdfevd
-         WRITE(numout,*) '            applied on momentum (=1/0)             nn_evdm = ', nn_evdm
-         WRITE(numout,*) '            vertical coefficient for evd           rn_evd  = ', rn_evd
-         WRITE(numout,*) '         non-penetrative convection (npc)        ln_zdfnpc = ', ln_zdfnpc
-         WRITE(numout,*) '            npc call  frequency                    nn_npc  = ', nn_npc
-         WRITE(numout,*) '            npc print frequency                    nn_npcp = ', nn_npcp
-         WRITE(numout,*) '      double diffusive mixing                    ln_zdfddm = ', ln_zdfddm
-         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,*) '         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
-
-      IF(ln_zdfddm) THEN                    ! double diffusive mixing'
-         ALLOCATE( avs(jpi,jpj,jpk) )
-         avs(:,:,:) = rn_avt0 * wmask(:,:,:) 
-      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( lk_zdfcst ) THEN
-         IF(lwp) WRITE(numout,*) '      constant eddy diffusion coefficients'
-         ioptio = ioptio+1
-      ENDIF
-      IF( lk_zdfric ) THEN
-         IF(lwp) WRITE(numout,*) '      Richardson dependent eddy coefficients'
-         ioptio = ioptio+1
-      ENDIF
-      IF( lk_zdftke ) THEN
-         IF(lwp) WRITE(numout,*) '      TKE dependent eddy coefficients'
-         ioptio = ioptio+1
-      ENDIF
-      IF( lk_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( ( lk_zdfric .OR. lk_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_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( lk_zdftke ) THEN
-         IF(lwp) WRITE(numout,*) '      use the 1.5 turbulent closure'
-      ENDIF
-      IF( lk_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.( lk_zdftke .OR. lk_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
-         avmb(:) = rn_avm0
-         avtb(:) = rn_avt0                     
-      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
-           WHERE(-15. <= gphit .AND. gphit < -5 )   avtb_2d = (1.  - 0.09 * (gphit + 15.))
-           WHERE( -5. <= gphit .AND. gphit <  5 )   avtb_2d =  0.1
-           WHERE(  5. <= gphit .AND. gphit < 15 )   avtb_2d = (0.1 + 0.09 * (gphit -  5.))
-      ENDIF
-      !
-   END SUBROUTINE zdf_init
-
-   !!======================================================================
-END MODULE zdfini
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 7953)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfphy.F90	(revision 7953)
@@ -0,0 +1,282 @@
+MODULE zdfphy
+   !!======================================================================
+   !!                      ***  MODULE  zdfphy  ***
+   !! Ocean physics :   manager of vertical mixing parametrizations
+   !!======================================================================
+   !! History :  4.0  !  2017-04  (G. Madec)  original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   zdf_phy_init  : initialization of all vertical physics pakages
+   !!   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 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 lib_mpp        ! distribued memory computing
+
+   IMPLICIT NONE
+   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)
+   !! $Id$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE zdf_phy_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE zdf_phy_init  ***
+      !! 
+      !! ** Purpose :   initializations of the vertical ocean physics
+      !!
+      !! ** Method  :   Read namelist namzdf, control logicals 
+      !!----------------------------------------------------------------------
+      INTEGER ::   ioptio, ios   ! local integers
+      !!
+      NAMELIST/namzdf/ ln_zdfcst, ln_zdfric, ln_zdftke, ln_zdfgls,   &     ! type of closure scheme
+         &             ln_zdfevd, nn_evdm, rn_evd ,                  &     ! convection : evd
+         &             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_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
+      !!----------------------------------------------------------------------
+
+      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 )
+
+      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
+         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'
+         WRITE(numout,*) '         constant vertical mixing coefficient    ln_zdfcst = ', ln_zdfcst
+         WRITE(numout,*) '         constant vertical mixing coefficient    ln_zdfric = ', ln_zdfric
+         WRITE(numout,*) '         constant vertical mixing coefficient    ln_zdftke = ', ln_zdftke
+         WRITE(numout,*) '         constant vertical mixing coefficient    ln_zdfgls = ', ln_zdfgls
+         WRITE(numout,*) '      convection: '
+         WRITE(numout,*) '         enhanced vertical diffusion             ln_zdfevd = ', ln_zdfevd
+         WRITE(numout,*) '            applied on momentum (=1/0)             nn_evdm = ', nn_evdm
+         WRITE(numout,*) '            vertical coefficient for evd           rn_evd  = ', rn_evd
+         WRITE(numout,*) '         non-penetrative convection (npc)        ln_zdfnpc = ', ln_zdfnpc
+         WRITE(numout,*) '            npc call  frequency                    nn_npc  = ', nn_npc
+         WRITE(numout,*) '            npc print frequency                    nn_npcp = ', nn_npcp
+         WRITE(numout,*) '      double diffusive mixing                    ln_zdfddm = ', ln_zdfddm
+         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,*) '         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
+         avmb(:) = rn_avm0
+         avtb(:) = rn_avt0                     
+      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
+           WHERE(-15. <= gphit .AND. gphit < -5 )   avtb_2d = (1.  - 0.09 * (gphit + 15.))
+           WHERE( -5. <= gphit .AND. gphit <  5 )   avtb_2d =  0.1
+           WHERE(  5. <= gphit .AND. gphit < 15 )   avtb_2d = (0.1 + 0.09 * (gphit -  5.))
+      ENDIF
+      !
+
+!!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
+      !
+   END SUBROUTINE zdf_phy_init
+
+
+   SUBROUTINE zdf_phy( kstp )
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE zdf_phy  ***
+      !!
+      !! ** Purpose :  Update ocean physics at each time-step
+      !!
+      !! ** Method  : 
+      !!
+      !! ** Action  :   avm, avt vertical eddy viscosity and diffusivity at w-points
+      !!                nmld ??? mixed layer depth in level and meters   <<<<====verifier !
+      !!                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
+
+                                                      ! 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' )
+      !
+   END SUBROUTINE zdf_phy
+
+   !!======================================================================
+END MODULE zdfphy
Index: /branches/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 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfqiao.F90	(revision 7953)
@@ -67,4 +67,6 @@
       !---------------------------------------------------------------------------------
       !
+!!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
@@ -101,4 +103,6 @@
       !-------------------------------
       !
+!!gm  with double diffusion activated, avs is not updated...
+!!gm                    =====>>> BUG
       DO jk = 1, jpkm1
          DO jj = 1, jpj
@@ -112,4 +116,5 @@
       !
    END SUBROUTINE zdf_qiao
+
 
    INTEGER FUNCTION zdf_qiao_alloc()
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 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90	(revision 7953)
@@ -13,10 +13,7 @@
    !!            3.3.1! 2011-09  (P. Oddo) Mixed layer depth parameterization
    !!----------------------------------------------------------------------
-#if defined key_zdfric
-   !!----------------------------------------------------------------------
-   !!   'key_zdfric'                                             Kz = f(Ri)
-   !!----------------------------------------------------------------------
-   !!   zdf_ric       : update momentum and tracer Kz from the Richardson
-   !!                  number computation
+
+   !!----------------------------------------------------------------------
+   !!   zdf_ric       : update momentum and tracer Kz from the Richardson number
    !!   zdf_ric_init  : initialization, namelist read, & parameters control
    !!----------------------------------------------------------------------
@@ -38,6 +35,4 @@
    PUBLIC   zdf_ric         ! called by step.F90
    PUBLIC   zdf_ric_init    ! called by opa.F90
-
-   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfric = .TRUE.   !: Richardson vertical mixing flag
 
    !                           !!* Namelist namzdf_ric : Richardson number dependent Kz *
@@ -108,5 +103,5 @@
       !!      namelist
       !!        N.B. the mask are required for implicit scheme, and surface
-      !!      and bottom value already set in zdfini.F90
+      !!      and bottom value already set in zdfphy.F90
       !!
       !! References : Pacanowski & Philander 1981, JPO, 1441-1451.
@@ -183,5 +178,5 @@
             zrhos          = rhop(ji,jj,1) + zflageos * ( 1. - tmask(ji,jj,1) )
             zustar         = SQRT( taum(ji,jj) / ( zrhos +  rsmall ) )
-            ekm_dep(ji,jj) = rn_ekmfc * zustar / ( ABS( ff(ji,jj) ) + rsmall )
+            ekm_dep(ji,jj) = 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
@@ -303,17 +298,4 @@
    END SUBROUTINE zdf_ric_init
 
-#else
-   !!----------------------------------------------------------------------
-   !!   Dummy module :              NO Richardson dependent vertical mixing
-   !!----------------------------------------------------------------------
-   LOGICAL, PUBLIC, PARAMETER ::   lk_zdfric = .FALSE.   !: Richardson mixing flag
-CONTAINS
-   SUBROUTINE zdf_ric_init         ! Dummy routine
-   END SUBROUTINE zdf_ric_init
-   SUBROUTINE zdf_ric( kt )        ! Dummy routine
-      WRITE(*,*) 'zdf_ric: You should not have seen this print! error?', kt
-   END SUBROUTINE zdf_ric
-#endif
-
    !!======================================================================
 END MODULE zdfric
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 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 7953)
@@ -27,8 +27,7 @@
    !!            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
    !!----------------------------------------------------------------------
-#if defined key_zdftke
-   !!----------------------------------------------------------------------
-   !!   'key_zdftke'                                   TKE vertical physics
+
    !!----------------------------------------------------------------------
    !!   zdf_tke       : update momentum and tracer Kz from a tke scheme
@@ -65,6 +64,4 @@
    PUBLIC   tke_rst        ! routine called in step module
 
-   LOGICAL , PUBLIC, PARAMETER ::   lk_zdftke = .TRUE.  !: TKE vertical mixing flag
-
    !                      !!** Namelist  namzdf_tke  **
    LOGICAL  ::   ln_mxl0   ! mixing length scale surface value as function of wind stress or not
@@ -376,16 +373,11 @@
             DO ji = fs_2, fs_jpim1   ! vector opt.
                zcof   = zfact1 * tmask(ji,jj,jk)
-# if defined key_zdftmx_new
-               ! key_zdftmx_new: New internal wave-driven param: set a minimum value for Kz on TKE (ensure numerical stability)
-               zzd_up = zcof * MAX( avm(ji,jj,jk+1) + avm(ji,jj,jk), 2.e-5_wp )   &  ! upper diagonal
-                  &          / (  e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk  )  )
-               zzd_lw = zcof * MAX( avm(ji,jj,jk) + avm(ji,jj,jk-1), 2.e-5_wp )   &  ! lower diagonal
-                  &          / (  e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk  )  )
-# else
-               zzd_up = zcof * ( avm  (ji,jj,jk+1) + avm  (ji,jj,jk  ) )   &  ! upper diagonal
-                  &          / ( e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk  ) )
-               zzd_lw = zcof * ( avm  (ji,jj,jk  ) + avm  (ji,jj,jk-1) )   &  ! lower diagonal
-                  &          / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk  ) )
-# endif
+               !                                   ! A minimum of 2.e-5 m2/s is imposed on TKE vertical
+               !                                   ! eddy coefficient (ensure numerical stability)
+               zzd_up = zcof * MAX(   avm(ji,jj,jk+1) +   avm(ji,jj,jk  ) , 2.e-5_wp  )   &  ! upper diagonal
+                  &          /    ( e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk  )  )
+               zzd_lw = zcof * MAX(   avm(ji,jj,jk  ) +   avm(ji,jj,jk-1) , 2.e-5_wp  )   &  ! lower diagonal
+                  &          /    ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk  )  )
+               !
                !                                   ! shear prod. at w-point weightened by mask
                zesh2  =  ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
@@ -741,16 +733,4 @@
       !
       ri_cri   = 2._wp    / ( 2._wp + rn_ediss / rn_ediff )   ! resulting critical Richardson number
-# if defined key_zdftmx_new
-      ! key_zdftmx_new: New internal wave-driven param: specified value of rn_emin & rmxl_min are used
-      rn_emin  = 1.e-10_wp
-      rmxl_min = 1.e-03_wp
-      IF(lwp) THEN                  ! Control print
-         WRITE(numout,*)
-         WRITE(numout,*) 'zdf_tke_init :  New tidal mixing case: force rn_emin = 1.e-10 and rmxl_min = 1.e-3 '
-         WRITE(numout,*) '~~~~~~~~~~~~'
-      ENDIF
-# else
-      rmxl_min = 1.e-6_wp / ( rn_ediff * SQRT( rn_emin ) )    ! resulting minimum length to recover molecular viscosity
-# endif
       !
       IF(lwp) THEN                    !* Control print
@@ -776,4 +756,15 @@
          WRITE(numout,*)
          WRITE(numout,*) '      critical Richardson nb with your parameters  ri_cri = ', ri_cri
+         WRITE(numout,*)
+      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(lwp) WRITE(numout,*) '      Internal wave-driven mixing case:   force   rn_emin = 1.e-10 and rmxl_min = 1.e-3 '
+      ELSE
+         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
       ENDIF
       !
@@ -891,21 +882,4 @@
    END SUBROUTINE tke_rst
 
-#else
-   !!----------------------------------------------------------------------
-   !!   Dummy module :                                        NO TKE scheme
-   !!----------------------------------------------------------------------
-   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftke = .FALSE.   !: TKE flag
-CONTAINS
-   SUBROUTINE zdf_tke_init           ! Dummy routine
-   END SUBROUTINE zdf_tke_init
-   SUBROUTINE zdf_tke( kt )          ! Dummy routine
-      WRITE(*,*) 'zdf_tke: You should not have seen this print! error?', kt
-   END SUBROUTINE zdf_tke
-   SUBROUTINE tke_rst( kt, cdrw )
-     CHARACTER(len=*) ::   cdrw
-     WRITE(*,*) 'tke_rst: You should not have seen this print! error?', kt, cdwr
-   END SUBROUTINE tke_rst
-#endif
-
    !!======================================================================
 END MODULE zdftke
Index: /branches/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 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90	(revision 7953)
@@ -2,546 +2,13 @@
    !!========================================================================
    !!                       ***  MODULE  zdftmx  ***
-   !! Ocean physics: vertical tidal mixing coefficient
+   !! 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
+   !!             -   !  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)
    !!----------------------------------------------------------------------
-#if defined key_zdftmx
-   !!----------------------------------------------------------------------
-   !!   'key_zdftmx'                                  Tidal vertical mixing
-   !!----------------------------------------------------------------------
-   !!   zdf_tmx       : global     momentum & tracer Kz with tidal induced Kz
-   !!   tmx_itf       : Indonesian momentum & tracer Kz with tidal 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 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 opa module 
-   PUBLIC   zdf_tmx_alloc   ! called in nemogcm module
-
-   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftmx = .TRUE.    !: tidal mixing flag
-
-   !                       !!* Namelist  namzdf_tmx : tidal mixing *
-   REAL(wp) ::  rn_htmx     ! vertical decay scale for turbulence (meters)
-   REAL(wp) ::  rn_n2min    ! threshold of the Brunt-Vaisala frequency (s-1)
-   REAL(wp) ::  rn_tfe      ! tidal dissipation efficiency (St Laurent et al. 2002)
-   REAL(wp) ::  rn_me       ! mixing efficiency (Osborn 1980)
-   LOGICAL  ::  ln_tmx_itf  ! Indonesian Through Flow (ITF): Koch-Larrouy et al. (2007) parameterization
-   REAL(wp) ::  rn_tfe_itf  ! ITF tidal dissipation efficiency (St Laurent et al. 2002)
-
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   en_tmx     ! energy available for tidal mixing (W/m2)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   mask_itf   ! mask to use over Indonesian area
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   az_tmx     ! coefficient used to evaluate the tidal induced Kz
-
-   !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
-   !! $Id$
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   INTEGER FUNCTION zdf_tmx_alloc()
-      !!----------------------------------------------------------------------
-      !!                ***  FUNCTION zdf_tmx_alloc  ***
-      !!----------------------------------------------------------------------
-      ALLOCATE(en_tmx(jpi,jpj), mask_itf(jpi,jpj), az_tmx(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
-      !!              tidal mixing (Simmons et al 2004).
-      !!
-      !! ** Method  : - tidal-induced vertical mixing is given by:
-      !!                  Kz_tides = az_tmx / max( rn_n2min, N^2 )
-      !!              where az_tmx is a coefficient that specified the 3D space 
-      !!              distribution of the faction of tidal energy taht is used
-      !!              for mixing. Its expression is set in zdf_tmx_init routine,
-      !!              following Simmons et al. 2004.
-      !!                NB: a specific bounding procedure is performed on av_tide
-      !!              so that the input tidal energy is actually almost used. The
-      !!              basic maximum value is 60 cm2/s, but values of 300 cm2/s 
-      !!              can be reached in area where bottom stratification is too 
-      !!              weak.
-      !!
-      !!              - update av_tide in the Indonesian Through Flow area
-      !!              following Koch-Larrouy et al. (2007) parameterisation
-      !!              (see tmx_itf routine).
-      !!
-      !!              - update the model vertical eddy viscosity and diffusivity: 
-      !!                     avt  = avt  +    av_tides
-      !!                     avm  = avm  +    av_tides
-      !!                     avmu = avmu + mi(av_tides)
-      !!                     avmv = avmv + mj(av_tides)
-      !!
-      !! ** Action  :   avt, avm, avmu, avmv   increased by tidal mixing
-      !!
-      !! References : Simmons et al. 2004, Ocean Modelling, 6, 3-4, 245-263.
-      !!              Koch-Larrouy et al. 2007, GRL.
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step 
-      !
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   ztpc         ! scalar workspace
-      REAL(wp), POINTER, DIMENSION(:,:)   ::   zkz
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zav_tide
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('zdf_tmx')
-      !
-      CALL wrk_alloc( jpi,jpj,       zkz )
-      CALL wrk_alloc( jpi,jpj,jpk,   zav_tide )
-      !
-      !                          ! ----------------------- !
-      !                          !  Standard tidal mixing  !  (compute zav_tide)
-      !                          ! ----------------------- !
-      !                             !* First estimation (with n2 bound by rn_n2min) bounded by 60 cm2/s
-      zav_tide(:,:,:) = MIN(  60.e-4, az_tmx(:,:,:) / MAX( rn_n2min, rn2(:,:,:) )  )
-
-      zkz(:,:) = 0.e0               !* Associated potential energy consummed over the whole water column
-      DO jk = 2, jpkm1
-         zkz(:,:) = zkz(:,:) + e3w_n(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zav_tide(:,:,jk) * wmask(:,:,jk)
-      END DO
-
-      DO jj = 1, jpj                !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx
-         DO ji = 1, jpi
-            IF( zkz(ji,jj) /= 0.e0 )   zkz(ji,jj) = en_tmx(ji,jj) / zkz(ji,jj)
-         END DO
-      END DO
-
-      DO jk = 2, jpkm1     !* Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zav_tide bound by 300 cm2/s
-         zav_tide(:,:,jk) = zav_tide(:,:,jk) * MIN( zkz(:,:), 30./6. ) * wmask(:,:,jk)  !kz max = 300 cm2/s
-      END DO
-
-      IF( kt == nit000 ) THEN       !* check at first time-step: diagnose the energy consumed by zav_tide
-         ztpc = 0._wp
-         DO jk= 1, jpk
-            DO jj= 1, jpj
-               DO ji= 1, jpi
-                  ztpc = ztpc + e3w_n(ji,jj,jk) * e1e2t(ji,jj)                  &
-                     &        * MAX( 0.e0, rn2(ji,jj,jk) ) * zav_tide(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
-               END DO
-            END DO
-         END DO
-         ztpc= rau0 / ( rn_tfe * rn_me ) * ztpc
-         IF( lk_mpp )   CALL mpp_sum( ztpc )
-         IF(lwp) WRITE(numout,*) 
-         IF(lwp) WRITE(numout,*) '          N Total power consumption by av_tide    : ztpc = ', ztpc * 1.e-12 ,'TW'
-      ENDIF
-       
-      !                          ! ----------------------- !
-      !                          !    ITF  tidal mixing    !  (update zav_tide)
-      !                          ! ----------------------- !
-      IF( ln_tmx_itf )   CALL tmx_itf( kt, zav_tide )
-
-      !                          ! ----------------------- !
-      !                          !   Update  mixing coefs  !                          
-      !                          ! ----------------------- !
-      DO jk = 2, jpkm1              !* update momentum & tracer diffusivity with tidal mixing
-         avt(:,:,jk) = avs(:,:,jk) + zav_tide(:,:,jk) * wmask(:,:,jk)
-         avm(:,:,jk) = avm(:,:,jk) + zav_tide(:,:,jk) * wmask(:,:,jk)
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1  ! vector opt.
-               avmu(ji,jj,jk) = avmu(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji+1,jj  ,jk) ) * wumask(ji,jj,jk)
-               avmv(ji,jj,jk) = avmv(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(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 tidal mixing coefficient
-      CALL iom_put( "av_tide", zav_tide )
-
-      IF(ln_ctl)   CALL prt_ctl(tab3d_1=zav_tide , clinfo1=' tmx - av_tide: ', tab3d_2=avt, clinfo2=' avt: ', ovlap=1, kdim=jpk)
-      !
-      CALL wrk_dealloc( jpi,jpj,       zkz )
-      CALL wrk_dealloc( jpi,jpj,jpk,   zav_tide )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('zdf_tmx')
-      !
-   END SUBROUTINE zdf_tmx
-
-
-   SUBROUTINE tmx_itf( kt, pav )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE tmx_itf  ***
-      !!                   
-      !! ** Purpose :   modify the vertical eddy diffusivity coefficients 
-      !!              (pav) in the Indonesian Through Flow area (ITF).
-      !!
-      !! ** Method  : - Following Koch-Larrouy et al. (2007), in the ITF defined
-      !!                by msk_itf (read in a file, see tmx_init), the tidal
-      !!                mixing coefficient is computed with :
-      !!                  * q=1 (i.e. all the tidal energy remains trapped in
-      !!                         the area and thus is used for mixing)
-      !!                  * the vertical distribution of the tifal energy is a
-      !!                    proportional to N above the thermocline (d(N^2)/dz > 0)
-      !!                    and to N^2 below the thermocline (d(N^2)/dz < 0)
-      !!
-      !! ** Action  :   av_tide   updated in the ITF area (msk_itf)
-      !!
-      !! References :  Koch-Larrouy et al. 2007, GRL 
-      !!----------------------------------------------------------------------
-      INTEGER , INTENT(in   )                         ::   kt   ! ocean time-step
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   pav  ! Tidal mixing coef.
-      !! 
-      INTEGER  ::   ji, jj, jk    ! dummy loop indices
-      REAL(wp) ::   zcoef, ztpc   ! temporary scalar
-      REAL(wp), DIMENSION(:,:)  , POINTER ::   zkz                        ! 2D workspace
-      REAL(wp), DIMENSION(:,:)  , POINTER ::   zsum1 , zsum2 , zsum       !  -      -
-      REAL(wp), DIMENSION(:,:,:), POINTER ::   zempba_3d_1, zempba_3d_2   ! 3D workspace
-      REAL(wp), DIMENSION(:,:,:), POINTER ::   zempba_3d  , zdn2dz        !  -      -
-      REAL(wp), DIMENSION(:,:,:), POINTER ::   zavt_itf                   !  -      -
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('tmx_itf')
-      !
-      CALL wrk_alloc( jpi,jpj, zkz, zsum1 , zsum2 , zsum )
-      CALL wrk_alloc( jpi,jpj,jpk, zempba_3d_1, zempba_3d_2, zempba_3d, zdn2dz, zavt_itf )
-
-      !                             ! compute the form function using N2 at each time step
-      zempba_3d_1(:,:,jpk) = 0.e0
-      zempba_3d_2(:,:,jpk) = 0.e0
-      DO jk = 1, jpkm1             
-         zdn2dz     (:,:,jk) = rn2(:,:,jk) - rn2(:,:,jk+1)           ! Vertical profile of dN2/dz
-         zempba_3d_1(:,:,jk) = SQRT(  MAX( 0.e0, rn2(:,:,jk) )  )    !    -        -    of N
-         zempba_3d_2(:,:,jk) =        MAX( 0.e0, rn2(:,:,jk) )       !    -        -    of N^2
-      END DO
-      !
-      zsum (:,:) = 0.e0
-      zsum1(:,:) = 0.e0
-      zsum2(:,:) = 0.e0
-      DO jk= 2, jpk
-         zsum1(:,:) = zsum1(:,:) + zempba_3d_1(:,:,jk) * e3w_n(:,:,jk) * wmask(:,:,jk)
-         zsum2(:,:) = zsum2(:,:) + zempba_3d_2(:,:,jk) * e3w_n(:,:,jk) * wmask(:,:,jk)               
-      END DO
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            IF( zsum1(ji,jj) /= 0.e0 )   zsum1(ji,jj) = 1.e0 / zsum1(ji,jj)
-            IF( zsum2(ji,jj) /= 0.e0 )   zsum2(ji,jj) = 1.e0 / zsum2(ji,jj)                
-         END DO
-      END DO
-
-      DO jk= 1, jpk
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zcoef = 0.5 - SIGN( 0.5, zdn2dz(ji,jj,jk) )       ! =0 if dN2/dz > 0, =1 otherwise 
-               ztpc  = zempba_3d_1(ji,jj,jk) * zsum1(ji,jj) *        zcoef     &
-                  &  + zempba_3d_2(ji,jj,jk) * zsum2(ji,jj) * ( 1. - zcoef )
-               !
-               zempba_3d(ji,jj,jk) =               ztpc 
-               zsum     (ji,jj)    = zsum(ji,jj) + ztpc * e3w_n(ji,jj,jk)
-            END DO
-         END DO
-       END DO
-       DO jj = 1, jpj
-          DO ji = 1, jpi
-             IF( zsum(ji,jj) > 0.e0 )   zsum(ji,jj) = 1.e0 / zsum(ji,jj)                
-          END DO
-       END DO
-
-      !                             ! first estimation bounded by 10 cm2/s (with n2 bounded by rn_n2min) 
-      zcoef = rn_tfe_itf / ( rn_tfe * rau0 )
-      DO jk = 1, jpk
-         zavt_itf(:,:,jk) = MIN(  10.e-4, zcoef * en_tmx(:,:) * zsum(:,:) * zempba_3d(:,:,jk)   &
-            &                                      / MAX( rn_n2min, rn2(:,:,jk) ) * tmask(:,:,jk)  )
-      END DO           
-
-      zkz(:,:) = 0.e0               ! Associated potential energy consummed over the whole water column
-      DO jk = 2, jpkm1
-         zkz(:,:) = zkz(:,:) + e3w_n(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zavt_itf(:,:,jk) * wmask(:,:,jk)
-      END DO
-
-      DO jj = 1, jpj                ! Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx
-         DO ji = 1, jpi
-            IF( zkz(ji,jj) /= 0.e0 )   zkz(ji,jj) = en_tmx(ji,jj) * rn_tfe_itf / rn_tfe / zkz(ji,jj)
-         END DO
-      END DO
-
-      DO jk = 2, jpkm1              ! Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zavt_itf bound by 300 cm2/s
-         zavt_itf(:,:,jk) = zavt_itf(:,:,jk) * MIN( zkz(:,:), 120./10. ) * wmask(:,:,jk)   ! kz max = 120 cm2/s
-      END DO
-
-      IF( kt == nit000 ) THEN       ! diagnose the nergy consumed by zavt_itf
-         ztpc = 0.e0
-         DO jk= 1, jpk
-            DO jj= 1, jpj
-               DO ji= 1, jpi
-                  ztpc = ztpc + e1e2t(ji,jj) * e3w_n(ji,jj,jk) * MAX( 0.e0, rn2(ji,jj,jk) )   &
-                     &                       * zavt_itf(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
-               END DO
-            END DO
-         END DO
-         IF( lk_mpp )   CALL mpp_sum( ztpc )
-         ztpc= rau0 * ztpc / ( rn_me * rn_tfe_itf )
-         IF(lwp) WRITE(numout,*) '          N Total power consumption by zavt_itf: ztpc = ', ztpc * 1.e-12 ,'TW'
-      ENDIF
-
-      !                             ! Update pav with the ITF mixing coefficient
-      DO jk = 2, jpkm1
-         pav(:,:,jk) = pav     (:,:,jk) * ( 1.e0 - mask_itf(:,:) )   &
-            &        + zavt_itf(:,:,jk) *          mask_itf(:,:) 
-      END DO
-      !
-      CALL wrk_dealloc( jpi,jpj, zkz, zsum1 , zsum2 , zsum )
-      CALL wrk_dealloc( jpi,jpj,jpk, zempba_3d_1, zempba_3d_2, zempba_3d, zdn2dz, zavt_itf )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('tmx_itf')
-      !
-   END SUBROUTINE tmx_itf
-
-
-   SUBROUTINE zdf_tmx_init
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE zdf_tmx_init  ***
-      !!                     
-      !! ** Purpose :   Initialization of the vertical tidal mixing, Reading
-      !!              of M2 and K1 tidal energy in nc files
-      !!
-      !! ** Method  : - Read the namtmx namelist and check the parameters
-      !!
-      !!              - Read the input data in NetCDF files :
-      !!              M2 and K1 tidal energy. The total tidal energy, en_tmx, 
-      !!              is the sum of M2, K1 and S2 energy where S2 is assumed 
-      !!              to be: S2=(1/2)^2 * M2
-      !!              mask_itf, a mask array that determine where substituing 
-      !!              the standard Simmons et al. (2005) formulation with the
-      !!              one of Koch_Larrouy et al. (2007).
-      !!
-      !!              - Compute az_tmx, a 3D coefficient that allows to compute
-      !!             the standard tidal-induced vertical mixing as follows:
-      !!                  Kz_tides = az_tmx / max( rn_n2min, N^2 )
-      !!             with az_tmx a bottom intensified coefficient is given by:
-      !!                 az_tmx(z) = en_tmx / ( rau0 * rn_htmx ) * EXP( -(H-z)/rn_htmx )
-      !!                                                  / ( 1. - EXP( - H   /rn_htmx ) ) 
-      !!             where rn_htmx the characteristic length scale of the bottom 
-      !!             intensification, en_tmx the tidal energy, and H the ocean depth
-      !!
-      !! ** input   :   - Namlist namtmx
-      !!                - NetCDF file : M2_ORCA2.nc, K1_ORCA2.nc, and mask_itf.nc
-      !!
-      !! ** Action  : - Increase by 1 the nstop flag is setting problem encounter
-      !!              - defined az_tmx used to compute tidal-induced mixing
-      !!
-      !! References : Simmons et al. 2004, Ocean Modelling, 6, 3-4, 245-263.
-      !!              Koch-Larrouy et al. 2007, GRL.
-      !!----------------------------------------------------------------------
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      INTEGER  ::   inum         ! local integer
-      INTEGER  ::   ios
-      REAL(wp) ::   ztpc, ze_z   ! local scalars
-      REAL(wp), DIMENSION(:,:)  , POINTER ::  zem2, zek1     ! read M2 and K1 tidal energy
-      REAL(wp), DIMENSION(:,:)  , POINTER ::  zkz            ! total M2, K1 and S2 tidal energy
-      REAL(wp), DIMENSION(:,:)  , POINTER ::  zfact          ! used for vertical structure function
-      REAL(wp), DIMENSION(:,:)  , POINTER ::  zhdep          ! Ocean depth 
-      REAL(wp), DIMENSION(:,:,:), POINTER ::  zpc, zav_tide  ! power consumption
-      !!
-      NAMELIST/namzdf_tmx/ rn_htmx, rn_n2min, rn_tfe, rn_me, ln_tmx_itf, rn_tfe_itf
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('zdf_tmx_init')
-      !
-      CALL wrk_alloc( jpi,jpj,       zem2, zek1, zkz, zfact, zhdep )
-      CALL wrk_alloc( jpi,jpj,jpk,   zpc, zav_tide )
-      !
-      REWIND( numnam_ref )             ! Namelist namzdf_tmx in reference namelist : Tidal Mixing
-      READ  ( numnam_ref, namzdf_tmx, 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 : Tidal Mixing
-      READ  ( numnam_cfg, namzdf_tmx, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tmx in configuration namelist', lwp )
-      IF(lwm) WRITE ( numond, namzdf_tmx )
-      !
-      IF(lwp) THEN                     ! Control print
-         WRITE(numout,*)
-         WRITE(numout,*) 'zdf_tmx_init : tidal mixing'
-         WRITE(numout,*) '~~~~~~~~~~~~'
-         WRITE(numout,*) '   Namelist namzdf_tmx : set tidal mixing parameters'
-         WRITE(numout,*) '      Vertical decay scale for turbulence   = ', rn_htmx 
-         WRITE(numout,*) '      Brunt-Vaisala frequency threshold     = ', rn_n2min
-         WRITE(numout,*) '      Tidal dissipation efficiency          = ', rn_tfe
-         WRITE(numout,*) '      Mixing efficiency                     = ', rn_me
-         WRITE(numout,*) '      ITF specific parameterisation         = ', ln_tmx_itf
-         WRITE(numout,*) '      ITF tidal dissipation efficiency      = ', rn_tfe_itf
-      ENDIF
-      !                                ! allocate tmx arrays
-      IF( zdf_tmx_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_tmx_init : unable to allocate tmx arrays' )
-
-      IF( ln_tmx_itf ) THEN            ! read the Indonesian Through Flow mask
-         CALL iom_open('mask_itf',inum)
-         CALL iom_get (inum, jpdom_data, 'tmaskitf',mask_itf,1) ! 
-         CALL iom_close(inum)
-      ENDIF
-      !                                ! read M2 tidal energy flux : W/m2  ( zem2 < 0 )
-      CALL iom_open('M2rowdrg',inum)
-      CALL iom_get (inum, jpdom_data, 'field',zem2,1) ! 
-      CALL iom_close(inum)
-      !                                ! read K1 tidal energy flux : W/m2  ( zek1 < 0 )
-      CALL iom_open('K1rowdrg',inum)
-      CALL iom_get (inum, jpdom_data, 'field',zek1,1) ! 
-      CALL iom_close(inum)
-      !                                ! Total tidal energy ( M2, S2 and K1  with S2=(1/2)^2 * M2 )
-      !                                ! only the energy available for mixing is taken into account,
-      !                                ! (mixing efficiency tidal dissipation efficiency)
-      en_tmx(:,:) = - rn_tfe * rn_me * ( zem2(:,:) * 1.25 + zek1(:,:) ) * ssmask(:,:)
-
-!============
-!TG: Bug for VVL? Should this section be moved out of _init and be updated at every timestep?
-!!gm : you are right, but tidal mixing acts in deep ocean (H>500m) where e3 is O(100m)
-!!     the error is thus ~1% which I feel comfortable with, compared to uncertainties in tidal energy dissipation.
-      !                                ! Vertical structure (az_tmx)
-      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 * rn_htmx * ( 1. - EXP( -zhdep(ji,jj) / rn_htmx ) )
-            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = en_tmx(ji,jj) / zfact(ji,jj)
-         END DO
-      END DO
-      DO jk= 1, jpk                          ! complete with the level-dependent part
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               az_tmx(ji,jj,jk) = zfact(ji,jj) * EXP( -( zhdep(ji,jj)-gdepw_0(ji,jj,jk) ) / rn_htmx ) * tmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-!===========
-      !
-      IF( nprint == 1 .AND. lwp ) THEN
-         ! Control print
-         ! Total power consumption due to vertical mixing
-         ! zpc = rau0 * 1/rn_me * rn2 * zav_tide
-         zav_tide(:,:,:) = 0.e0
-         DO jk = 2, jpkm1
-            zav_tide(:,:,jk) = az_tmx(:,:,jk) / MAX( rn_n2min, rn2(:,:,jk) )
-         END DO
-         !
-         ztpc = 0._wp
-         zpc(:,:,:) = MAX(rn_n2min,rn2(:,:,:)) * zav_tide(:,:,:)
-         DO jk= 2, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ztpc = ztpc + e3w_n(ji,jj,jk) * e1e2t(ji,jj) * zpc(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 * 1/(rn_tfe * rn_me) * ztpc
-         !
-         WRITE(numout,*) 
-         WRITE(numout,*) '          Total power consumption of the tidally driven part of Kz : ztpc = ', ztpc * 1.e-12 ,'TW'
-         !
-         ! control print 2
-         zav_tide(:,:,:) = MIN( zav_tide(:,:,:), 60.e-4 )   
-         zkz(:,:) = 0._wp
-         DO jk = 2, jpkm1
-               zkz(:,:) = zkz(:,:) + e3w_n(:,:,jk) * MAX(0.e0, rn2(:,:,jk)) * rau0 * zav_tide(:,:,jk) * wmask(:,:,jk)
-         END DO
-         ! Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF( zkz(ji,jj) /= 0.e0 )   THEN
-                   zkz(ji,jj) = en_tmx(ji,jj) / zkz(ji,jj)
-               ENDIF
-            END DO
-         END DO
-         ztpc = 1.e50
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF( zkz(ji,jj) /= 0.e0 )   THEN
-                   ztpc = Min( zkz(ji,jj), ztpc)
-               ENDIF
-            END DO
-         END DO
-         WRITE(numout,*) '          Min de zkz ', ztpc, ' Max = ', maxval(zkz(:,:) )
-         !
-         DO jk = 2, jpkm1
-            zav_tide(:,:,jk) = zav_tide(:,:,jk) * MIN( zkz(:,:), 30./6. ) * wmask(:,:,jk)  !kz max = 300 cm2/s
-         END DO
-         ztpc = 0._wp
-         zpc(:,:,:) = Max(0.e0,rn2(:,:,:)) * zav_tide(:,:,:)
-         DO jk= 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ztpc = ztpc + e3w_n(ji,jj,jk) * e1e2t(ji,jj) * zpc(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 * 1/(rn_tfe * rn_me) * ztpc
-         WRITE(numout,*) '          2 Total power consumption of the tidally driven part of Kz : ztpc = ', ztpc * 1.e-12 ,'TW'
-!!gm bug mpp  in these diagnostics
-         DO jk = 1, jpk
-            ze_z =                  SUM( e1e2t(:,:) * zav_tide(:,:,jk) * tmask_i(:,:) )   &
-               &     / MAX( 1.e-20, SUM( e1e2t(:,:) * wmask   (:,:,jk) * tmask_i(:,:) ) )
-            ztpc = 1.e50
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  IF( zav_tide(ji,jj,jk) /= 0.e0 )   ztpc = MIN( ztpc, zav_tide(ji,jj,jk) )
-               END DO
-            END DO
-            WRITE(numout,*) '            N2 min - jk= ', jk,'   ', ze_z * 1.e4,' cm2/s min= ',ztpc*1.e4,   &
-               &       'max= ', MAXVAL(zav_tide(:,:,jk) )*1.e4, ' cm2/s'
-         END DO
-
-         WRITE(numout,*) '          e_tide : ', SUM( e1e2t*en_tmx ) / ( rn_tfe * rn_me ) * 1.e-12, 'TW'
-         WRITE(numout,*) 
-         WRITE(numout,*) '          Initial profile of tidal vertical mixing'
-         DO jk = 1, jpk
-            DO jj = 1,jpj
-               DO ji = 1,jpi
-                  zkz(ji,jj) = az_tmx(ji,jj,jk) /MAX( rn_n2min, rn2(ji,jj,jk) )
-               END DO
-            END DO
-            ze_z =                  SUM( e1e2t(:,:) * zkz  (:,:)    * tmask_i(:,:) )   &
-               &     / MAX( 1.e-20, SUM( e1e2t(:,:) * wmask(:,:,jk) * tmask_i(:,:) ) )
-            WRITE(numout,*) '                jk= ', jk,'   ', ze_z * 1.e4,' cm2/s'
-         END DO
-         DO jk = 1, jpk
-            zkz(:,:) = az_tmx(:,:,jk) /rn_n2min
-            ze_z =                  SUM( e1e2t(:,:) * zkz  (:,:)    * tmask_i(:,:) )   &
-               &     / MAX( 1.e-20, SUM( e1e2t(:,:) * wmask(:,:,jk) * tmask_i(:,:) ) )
-            WRITE(numout,*) 
-            WRITE(numout,*) '          N2 min - jk= ', jk,'   ', ze_z * 1.e4,' cm2/s min= ',MINVAL(zkz)*1.e4,   &
-               &       'max= ', MAXVAL(zkz)*1.e4, ' cm2/s'
-         END DO
-!!gm  end bug mpp
-         !
-      ENDIF
-      !
-      CALL wrk_dealloc( jpi,jpj, zem2, zek1, zkz, zfact, zhdep )
-      CALL wrk_dealloc( jpi,jpj,jpk, zpc, zav_tide )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('zdf_tmx_init')
-      !
-   END SUBROUTINE zdf_tmx_init
-
-#elif defined key_zdftmx_new
-   !!----------------------------------------------------------------------
-   !!   'key_zdftmx_new'               Internal wave-driven vertical mixing
+
    !!----------------------------------------------------------------------
    !!   zdf_tmx       : global     momentum & tracer Kz with wave induced Kz
@@ -569,6 +36,4 @@
    PUBLIC   zdf_tmx_init    ! called in nemogcm module 
    PUBLIC   zdf_tmx_alloc   ! called in nemogcm module
-
-   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftmx = .TRUE.    !: wave-driven mixing flag
 
    !                       !!* Namelist  namzdf_tmx : internal wave-driven mixing *
@@ -1027,18 +492,4 @@
    END SUBROUTINE zdf_tmx_init
 
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option          Dummy module                NO Tidal MiXing
-   !!----------------------------------------------------------------------
-   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftmx = .FALSE.   !: tidal mixing flag
-CONTAINS
-   SUBROUTINE zdf_tmx_init           ! Dummy routine
-      WRITE(*,*) 'zdf_tmx: You should not have seen this print! error?'
-   END SUBROUTINE zdf_tmx_init
-   SUBROUTINE zdf_tmx( kt )          ! Dummy routine
-      WRITE(*,*) 'zdf_tmx: You should not have seen this print! error?', kt
-   END SUBROUTINE zdf_tmx
-#endif
-
    !!======================================================================
 END MODULE zdftmx
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 7953)
@@ -55,5 +55,5 @@
    USE ldfdyn         ! lateral viscosity setting      (ldfdyn_init routine)
    USE ldftra         ! lateral diffusivity setting    (ldftra_init routine)
-   USE zdfini         ! vertical physics setting          (zdf_init routine)
+!!gm   USE zdfphy         ! vertical physics manager      (zdf_phy_init routine)
    USE trdini         ! dyn/tra trends initialization     (trd_init routine)
    USE asminc         ! assimilation increments     
@@ -429,5 +429,5 @@
       IF( ln_ctl        )   CALL prt_ctl_init   ! Print control
       
-      CALL diurnal_sst_bulk_init            ! diurnal sst
+      CALL diurnal_sst_bulk_init             ! diurnal sst
       IF ( ln_diurnal ) CALL diurnal_sst_coolskin_init   ! cool skin   
       
@@ -455,14 +455,8 @@
                             CALL     sbc_init   ! surface boundary conditions (including sea-ice)
                             CALL     bdy_init   ! Open boundaries initialisation
+
       !                                      ! Ocean physics
-      !                                         ! Vertical physics
-                            CALL     zdf_init      ! namelist read
-                            CALL zdf_bfr_init      ! bottom friction
-      IF( lk_zdfric     )   CALL zdf_ric_init      ! Richardson number dependent Kz
-      IF( lk_zdftke     )   CALL zdf_tke_init      ! TKE closure scheme
-      IF( lk_zdfgls     )   CALL zdf_gls_init      ! GLS closure scheme
-      IF( lk_zdftmx     )   CALL zdf_tmx_init      ! tidal vertical mixing
-!!gm      IF( ln_zdfddm     )   CALL zdf_ddm_init      ! double diffusive mixing
-         
+                            CALL zdf_phy_init   ! Vertical physics
+                                     
       !                                         ! Lateral physics
                             CALL ldf_tra_init      ! Lateral ocean tracer physics
@@ -470,5 +464,5 @@
                             CALL ldf_dyn_init      ! Lateral ocean momentum physics
 
-      !                                         ! Active tracers
+      !                                      ! Active tracers
                             CALL tra_qsr_init      ! penetrative solar radiation qsr
                             CALL tra_bbc_init      ! bottom heat flux
@@ -479,5 +473,5 @@
                             CALL tra_zdf_init      ! vertical mixing and after tracer fields
 
-      !                                         ! Dynamics
+      !                                      ! Dynamics
       IF( lk_c1d        )   CALL dyn_dmp_init      ! internal momentum damping
                             CALL dyn_adv_init      ! advection (vector or flux form)
@@ -511,5 +505,5 @@
       IF( ln_diaobs     )   CALL dia_obs( nit000 - 1 )   ! Observation operator for restart
 
-      !                                         ! Assimilation increments
+      !                                      ! Assimilation increments
       IF( lk_asminc     )   CALL asm_inc_init   ! Initialize assimilation increments
       IF(lwp) WRITE(numout,*) 'Euler time step switch is ', neuler
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 7953)
@@ -74,7 +74,7 @@
       !!              -8- Outputs and diagnostics
       !!----------------------------------------------------------------------
-      INTEGER ::   ji,jj,jk ! dummy loop indice
-      INTEGER ::   indic    ! error indicator if < 0
-      INTEGER ::   kcall    ! optional integer argument (dom_vvl_sf_nxt)
+      INTEGER ::   ji, jj, jk   ! dummy loop indice
+      INTEGER ::   indic        ! error indicator if < 0
+      INTEGER ::   kcall        ! optional integer argument (dom_vvl_sf_nxt)
       !! ---------------------------------------------------------------------
 #if defined key_agrif
@@ -125,40 +125,7 @@
                          CALL bn2    ( tsn, rab_n, rn2  ) ! now    Brunt-Vaisala frequency
 
-      !
       !  VERTICAL PHYSICS
-                         CALL zdf_bfr( kstp )         ! bottom friction (if quadratic)
-      !                                               ! Vertical eddy viscosity and diffusivity coefficients
-      IF( lk_zdfric  )   CALL zdf_ric ( kstp )             ! Richardson number dependent Kz
-      IF( lk_zdftke  )   CALL zdf_tke ( kstp )             ! TKE closure scheme for Kz
-      IF( lk_zdfgls  )   CALL zdf_gls ( kstp )             ! GLS closure scheme for Kz
-      IF( ln_zdfqiao )   CALL zdf_qiao( kstp )             ! Qiao vertical mixing 
-      !
-      IF( lk_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( lk_zdftmx  )   CALL zdf_tmx( kstp )         ! tidal vertical mixing
-
-                         CALL zdf_mxl( kstp )         ! mixed layer depth
-
-                                                      ! write TKE or GLS information in the restart file
-      IF( lrst_oce .AND. lk_zdftke )   CALL tke_rst( kstp, 'WRITE' )
-      IF( lrst_oce .AND. lk_zdfgls )   CALL gls_rst( kstp, 'WRITE' )
-      !
+                         CALL zdf_phy( kstp )         ! vertical physics update (bfr, avt, avs, avm + MLD)
+
       !  LATERAL  PHYSICS
       !
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 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 7953)
@@ -63,4 +63,6 @@
    USE ldftra           ! lateral eddy diffusive coef.     (ldf_tra routine)
 
+   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)
@@ -72,4 +74,5 @@
    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/TRP/trcdmp.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/TOP_SRC/TRP/trcdmp.F90	(revision 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/TOP_SRC/TRP/trcdmp.F90	(revision 7953)
@@ -121,5 +121,5 @@
                      DO jj = 2, jpjm1
                         DO ji = fs_2, fs_jpim1   ! vector opt.
-                           IF( avt(ji,jj,jk) <= 5.e-4_wp )  THEN 
+                           IF( avs(ji,jj,jk) <= 5.e-4_wp )  THEN 
                               tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + restotr(ji,jj,jk) * ( ztrcdta(ji,jj,jk) - trb(ji,jj,jk,jn) )
                            ENDIF
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/TOP_SRC/TRP/trczdf.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/TOP_SRC/TRP/trczdf.F90	(revision 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/TOP_SRC/TRP/trczdf.F90	(revision 7953)
@@ -4,6 +4,7 @@
    !! Ocean Passive tracers : vertical diffusive trends 
    !!=====================================================================
-   !! History :  9.0  ! 2005-11 (G. Madec)  Original code
+   !! History :  9.0  ! 2005-11  (G. Madec)  Original code
    !!       NEMO 3.0  ! 2008-01  (C. Ethe, G. Madec)  merge TRC-TRA 
+   !!            4.0  ! 2017-04  (G. Madec)  remove the explicit case
    !!----------------------------------------------------------------------
 #if defined key_top
@@ -11,6 +12,5 @@
    !!   'key_top'                                                TOP models
    !!----------------------------------------------------------------------
-   !!   trc_zdf      : update the tracer trend with the lateral diffusion
-   !!   trc_zdf_ini  : initialization, namelist read, and parameters control
+   !!   trc_zdf      : update the tracer trend with the vertical diffusion
    !!----------------------------------------------------------------------
    USE trc           ! ocean passive tracers variables
@@ -27,14 +27,5 @@
 
    PUBLIC   trc_zdf         ! called by step.F90 
-   PUBLIC   trc_zdf_ini     ! called by nemogcm.F90 
    
-   !                                        !!** Vertical diffusion (nam_trczdf) **
-   LOGICAL , PUBLIC ::   ln_trczdf_exp       !: explicit vertical diffusion scheme flag
-   INTEGER , PUBLIC ::   nn_trczdf_exp       !: number of sub-time step (explicit time stepping)
-
-   INTEGER ::   nzdf = 0               ! type vertical diffusion algorithm used
-      !                                ! defined from ln_zdf...  namlist logicals)
-   !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 3.7 , NEMO Consortium (2015)
@@ -48,5 +39,6 @@
       !!                  ***  ROUTINE trc_zdf  ***
       !!
-      !! ** Purpose :   compute the vertical ocean tracer physics.
+      !! ** Purpose :   compute the vertical ocean tracer physics using
+      !!              an implicit time-stepping scheme.
       !!---------------------------------------------------------------------
       INTEGER, INTENT( in ) ::  kt      ! ocean time-step index
@@ -54,19 +46,13 @@
       INTEGER               ::  jk, jn
       CHARACTER (len=22)    :: charout
-      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   ztrtrd   ! 4D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra) ::   ztrtrd   ! 4D workspace
       !!---------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('trc_zdf')
       !
-      IF( l_trdtrc )  THEN
-         CALL wrk_alloc( jpi, jpj, jpk, jptra, ztrtrd )
-         ztrtrd(:,:,:,:)  = tra(:,:,:,:)
-      ENDIF
-
-      SELECT CASE ( nzdf )                       ! compute lateral mixing trend and add it to the general trend
-      CASE ( 0 ) ;  CALL tra_zdf_exp( kt, nittrc000, 'TRC', r2dttrc, nn_trczdf_exp, trb, tra, jptra )    !   explicit scheme 
-      CASE ( 1 ) ;  CALL tra_zdf_imp( kt, nittrc000, 'TRC', r2dttrc,                trb, tra, jptra )    !   implicit scheme          
-      END SELECT
-
+      IF( l_trdtrc )   ztrtrd(:,:,:,:)  = tra(:,:,:,:)
+      !
+      CALL tra_zdf_imp( kt, nittrc000, 'TRC', r2dttrc, trb, tra, jptra )    !   implicit scheme          
+      !
       IF( l_trdtrc )   THEN                      ! save the vertical diffusive trends for further diagnostics
          DO jn = 1, jptra
@@ -76,10 +62,10 @@
             CALL trd_tra( kt, 'TRC', jn, jptra_zdf, ztrtrd(:,:,:,jn) )
          END DO
-         CALL wrk_dealloc( jpi, jpj, jpk, jptra, ztrtrd )
       ENDIF
       !                                          ! print mean trends (used for debugging)
       IF( ln_ctl )   THEN
-         WRITE(charout, FMT="('zdf ')") ;  CALL prt_ctl_trc_info(charout)
-                                           CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
+         WRITE(charout, FMT="('zdf ')")
+         CALL prt_ctl_trc_info(charout)
+         CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
       END IF
       !
@@ -87,64 +73,4 @@
       !
    END SUBROUTINE trc_zdf
-
-
-   SUBROUTINE trc_zdf_ini
-      !!----------------------------------------------------------------------
-      !!                 ***  ROUTINE trc_zdf_ini  ***
-      !!
-      !! ** Purpose :   Choose the vertical mixing scheme
-      !!
-      !! ** Method  :   Set nzdf from ln_zdfexp
-      !!      nzdf = 0   explicit (time-splitting) scheme (ln_trczdf_exp=T)
-      !!           = 1   implicit (euler backward) scheme (ln_trczdf_exp=F)
-      !!      NB: The implicit scheme is required when using : 
-      !!             - rotated lateral mixing operator
-      !!             - TKE, GLS vertical mixing scheme
-      !!----------------------------------------------------------------------
-      INTEGER ::  ios                 ! Local integer output status for namelist read
-      !!
-      NAMELIST/namtrc_zdf/ ln_trczdf_exp  , nn_trczdf_exp
-      !!----------------------------------------------------------------------
-      !
-      REWIND( numnat_ref )             ! namtrc_zdf in reference namelist 
-      READ  ( numnat_ref, namtrc_zdf, IOSTAT = ios, ERR = 905)
-905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrc_zdf in reference namelist', lwp )
-      !
-      REWIND( numnat_cfg )             ! namtrc_zdf in configuration namelist 
-      READ  ( numnat_cfg, namtrc_zdf, IOSTAT = ios, ERR = 906 )
-906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrc_zdf in configuration namelist', lwp )
-      IF(lwm) WRITE ( numont, namtrc_zdf )
-      !
-      IF(lwp) THEN                     ! Control print
-         WRITE(numout,*)
-         WRITE(numout,*) '   Namelist namtrc_zdf : set vertical diffusion  parameters'
-         WRITE(numout,*) '      time splitting / backward scheme ln_trczdf_exp = ', ln_trczdf_exp
-         WRITE(numout,*) '      number of time step              nn_trczdf_exp = ', nn_trczdf_exp
-      ENDIF
-
-      !                                ! Define the vertical tracer physics scheme
-      IF( ln_trczdf_exp ) THEN   ;   nzdf = 0     ! explicit scheme
-      ELSE                       ;   nzdf = 1     ! implicit scheme
-      ENDIF
-
-      !                                ! Force implicit schemes
-      IF( ln_trcldf_iso              )   nzdf = 1      ! iso-neutral lateral physics
-      IF( ln_trcldf_hor .AND. ln_sco )   nzdf = 1      ! horizontal lateral physics in s-coordinate
-#if defined key_zdftke || defined key_zdfgls 
-                                         nzdf = 1      ! TKE or GLS physics       
-#endif
-      IF( ln_trczdf_exp .AND. nzdf == 1 )  & 
-         CALL ctl_stop( 'trc_zdf : If using the rotated lateral mixing operator or TKE, GLS vertical scheme ', &
-            &           '          the implicit scheme is required, set ln_trczdf_exp = .false.' )
-
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'trc:zdf_ctl : vertical passive tracer physics scheme'
-         WRITE(numout,*) '~~~~~~~~~~~'
-         IF( nzdf ==  0 )   WRITE(numout,*) '              Explicit time-splitting scheme'
-         IF( nzdf ==  1 )   WRITE(numout,*) '              Implicit (euler backward) scheme'
-      ENDIF
-      !
-   END SUBROUTINE trc_zdf_ini
    
 #else
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/TOP_SRC/oce_trc.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/TOP_SRC/oce_trc.F90	(revision 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/TOP_SRC/oce_trc.F90	(revision 7953)
@@ -101,8 +101,5 @@
 
    !* vertical diffusion *
-   USE zdf_oce , ONLY :   avt        =>   avt         !: vert. diffusivity coef. at w-point for temp  
-# if defined key_zdfddm
-   USE zdfddm  , ONLY :   avs        =>   avs         !: salinity vertical diffusivity coeff. at w-point
-# endif
+   USE zdf_oce , ONLY :   avs        =>   avs         !: vert. diffusivity coef. for salinity (w-point)
 
    !* mixing & mixed layer depth *
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/TOP_SRC/trc.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/TOP_SRC/trc.F90	(revision 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/TOP_SRC/trc.F90	(revision 7953)
@@ -126,9 +126,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::  vn_tm       !: j-horizontal velocity average     [m/s]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::  tsn_tm      !: t/s average     [m/s]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::  avt_tm      !: vertical diffusivity coeff. at  w-point   [m2/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::  avs_tm      !: vertical diffusivity coeff. at  w-point   [m2/s]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::  rhop_tm     !: 
-# if defined key_zdfddm
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::  avs_tm      !: vertical double diffusivity coeff. at w-point   [m/s]
-# endif
 #if defined key_trabbl
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::  ahu_bbl_tm  !: u-, w-points
@@ -154,5 +151,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::  tsn_temp
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::  un_temp,vn_temp,wn_temp     !: hold current values of avt, un, vn, wn
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::  avt_temp, rhop_temp     !: hold current values of avt, un, vn, wn
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::  avs_temp, rhop_temp     !: hold current values of avt, un, vn, wn
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::  sshn_temp, sshb_temp, ssha_temp, rnf_temp,h_rnf_temp
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::  hdivn_temp, rotn_temp
@@ -165,7 +162,4 @@
 #endif
    !
-# if defined key_zdfddm
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::  avs_temp      !: salinity vertical diffusivity coeff. at w-point   [m/s]
-# endif
    !
    CHARACTER(len=20), PUBLIC, ALLOCATABLE,  SAVE,  DIMENSION(:)   ::  cn_trc_dflt          ! Default OBC condition for all tracers
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/TOP_SRC/trcini.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 7953)
@@ -196,5 +196,4 @@
       USE trcadv , ONLY:  trc_adv_ini
       USE trcldf , ONLY:  trc_ldf_ini
-      USE trczdf , ONLY:  trc_zdf_ini
       USE trcrad , ONLY:  trc_rad_ini
       !
@@ -205,5 +204,5 @@
                        CALL  trc_adv_ini          ! advection
                        CALL  trc_ldf_ini          ! lateral diffusion
-                       CALL  trc_zdf_ini          ! vertical diffusion
+                       !                          ! vertical diffusion: always implicit time stepping scheme
                        CALL  trc_rad_ini          ! positivity of passive tracers 
       !
@@ -223,16 +222,16 @@
       !!----------------------------------------------------------------------
       !
-      ! Initialisation of tracers Initial Conditions
-      IF( ln_trcdta )      CALL trc_dta_ini(jptra)
-
-      ! Initialisation of tracers Boundary Conditions
-      IF( ln_my_trc )     CALL trc_bc_ini(jptra)
-
-      IF( ln_rsttr ) THEN
+
+      IF( ln_trcdta )   CALL trc_dta_ini( jptra )      ! set initial tracers values
+
+      IF( ln_my_trc )   CALL trc_bc_ini ( jptra )      ! set tracers Boundary Conditions
+
+
+      IF( ln_rsttr ) THEN              ! restart from a file
         !
-        CALL trc_rst_read              ! restart from a file
+        CALL trc_rst_read
         !
-      ELSE
-        ! Initialisation of tracer from a file that may also be used for damping
+      ELSE                             ! Initialisation of tracer from a file that may also be used for damping
+!!gm BUG ?   if damping and restart, what's happening ?
         IF( ln_trcdta .AND. nb_trcdta > 0 ) THEN
             ! update passive tracers arrays with input data read from file
@@ -250,5 +249,5 @@
                   ENDIF
                ENDIF
-            ENDDO
+            END DO
             !
         ENDIF
@@ -262,4 +261,5 @@
    END SUBROUTINE trc_ini_state
 
+
    SUBROUTINE top_alloc
       !!----------------------------------------------------------------------
Index: /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/TOP_SRC/trcsub.F90
===================================================================
--- /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/TOP_SRC/trcsub.F90	(revision 7952)
+++ /branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/TOP_SRC/trcsub.F90	(revision 7953)
@@ -84,8 +84,5 @@
           tsn_tm  (:,:,:,jp_sal) = tsn_tm  (:,:,:,jp_sal) + tsn  (:,:,:,jp_sal) * e3t_n(:,:,:)  
           rhop_tm (:,:,:)        = rhop_tm (:,:,:)        + rhop (:,:,:)        * e3t_n(:,:,:)  
-          avt_tm  (:,:,:)        = avt_tm  (:,:,:)        + avt  (:,:,:)        * e3w_n(:,:,:)  
-# if defined key_zdfddm
           avs_tm  (:,:,:)        = avs_tm  (:,:,:)        + avs  (:,:,:)        * e3w_n(:,:,:)  
-# endif
          IF( l_ldfslp ) THEN
             uslp_tm (:,:,:)      = uslp_tm (:,:,:)        + uslp (:,:,:)
@@ -122,8 +119,5 @@
          tsn_temp   (:,:,:,:)    = tsn   (:,:,:,:)
          rhop_temp  (:,:,:)      = rhop  (:,:,:)    
-         avt_temp   (:,:,:)      = avt   (:,:,:)
-# if defined key_zdfddm
          avs_temp   (:,:,:)      = avs   (:,:,:)
-# endif
          IF( l_ldfslp ) THEN
             uslp_temp  (:,:,:)   = uslp  (:,:,:)   ;   wslpi_temp (:,:,:)   = wslpi (:,:,:)
@@ -161,8 +155,5 @@
          tsn_tm   (:,:,:,jp_sal) = tsn_tm  (:,:,:,jp_sal) + tsn  (:,:,:,jp_sal) * e3t_n(:,:,:)  
          rhop_tm (:,:,:)         = rhop_tm (:,:,:)        + rhop (:,:,:)        * e3t_n(:,:,:)  
-         avt_tm   (:,:,:)        = avt_tm  (:,:,:)        + avt  (:,:,:)        * e3w_n(:,:,:)  
-# if defined key_zdfddm
          avs_tm   (:,:,:)        = avs_tm  (:,:,:)        + avs  (:,:,:)        * e3w_n(:,:,:)  
-# endif
          IF( l_ldfslp ) THEN
             uslp_tm  (:,:,:)     = uslp_tm (:,:,:)        + uslp (:,:,:) 
@@ -245,9 +236,6 @@
                   tsn  (ji,jj,jk,jp_sal) = tsn_tm  (ji,jj,jk,jp_sal) * z1_ne3t
                   rhop (ji,jj,jk)        = rhop_tm (ji,jj,jk)        * z1_ne3t
-!!gm : BUG? ==>> for avt & avs I don't understand the division by e3w
-                  avt  (ji,jj,jk)        = avt_tm  (ji,jj,jk)        * z1_ne3w
-# if defined key_zdfddm
+!!gm : BUG ==>> for avs I don't understand the division by e3w
                   avs  (ji,jj,jk)        = avs_tm  (ji,jj,jk)        * z1_ne3w
-# endif
                END DO
             END DO
@@ -297,8 +285,5 @@
       rhop_tm (:,:,:)        = rhop (:,:,:)        * e3t_n(:,:,:)  
 !!gm : BUG? ==>> for avt & avs I don't understand the division by e3w
-      avt_tm  (:,:,:)        = avt  (:,:,:)        * e3w_n(:,:,:)  
-# if defined key_zdfddm
       avs_tm  (:,:,:)        = avs  (:,:,:)        * e3w_n(:,:,:)  
-# endif
       IF( l_ldfslp ) THEN
          wslpi_tm(:,:,:)     = wslpi(:,:,:)
@@ -354,8 +339,5 @@
       tsn   (:,:,:,:) =  tsn_temp   (:,:,:,:)
       rhop  (:,:,:)   =  rhop_temp  (:,:,:)
-      avt   (:,:,:)   =  avt_temp   (:,:,:)
-# if defined key_zdfddm
       avs   (:,:,:)   =  avs_temp   (:,:,:)
-# endif
       IF( l_ldfslp ) THEN
          wslpi (:,:,:)=  wslpi_temp (:,:,:)
@@ -396,8 +378,5 @@
          tsn_tm  (:,:,:,jp_sal) = tsn  (:,:,:,jp_sal) * e3t_n(:,:,:)  
          rhop_tm (:,:,:)        = rhop (:,:,:)        * e3t_n(:,:,:)  
-         avt_tm  (:,:,:)        = avt  (:,:,:)        * e3w_n(:,:,:)  
-# if defined key_zdfddm
          avs_tm  (:,:,:)        = avs  (:,:,:)        * e3w_n(:,:,:)  
-# endif
       IF( l_ldfslp ) THEN
          uslp_tm (:,:,:)        = uslp (:,:,:)
@@ -534,5 +513,5 @@
       !
       ALLOCATE( un_temp(jpi,jpj,jpk)        ,  vn_temp(jpi,jpj,jpk)  ,   &
-         &      wn_temp(jpi,jpj,jpk)        ,  avt_temp(jpi,jpj,jpk) ,   &
+         &      wn_temp(jpi,jpj,jpk)        ,                            &
          &      rhop_temp(jpi,jpj,jpk)      ,  rhop_tm(jpi,jpj,jpk) ,   &
          &      sshn_temp(jpi,jpj)          ,  sshb_temp(jpi,jpj) ,      &
@@ -548,10 +527,7 @@
          &      fr_i_temp(jpi,jpj)          ,  fr_i_tm(jpi,jpj) ,        &
          &      wndm_temp(jpi,jpj)          ,  wndm_tm(jpi,jpj) ,        &
-# if defined key_zdfddm
          &      avs_tm(jpi,jpj,jpk)         ,  avs_temp(jpi,jpj,jpk) ,   &
-# endif
          &      hdivn_temp(jpi,jpj,jpk)     ,  hdivb_temp(jpi,jpj,jpk),  &
          &      un_tm(jpi,jpj,jpk)          ,  vn_tm(jpi,jpj,jpk)  ,     &
-         &      avt_tm(jpi,jpj,jpk)                                ,     &
          &      sshn_tm(jpi,jpj)            ,  sshb_hold(jpi,jpj) ,      &
          &      tsn_tm(jpi,jpj,jpk,2)       ,                            &
