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 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/AMM12/EXP00/namelist_cfg	(revision 8143)
@@ -174,13 +174,25 @@
    filtide      = 'bdydta/amm12_bdytide_'         !  file name root of tidal forcing files
 /
-!-----------------------------------------------------------------------
-&nambfr        !   bottom friction
-!-----------------------------------------------------------------------
-   nn_bfr      =    2      !  type of bottom friction :   = 0 : free slip,  = 1 : linear friction
-                           !                              = 2 : nonlinear friction
-   rn_bfri2    =    2.5e-3 !  bottom drag coefficient (non linear case)
-   rn_bfeb2    =    0.0e0  !  bottom turbulent kinetic energy background  (m2/s2)
-   ln_loglayer =    .true. !  loglayer bottom friction (only effect when nn_bfr = 2)
-   rn_bfrz0    =    0.003  !  bottom roughness (only effect when ln_loglayer = .true.)
+
+!-----------------------------------------------------------------------
+&namdrg            !   top/bottom drag coefficient                      (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_NONE    = .false.    !  free-slip       : Cd = 0                  (F => fill namdrg_bot
+   ln_lin     = .false.    !      linear  drag: Cd = Cd0 Uc0                   &   namdrg_top)
+   ln_non_lin = .false.    !  non-linear  drag: Cd = Cd0 |U|
+   ln_loglayer= .true.     !  logarithmic drag: Cd = vkarmn/log(z/z0) |U|
+   !
+   ln_drgimp  = .true.     !  implicit top/bottom friction flag
+/
+!-----------------------------------------------------------------------
+&namdrg_bot        !   BOTTOM friction                                  
+!-----------------------------------------------------------------------
+   rn_Cd0     =  2.5e-3   !  drag coefficient [-]
+   rn_Uc0     =  0.4      !  ref. velocity [m/s] (linear drag=Cd0*Uc0) 
+   rn_Cdmax   =  0.1      !  drag value maximum [-] (logarithmic drag)
+   rn_ke0     =  0.0e0    !  background kinetic energy  [m2/s2] (non-linear cases)
+   rn_z0      =  0.003    !  roughness [m]  (ln_loglayer=T)
+   ln_boost   = .false.   !  =T regional boost of Cd0 ; =F constant
+      rn_boost=  50.         !  local boost factor  [-]
 /
 !-----------------------------------------------------------------------
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/C1D_PAPA/EXP00/namelist_cfg
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/C1D_PAPA/EXP00/namelist_cfg	(revision 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/C1D_PAPA/EXP00/namelist_cfg	(revision 8143)
@@ -183,7 +183,7 @@
 /
 !-----------------------------------------------------------------------
-&nambfr        !   bottom friction
-!-----------------------------------------------------------------------
-   nn_bfr      =    2      !  type of bottom friction :   = 0 : free slip,  = 1 : linear friction
+&namdrg            !   top/bottom drag coefficient                      (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_non_lin = .true.     !  non-linear  drag: Cd = Cd0 |U|
 /
 !-----------------------------------------------------------------------
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/GYRE_BFM/EXP00/namelist_cfg
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/GYRE_BFM/EXP00/namelist_cfg	(revision 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/GYRE_BFM/EXP00/namelist_cfg	(revision 8143)
@@ -124,7 +124,7 @@
 /
 !-----------------------------------------------------------------------
-&nambfr        !   bottom friction
-!-----------------------------------------------------------------------
-   nn_bfr      =    2      !  type of bottom friction :   = 0 : free slip,  = 1 : linear friction
+&namdrg            !   top/bottom drag coefficient                      (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_non_lin = .true.     !  non-linear  drag: Cd = Cd0 |U|
 /
 !-----------------------------------------------------------------------
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 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/GYRE_PISCES/EXP00/namelist_cfg	(revision 8143)
@@ -78,7 +78,7 @@
 /
 !-----------------------------------------------------------------------
-&nambfr        !   bottom friction
-!-----------------------------------------------------------------------
-   nn_bfr      =    2      !  type of bottom friction :   = 0 : free slip,  = 1 : linear friction
+&namdrg        !   top/bottom friction
+!-----------------------------------------------------------------------
+   ln_non_lin = .true.     !  non-linear  drag: Cd = Cd0 |U|
 /
 !-----------------------------------------------------------------------
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/ORCA2_LIM3_PISCES/EXP00/1_namelist_cfg
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/ORCA2_LIM3_PISCES/EXP00/1_namelist_cfg	(revision 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/ORCA2_LIM3_PISCES/EXP00/1_namelist_cfg	(revision 8143)
@@ -86,6 +86,7 @@
 /
 !-----------------------------------------------------------------------
-&nambfr        !   bottom friction
-!-----------------------------------------------------------------------
+&namdrg        !   bottom friction
+!-----------------------------------------------------------------------
+   ln_lin = .true.         !     linear  drag: Cd = Cd0 Uc0
 /
 !-----------------------------------------------------------------------
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 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/ORCA2_LIM3_PISCES/EXP00/namelist_cfg	(revision 8143)
@@ -99,6 +99,7 @@
 /
 !-----------------------------------------------------------------------
-&nambfr        !   bottom friction
-!-----------------------------------------------------------------------
+&namdrg        !   top/bottom friction
+!-----------------------------------------------------------------------
+   ln_lin = .true.         !     linear  drag: Cd = Cd0 Uc0
 /
 !-----------------------------------------------------------------------
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 8143)
@@ -8,5 +8,5 @@
 !!                                    namsbc_apr, namsbc_ssr, namsbc_alb, namsbc_wave)
 !!              4 - lateral boundary (namlbc, namagrif, nambdy, nambdy_tide)
-!!              5 - bottom  boundary (nambfr, nambbc, nambbl)
+!!              5 - bottom  boundary (namdrg, namdrg_top, namdrg_bot, nambbc, nambbl)
 !!              6 - Tracer           (nameos, namtra_adv, namtra_ldf, namtra_ldfeiv, namtra_dmp)
 !!              7 - dynamics         (namdyn_adv, namdyn_vor, namdyn_hpg, namdyn_spg, namdyn_ldf)
@@ -603,6 +603,7 @@
 !!                ***  top/Bottom boundary condition  ***             !!
 !!======================================================================
-!!   nambfr        bottom friction                                      (default: NONE)
-!!   namtfr        top    friction                                      (default: NONE)
+!!   namdrg        top/bottom drag coefficient                          (default: NONE)
+!!   namdrg_top    top    friction                                      (ln_isfcav=T)
+!!   namdrg_bot    bottom friction                                      
 !!   nambbc        bottom temperature boundary condition                (default: NO)
 !!   nambbl        bottom boundary layer scheme                         (default: NO)
@@ -610,30 +611,5 @@
 !
 !-----------------------------------------------------------------------
-&nambfr        !   bottom friction                                      (default: linear)
-!-----------------------------------------------------------------------
-   nn_bfr      =    1      !  type of top/bottom drag: free slip (=0), linear drag (=1), 
-   !                       !  nonlinear drag (=2), nonlinear with logarithmic formulation (=3)
-   ln_bfrimp   = .true.    !  implicit bottom friction (requires ln_zdfexp = .false. if true)
-   ln_loglayer = .false.   !  logarithmic formulation (non linear case only)
-                           !
-   rn_bfri1    =    4.e-4  !  bottom drag coefficient (linear case)
-   rn_bfri2    =    1.e-3  !  bottom drag coefficient (non linear case). Minimum coeft if ln_loglayer=T
-   rn_bfri2_max=    1.e-1  !  max. bottom drag coefficient (non linear case and ln_loglayer=T)
-   rn_bfeb2    =    2.5e-3 !  bottom turbulent kinetic energy background  (m2/s2)
-   rn_bfrz0    =    3.e-3  !  bottom roughness [m] if ln_loglayer=T
-   ln_bfr2d    = .false.   !  horizontal variation of the bottom friction coef (read a 2D mask file )
-      rn_bfrien   =    50.    !  local boost factor
-   !
-   rn_tfri1    =    4.e-4  !  top drag coefficient (linear case)
-   rn_tfri2    =    2.5e-3 !  top drag coefficient (non linear case). Minimum coeft if ln_loglayer=T
-   rn_tfri2_max=    1.e-1  !  max. top drag coefficient (non linear case and ln_loglayer=T)
-   rn_tfeb2    =    0.0    !  top turbulent kinetic energy background  (m2/s2)
-   rn_tfrz0    =    3.e-3  !  top roughness [m] if ln_loglayer=T
-   ln_tfr2d    = .false.   !  horizontal variation of the top friction coef (read a 2D mask file )
-      rn_tfrien   =   50.     !  local boost factor
-/
-
-!-----------------------------------------------------------------------
-&namdrg            !   top/bottom drag coeeficient                      (default: NO selection)
+&namdrg            !   top/bottom drag coefficient                      (default: NO selection)
 !-----------------------------------------------------------------------
    ln_NONE    = .false.    !  free-slip       : Cd = 0                  (F => fill namdrg_bot
@@ -644,5 +620,4 @@
    ln_drgimp  = .true.     !  implicit top/bottom friction flag
 /
-
 !-----------------------------------------------------------------------
 &namdrg_top        !   TOP friction                                     (ln_isfcav=T)
@@ -656,5 +631,4 @@
       rn_boost=  50.          !  local boost factor  [-]
 /
-
 !-----------------------------------------------------------------------
 &namdrg_bot        !   BOTTOM friction                                  
@@ -668,5 +642,4 @@
       rn_boost=  50.         !  local boost factor  [-]
 /
-
 !-----------------------------------------------------------------------
 &nambbc        !   bottom temperature boundary condition                (default: NO)
@@ -980,16 +953,15 @@
    ln_mxl0     = .true.    !  surface mixing length scale = F(wind stress) (T) or not (F)
    rn_mxl0     =   0.04    !  surface  buoyancy lenght scale minimum value
-!!gm   ln_drg      = .false.   !  top/bottom friction added as boundary condition of TKE
+   ln_drg      = .false.   !  top/bottom friction added as boundary condition of TKE
    ln_lc       = .true.    !  Langmuir cell parameterisation (Axell 2002)
-   rn_lc       =   0.15    !  coef. associated to Langmuir cells
+      rn_lc       =   0.15    !  coef. associated to Langmuir cells
    nn_etau     =   1       !  penetration of tke below the mixed layer (ML) due to NIWs
-                           !        = 0 no penetration
-                           !        = 1 add a tke source below the ML
-                           !        = 2 add a tke source just at the base of the ML
-                           !        = 3 as = 1 applied on HF part of the stress           (ln_cpl=T)
-   rn_efr      =   0.05    !  fraction of surface tke value which penetrates below the ML (nn_etau=1 or 2)
-   nn_htau     =   1       !  type of exponential decrease of tke penetration below the ML
-                           !        = 0  constant 10 m length scale
-                           !        = 1  0.5m at the equator to 30m poleward of 40 degrees
+                              !        = 0 none ; = 1 add a tke source below the ML
+                              !        = 2 add a tke source just at the base of the ML
+                              !        = 3 as = 1 applied on HF part of the stress           (ln_cpl=T)
+      rn_efr      =   0.05    !  fraction of surface tke value which penetrates below the ML (nn_etau=1 or 2)
+      nn_htau     =   1       !  type of exponential decrease of tke penetration below the ML
+                              !        = 0  constant 10 m length scale
+                              !        = 1  0.5m at the equator to 30m poleward of 40 degrees
 /
 !-----------------------------------------------------------------------
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/ISOMIP/EXP00/namelist_cfg
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/ISOMIP/EXP00/namelist_cfg	(revision 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/ISOMIP/EXP00/namelist_cfg	(revision 8143)
@@ -5,5 +5,7 @@
 &namusr_def    !   ISOMIP user defined namelist  
 !-----------------------------------------------------------------------
+   ln_zco      = .false.      ! z-coordinate
    ln_zps      = .true.       ! z-partial-step coordinate
+   ln_sco      = .false.      ! s-coordinate
    rn_lam0     =   0.0        !  longitude of first raw and column T-point (jphgr_msh = 1)
    rn_phi0     = -80.0        ! latitude  of first raw and column T-point (jphgr_msh = 1)
@@ -179,25 +181,18 @@
 /
 !-----------------------------------------------------------------------
-&nambfr        !   bottom friction
-!-----------------------------------------------------------------------
-   nn_bfr      =    2      !  type of bottom friction :   = 0 : free slip,  = 1 : linear friction
-                           !                              = 2 : nonlinear friction
-   rn_bfri1    =    4.e-4  !  bottom drag coefficient (linear case)
-   rn_bfri2    =    1.e-3  !  bottom drag coefficient (non linear case). Minimum coeft if ln_loglayer=T
-   rn_bfri2_max =   1.e-1  !  max. bottom drag coefficient (non linear case and ln_loglayer=T)
-   rn_bfeb2    =    2.5e-3 !  bottom turbulent kinetic energy background  (m2/s2)
-   rn_bfrz0    =    3.e-3  !  bottom roughness [m] if ln_loglayer=T
-   ln_bfr2d    = .false.   !  horizontal variation of the bottom friction coef (read a 2D mask file )
-   rn_bfrien   =    50.    !  local multiplying factor of bfr (ln_bfr2d=T)
-   rn_tfri1    =    4.e-4  !  top drag coefficient (linear case)
-   rn_tfri2    =    2.5e-3 !  top drag coefficient (non linear case). Minimum coeft if ln_loglayer=T
-   rn_tfri2_max =   1.e-1  !  max. top drag coefficient (non linear case and ln_loglayer=T)
-   rn_tfeb2    =    0.0    !  top turbulent kinetic energy background  (m2/s2)
-   rn_tfrz0    =    3.e-3  !  top roughness [m] if ln_loglayer=T
-   ln_tfr2d    = .false.   !  horizontal variation of the top friction coef (read a 2D mask file )
-   rn_tfrien   =    50.    !  local multiplying factor of tfr (ln_tfr2d=T)
-
-   ln_bfrimp   = .true.    !  implicit bottom friction (requires ln_zdfexp = .false. if true)
-   ln_loglayer = .false.   !  logarithmic formulation (non linear case)
+&namdrg            !   top/bottom drag coefficient                      (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_non_lin = .false.    !  non-linear  drag: Cd = Cd0 |U|
+/
+!-----------------------------------------------------------------------
+&namdrg_top        !   TOP friction                                     (ln_isfcav=T)
+!-----------------------------------------------------------------------
+   rn_Cd0     =  2.5e-3    !  drag coefficient [-]
+   rn_Uc0     =  1.6       !  ref. velocity [m/s] (linear drag=Cd0*Uc0) 
+   rn_Cdmax   =  0.1       !  drag value maximum [-] (logarithmic drag)
+   rn_ke0     =  0.0e-0    !  background kinetic energy  [m2/s2] (non-linear cases)
+   rn_z0      =  3.0e-3    !  roughness [m] (ln_loglayer=T)
+   ln_boost   = .false.    !  =T regional boost of Cd0 ; =F constant
+      rn_boost=  50.          !  local boost factor  [-]
 /
 !-----------------------------------------------------------------------
@@ -215,7 +210,5 @@
 &nameos        !   ocean physical parameters
 !-----------------------------------------------------------------------
-   ln_teos10   = .false.         !  = Use TEOS-10 equation of state
    ln_eos80    = .true.          !  = Use EOS80 equation of state
-   !                             !  rd(T,S,Z)*rau0 = -a0*(1+.5*lambda*dT+mu*Z+nu*dS)*dT+b0*dS
 /
 !-----------------------------------------------------------------------
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_FCT2_flux_cen2_cfg
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_FCT2_flux_cen2_cfg	(revision 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_FCT2_flux_cen2_cfg	(revision 8143)
@@ -68,8 +68,7 @@
 /
 !-----------------------------------------------------------------------
-&nambfr        !   bottom friction
-!-----------------------------------------------------------------------
-   nn_bfr      =    0      !  type of bottom friction :   = 0 : free slip,  = 1 : linear friction
-                           !                              = 2 : nonlinear friction
+&namdrg            !   top/bottom drag coefficient                      (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_NONE    = .false.    !  free-slip       : Cd = 0                  
 /
 !-----------------------------------------------------------------------
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_FCT2_flux_ubs_cfg
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_FCT2_flux_ubs_cfg	(revision 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_FCT2_flux_ubs_cfg	(revision 8143)
@@ -68,8 +68,7 @@
 /
 !-----------------------------------------------------------------------
-&nambfr        !   bottom friction
-!-----------------------------------------------------------------------
-   nn_bfr      =    0      !  type of bottom friction :   = 0 : free slip,  = 1 : linear friction
-                           !                              = 2 : nonlinear friction
+&namdrg            !   top/bottom drag coefficient                      (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_NONE    = .false.    !  free-slip       : Cd = 0                  
 /
 !-----------------------------------------------------------------------
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_FCT2_vect_eenH_cfg
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_FCT2_vect_eenH_cfg	(revision 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_FCT2_vect_eenH_cfg	(revision 8143)
@@ -68,8 +68,7 @@
 /
 !-----------------------------------------------------------------------
-&nambfr        !   bottom friction
-!-----------------------------------------------------------------------
-   nn_bfr      =    0      !  type of bottom friction :   = 0 : free slip,  = 1 : linear friction
-                           !                              = 2 : nonlinear friction
+&namdrg            !   top/bottom drag coefficient                      (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_NONE    = .false.    !  free-slip       : Cd = 0                  
 /
 !-----------------------------------------------------------------------
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_FCT2_vect_een_cfg
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_FCT2_vect_een_cfg	(revision 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_FCT2_vect_een_cfg	(revision 8143)
@@ -68,8 +68,7 @@
 /
 !-----------------------------------------------------------------------
-&nambfr        !   bottom friction
-!-----------------------------------------------------------------------
-   nn_bfr      =    0      !  type of bottom friction :   = 0 : free slip,  = 1 : linear friction
-                           !                              = 2 : nonlinear friction
+&namdrg            !   top/bottom drag coefficient                      (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_NONE    = .false.    !  free-slip       : Cd = 0                  
 /
 !-----------------------------------------------------------------------
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_FCT2_vect_ene_cfg
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_FCT2_vect_ene_cfg	(revision 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_FCT2_vect_ene_cfg	(revision 8143)
@@ -68,8 +68,7 @@
 /
 !-----------------------------------------------------------------------
-&nambfr        !   bottom friction
-!-----------------------------------------------------------------------
-   nn_bfr      =    0      !  type of bottom friction :   = 0 : free slip,  = 1 : linear friction
-                           !                              = 2 : nonlinear friction
+&namdrg            !   top/bottom drag coefficient                      (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_NONE    = .false.    !  free-slip       : Cd = 0                  
 /
 !-----------------------------------------------------------------------
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_FCT2_vect_ens_cfg
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_FCT2_vect_ens_cfg	(revision 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_FCT2_vect_ens_cfg	(revision 8143)
@@ -68,8 +68,7 @@
 /
 !-----------------------------------------------------------------------
-&nambfr        !   bottom friction
-!-----------------------------------------------------------------------
-   nn_bfr      =    0      !  type of bottom friction :   = 0 : free slip,  = 1 : linear friction
-                           !                              = 2 : nonlinear friction
+&namdrg            !   top/bottom drag coefficient                      (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_NONE    = .false.    !  free-slip       : Cd = 0                  
 /
 !-----------------------------------------------------------------------
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_FCT4_flux_cen2_cfg
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_FCT4_flux_cen2_cfg	(revision 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_FCT4_flux_cen2_cfg	(revision 8143)
@@ -68,8 +68,7 @@
 /
 !-----------------------------------------------------------------------
-&nambfr        !   bottom friction
-!-----------------------------------------------------------------------
-   nn_bfr      =    0      !  type of bottom friction :   = 0 : free slip,  = 1 : linear friction
-                           !                              = 2 : nonlinear friction
+&namdrg            !   top/bottom drag coefficient                      (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_NONE    = .false.    !  free-slip       : Cd = 0                  
 /
 !-----------------------------------------------------------------------
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_FCT4_flux_ubs_cfg
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_FCT4_flux_ubs_cfg	(revision 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_FCT4_flux_ubs_cfg	(revision 8143)
@@ -68,8 +68,7 @@
 /
 !-----------------------------------------------------------------------
-&nambfr        !   bottom friction
-!-----------------------------------------------------------------------
-   nn_bfr      =    0      !  type of bottom friction :   = 0 : free slip,  = 1 : linear friction
-                           !                              = 2 : nonlinear friction
+&namdrg            !   top/bottom drag coefficient                      (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_NONE    = .false.    !  free-slip       : Cd = 0                  
 /
 !-----------------------------------------------------------------------
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_FCT4_vect_eenH_cfg
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_FCT4_vect_eenH_cfg	(revision 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_FCT4_vect_eenH_cfg	(revision 8143)
@@ -68,8 +68,7 @@
 /
 !-----------------------------------------------------------------------
-&nambfr        !   bottom friction
-!-----------------------------------------------------------------------
-   nn_bfr      =    0      !  type of bottom friction :   = 0 : free slip,  = 1 : linear friction
-                           !                              = 2 : nonlinear friction
+&namdrg            !   top/bottom drag coefficient                      (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_NONE    = .false.    !  free-slip       : Cd = 0                  
 /
 !-----------------------------------------------------------------------
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_FCT4_vect_een_cfg
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_FCT4_vect_een_cfg	(revision 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_FCT4_vect_een_cfg	(revision 8143)
@@ -68,8 +68,7 @@
 /
 !-----------------------------------------------------------------------
-&nambfr        !   bottom friction
-!-----------------------------------------------------------------------
-   nn_bfr      =    0      !  type of bottom friction :   = 0 : free slip,  = 1 : linear friction
-                           !                              = 2 : nonlinear friction
+&namdrg            !   top/bottom drag coefficient                      (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_NONE    = .false.    !  free-slip       : Cd = 0                  
 /
 !-----------------------------------------------------------------------
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_FCT4_vect_ene_cfg
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_FCT4_vect_ene_cfg	(revision 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_FCT4_vect_ene_cfg	(revision 8143)
@@ -68,8 +68,7 @@
 /
 !-----------------------------------------------------------------------
-&nambfr        !   bottom friction
-!-----------------------------------------------------------------------
-   nn_bfr      =    0      !  type of bottom friction :   = 0 : free slip,  = 1 : linear friction
-                           !                              = 2 : nonlinear friction
+&namdrg            !   top/bottom drag coefficient                      (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_NONE    = .false.    !  free-slip       : Cd = 0                  
 /
 !-----------------------------------------------------------------------
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_FCT4_vect_ens_cfg
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_FCT4_vect_ens_cfg	(revision 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_FCT4_vect_ens_cfg	(revision 8143)
@@ -68,8 +68,7 @@
 /
 !-----------------------------------------------------------------------
-&nambfr        !   bottom friction
-!-----------------------------------------------------------------------
-   nn_bfr      =    0      !  type of bottom friction :   = 0 : free slip,  = 1 : linear friction
-                           !                              = 2 : nonlinear friction
+&namdrg            !   top/bottom drag coefficient                      (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_NONE    = .false.    !  free-slip       : Cd = 0                  
 /
 !-----------------------------------------------------------------------
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_cfg
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_cfg	(revision 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_cfg	(revision 8143)
@@ -68,8 +68,7 @@
 /
 !-----------------------------------------------------------------------
-&nambfr        !   bottom friction
-!-----------------------------------------------------------------------
-   nn_bfr      =    0      !  type of bottom friction :   = 0 : free slip,  = 1 : linear friction
-                           !                              = 2 : nonlinear friction
+&namdrg            !   top/bottom drag coefficient                      (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_NONE    = .false.    !  free-slip       : Cd = 0                  
 /
 !-----------------------------------------------------------------------
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/OVERFLOW/EXP00/namelist_cfg
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/OVERFLOW/EXP00/namelist_cfg	(revision 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/OVERFLOW/EXP00/namelist_cfg	(revision 8143)
@@ -62,8 +62,7 @@
 /
 !-----------------------------------------------------------------------
-&nambfr        !   bottom friction
-!-----------------------------------------------------------------------
-   nn_bfr      =    0      !  type of bottom friction :   = 0 : free slip,  = 1 : linear friction
-                           !                              = 2 : nonlinear friction
+&namdrg            !   top/bottom drag coefficient                      (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_NONE    = .false.    !  free-slip       : Cd = 0                  
 /
 !-----------------------------------------------------------------------
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/SAS_BIPER/EXP00/namelist_cfg
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/SAS_BIPER/EXP00/namelist_cfg	(revision 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/SAS_BIPER/EXP00/namelist_cfg	(revision 8143)
@@ -82,6 +82,7 @@
 /
 !-----------------------------------------------------------------------
-&nambfr        !   bottom friction
-!-----------------------------------------------------------------------
+&namdrg            !   top/bottom drag coefficient                      (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_NONE    = .false.    !  free-slip       : Cd = 0                  
 /
 !-----------------------------------------------------------------------
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/WAD/EXP00/namelist_cfg
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/WAD/EXP00/namelist_cfg	(revision 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/WAD/EXP00/namelist_cfg	(revision 8143)
@@ -173,14 +173,18 @@
 /
 !-----------------------------------------------------------------------
-&nambfr        !   bottom friction
-!-----------------------------------------------------------------------
-   nn_bfr      =    2      !  type of bottom friction :   = 0 : free slip,  = 1 : linear friction
-   !rn_bfri2    =    1.e-5  !  bottom drag coefficient (non linear case). Minimum coeft if ln_loglayer=T
-   !rn_bfri2_max =   1.e-4  !  max. bottom drag coefficient (non linear case and ln_loglayer=T)
-   rn_bfri2    =    1.e-5  !  bottom drag coefficient (non linear case). Minimum coeft if ln_loglayer=T
-   rn_bfri2_max =   1.e-4  !  max. bottom drag coefficient (non linear case and ln_loglayer=T)
-   !rn_bfeb2    =    2.5e-3 !  bottom turbulent kinetic energy background  (m2/s2)
-   !rn_bfrz0    =    3.e-3  !  bottom roughness [m] if ln_loglayer=T
-   ln_loglayer = .true.    !  logarithmic formulation (non linear case)
+&namdrg            !   top/bottom drag coefficient                      (default: NO selection)
+!-----------------------------------------------------------------------
+   ln_loglayer= .false.    !  logarithmic drag: Cd = vkarmn/log(z/z0) |U|
+/
+!-----------------------------------------------------------------------
+&namdrg_bot        !   BOTTOM friction                                  
+!-----------------------------------------------------------------------
+   rn_Cd0     =  1.e-4    !  drag coefficient [-]
+   rn_Uc0     =  0.1      !  ref. velocity [m/s] (linear drag=Cd0*Uc0) 
+   rn_Cdmax   =  1.e-4    !  drag value maximum [-] (logarithmic drag)
+   rn_ke0     =  2.5e-3   !  background kinetic energy  [m2/s2] (non-linear cases)
+   rn_z0      =  3.e-3    !  roughness [m] (ln_loglayer=T)
+   ln_boost   = .false.   !  =T regional boost of Cd0 ; =F constant
+      rn_boost=  50.         !  local boost factor  [-]
 /
 !-----------------------------------------------------------------------
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 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/CRS/crs.F90	(revision 8143)
@@ -151,5 +151,6 @@
 
       ! Vertical diffusion
-      REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)  ::  avs_crs           !: salinity vertical diffusivity coeff. [m2/s] at w-point
+      REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)  ::  avt_crs           !: temperature vertical diffusivity coeff. [m2/s] at w-point
+      REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)  ::  avs_crs           !: salinity    vertical diffusivity coeff. [m2/s] at w-point
 
       ! Mixing and Mixed Layer Depth
@@ -235,5 +236,6 @@
          &     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), avs_crs(jpi_crs,jpj_crs,jpk), STAT=ierr(13) )
+     ALLOCATE( tsn_crs(jpi_crs,jpj_crs,jpk,jpts), avt_crs(jpi_crs,jpj_crs,jpk),   &
+         &                                        avs_crs(jpi_crs,jpj_crs,jpk), STAT=ierr(13) )
 
       ALLOCATE( nmln_crs(jpi_crs,jpj_crs) , hmld_crs(jpi_crs,jpj_crs) , &
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 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/CRS/crsfld.F90	(revision 8143)
@@ -210,16 +210,19 @@
       !  free memory
 
-      !  avt, avs
-!!gm BUG   TOP always uses avs !!!
+      !  avs
       SELECT CASE ( nn_crs_kz )
          CASE ( 0 )
-            CALL crs_dom_ope( avt, 'VOL', 'W', tmask, avs_crs, p_e12=e1e2t, p_e3=ze3w, psgn=1.0 )
+            CALL crs_dom_ope( avt, 'VOL', 'W', tmask, avt_crs, p_e12=e1e2t, p_e3=ze3w, psgn=1.0 )
+            CALL crs_dom_ope( avs, 'VOL', 'W', tmask, avs_crs, p_e12=e1e2t, p_e3=ze3w, psgn=1.0 )
          CASE ( 1 )
-            CALL crs_dom_ope( avt, 'MAX', 'W', tmask, avs_crs, p_e12=e1e2t, p_e3=ze3w, psgn=1.0 )
+            CALL crs_dom_ope( avt, 'MAX', 'W', tmask, avt_crs, p_e12=e1e2t, p_e3=ze3w, psgn=1.0 )
+            CALL crs_dom_ope( avs, 'MAX', 'W', tmask, avs_crs, p_e12=e1e2t, p_e3=ze3w, psgn=1.0 )
          CASE ( 2 )
-            CALL crs_dom_ope( avt, 'MIN', 'W', tmask, avs_crs, p_e12=e1e2t, p_e3=ze3w, psgn=1.0 )
+            CALL crs_dom_ope( avt, 'MIN', 'W', tmask, avt_crs, p_e12=e1e2t, p_e3=ze3w, psgn=1.0 )
+            CALL crs_dom_ope( avs, 'MIN', 'W', tmask, avs_crs, p_e12=e1e2t, p_e3=ze3w, psgn=1.0 )
       END SELECT
       !
-      CALL iom_put( "avt", avs_crs )   !  Kz
+      CALL iom_put( "avt", avt_crs )   !  Kz on T
+      CALL iom_put( "avs", avs_crs )   !  Kz on S
       
       !  sbc fields  
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 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 8143)
@@ -61,5 +61,4 @@
    USE diurnal_bulk    ! diurnal warm layer
    USE cool_skin       ! Cool skin
-   USE wrk_nemo        ! working array
 
    IMPLICIT NONE
@@ -183,12 +182,4 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-!!gm old
-!!gm BUG  missing x 0.5
-               zztmpx = (  bfrua(ji  ,jj) * un(ji  ,jj,mbku(ji  ,jj))  &
-                      &  + bfrua(ji-1,jj) * un(ji-1,jj,mbku(ji-1,jj))  )      
-               zztmpy = (  bfrva(ji,  jj) * vn(ji,jj  ,mbkv(ji,jj  ))  &
-                      &  + bfrva(ji,jj-1) * vn(ji,jj-1,mbkv(ji,jj-1))  ) 
-               z2d(ji,jj) = rau0 * SQRT( zztmpx * zztmpx + zztmpy * zztmpy ) * tmask(ji,jj,1) 
-!!gm
                zztmp2 = (  ( rCdU_bot(ji+1,jj)+rCdU_bot(ji  ,jj) ) * un(ji  ,jj,mbku(ji  ,jj))  )**2   &
                   &   + (  ( rCdU_bot(ji  ,jj)+rCdU_bot(ji-1,jj) ) * un(ji-1,jj,mbku(ji-1,jj))  )**2   &
@@ -196,8 +187,7 @@
                   &   + (  ( rCdU_bot(ji,jj  )+rCdU_bot(ji,jj-1) ) * vn(ji,jj-1,mbkv(ji,jj-1))  )**2
                z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1) 
-!!gm new end
                !
-            ENDDO
-         ENDDO
+            END DO
+         END DO
          CALL lbc_lnk( z2d, 'T', 1. )
          CALL iom_put( "taubot", z2d )           
@@ -449,15 +439,11 @@
       REAL(wp) ::   zsto, zout, zmax, zjulian                ! local scalars
       !
-      REAL(wp), POINTER, DIMENSION(:,:)   :: zw2d       ! 2D workspace
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d       ! 3D workspace
+      REAL(wp), DIMENSION(jpi,jpj)   :: zw2d       ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d       ! 3D workspace
       !!----------------------------------------------------------------------
       ! 
       IF( nn_timing == 1 )   CALL timing_start('dia_wri')
       !
-                             CALL wrk_alloc( jpi,jpj      , zw2d )
-      IF( .NOT.ln_linssh )   CALL wrk_alloc( jpi,jpj,jpk  , zw3d )
-      !
-      ! Output the initial state and forcings
-      IF( ninist == 1 ) THEN                       
+      IF( ninist == 1 ) THEN     !==  Output the initial state and forcings  ==!
          CALL dia_wri_state( 'output.init', kt )
          ninist = 0
@@ -467,6 +453,5 @@
       ! -----------------
 
-      ! local variable for debugging
-      ll_print = .FALSE.
+      ll_print = .FALSE.                  ! local variable for debugging
       ll_print = ll_print .AND. lwp
 
@@ -891,7 +876,4 @@
       ENDIF
       !
-                             CALL wrk_dealloc( jpi , jpj        , zw2d )
-      IF( .NOT.ln_linssh )   CALL wrk_dealloc( jpi , jpj , jpk  , zw3d )
-      !
       IF( nn_timing == 1 )   CALL timing_stop('dia_wri')
       !
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90	(revision 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90	(revision 8143)
@@ -5,6 +5,6 @@
    !!==============================================================================
    !! History :  3.2  ! 2008-11  (A. C. Coward)  Original code
-   !!            3.4  ! 2011-09  (H. Liu) Make it consistent with semi-implicit
-   !!                            Bottom friction (ln_bfrimp = .true.) 
+   !!            3.4  ! 2011-09  (H. Liu) Make it consistent with semi-implicit Bottom friction (ln_drgimp =T) 
+   !!            4.0  ! 2017-05  (G. Madec)  drag coef. defined at t-point (zdfdrg.F90)
    !!----------------------------------------------------------------------
 
@@ -14,10 +14,6 @@
    USE oce            ! ocean dynamics and tracers variables
    USE dom_oce        ! ocean space and time domain variables 
-   USE zdf_oce        ! ocean vertical physics variables
-!!gm new
+   USE zdf_oce        ! vertical physics: variables
    USE zdfdrg         ! vertical physics: top/bottom drag coef.
-!!gm old
-   USE zdfbfr         ! ocean bottom friction variables
-!!gm
    USE trd_oce        ! trends: ocean variables
    USE trddyn         ! trend manager: dynamics
@@ -26,5 +22,4 @@
    USE prtctl         ! Print control
    USE timing         ! Timing
-   USE wrk_nemo       ! Memory Allocation
 
    IMPLICIT NONE
@@ -36,5 +31,5 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 4.0 , NEMO Consortium (2017)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -47,4 +42,7 @@
       !!
       !! ** Purpose :   compute the bottom friction ocean dynamics physics.
+      !!
+      !!              only for explicit bottom friction form
+      !!              implicit bfr is implemented in dynzdf_imp
       !!
       !! ** Action  :   (ua,va)   momentum trend increased by bottom friction trend
@@ -61,72 +59,53 @@
       IF( nn_timing == 1 )  CALL timing_start('dyn_bfr')
       !
-!!gm issue: better to put the logical in step to control the call of zdf_bfr
-!!          ==> change the logical from ln_bfrimp to ln_bfr_exp !!
-      IF( .NOT.ln_bfrimp) THEN     ! only for explicit bottom friction form
-                                    ! implicit bfr is implemented in dynzdf_imp
-
 !!gm bug : time step is only rdt (not 2 rdt if euler start !)
          zm1_2dt = - 1._wp / ( 2._wp * rdt )
 
-         IF( l_trddyn ) THEN      ! trends: store the input trends
-            ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) )
-            ztrdu(:,:,:) = ua(:,:,:)
+      IF( l_trddyn ) THEN      ! trends: store the input trends
+         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) )
+         ztrdu(:,:,:) = ua(:,:,:)
             ztrdv(:,:,:) = va(:,:,:)
-         ENDIF
+      ENDIF
 
 
+      DO jj = 2, jpjm1
+         DO ji = 2, jpim1
+            ikbu = mbku(ji,jj)          ! deepest wet ocean u- & v-levels
+            ikbv = mbkv(ji,jj)
+            !
+            ! Apply stability criteria on absolute value  : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt)
+            zCdu = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / e3u_n(ji,jj,ikbu)
+            zCdv = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / e3v_n(ji,jj,ikbv)
+            !
+            ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX(  zCdu , zm1_2dt  ) * ub(ji,jj,ikbu)
+            va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX(  zCdv , zm1_2dt  ) * vb(ji,jj,ikbv)
+         END DO
+      END DO
+      !
+      IF( ln_isfcav ) THEN        ! ocean cavities
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
-               ikbu = mbku(ji,jj)          ! deepest wet ocean u- & v-levels
-               ikbv = mbkv(ji,jj)
+               ikbu = miku(ji,jj)          ! first wet ocean u- & v-levels
+               ikbv = mikv(ji,jj)
                !
                ! Apply stability criteria on absolute value  : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt)
-!!gm old
-               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX(  bfrua(ji,jj) / e3u_n(ji,jj,ikbu) , zm1_2dt  ) * ub(ji,jj,ikbu)
-               va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX(  bfrva(ji,jj) / e3v_n(ji,jj,ikbv) , zm1_2dt  ) * vb(ji,jj,ikbv)
-!!gm new
-!               zCdu = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / e3u_n(ji,jj,ikbu)
-!               zCdv = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / e3v_n(ji,jj,ikbv)
-!               !
-!               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX(  zCdu , zm1_2dt  ) * ub(ji,jj,ikbu)
-!               va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX(  zCdv , zm1_2dt  ) * vb(ji,jj,ikbv)
-!!gm
-            END DO
+               zCdu = 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / e3u_n(ji,jj,ikbu)    ! NB: Cdtop masked
+               zCdv = 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / e3v_n(ji,jj,ikbv)
+               !
+               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX(  zCdu , zm1_2dt  ) * ub(ji,jj,ikbu)
+               va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX(  zCdv , zm1_2dt  ) * vb(ji,jj,ikbv)
+           END DO
          END DO
-         !
-         IF( ln_isfcav ) THEN        ! ocean cavities
-            DO jj = 2, jpjm1
-               DO ji = 2, jpim1
-                  ikbu = miku(ji,jj)          ! first wet ocean u- & v-levels
-                  ikbv = mikv(ji,jj)
-                  !
-                  ! Apply stability criteria on absolute value  : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt)
-!!gm old
-                 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX(  tfrua(ji,jj) / e3u_n(ji,jj,ikbu) , zm1_2dt  ) * ub(ji,jj,ikbu) &
-                    &             * (1.-umask(ji,jj,1))
-                 va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX(  tfrva(ji,jj) / e3v_n(ji,jj,ikbv) , zm1_2dt  ) * vb(ji,jj,ikbv) &
-                    &             * (1.-vmask(ji,jj,1))
-!!gm new
-!                  zCdu = 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / e3u_n(ji,jj,ikbu)    ! NB: Cdtop masked
-!                  zCdv = 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / e3v_n(ji,jj,ikbv)
-!                  !
-!                  ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX(  zCdu , zm1_2dt  ) * ub(ji,jj,ikbu)
-!                  va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX(  zCdv , zm1_2dt  ) * vb(ji,jj,ikbv)
-!!gm
-              END DO
-           END DO
-         END IF
-        !
-        IF( l_trddyn ) THEN      ! trends: send trends to trddyn for further diagnostics
-           ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
-           ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-           CALL trd_dyn( ztrdu(:,:,:), ztrdv(:,:,:), jpdyn_bfr, kt )
-           DEALLOCATE( ztrdu, ztrdv )
-        ENDIF
-        !                                          ! print mean trends (used for debugging)
-        IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' bfr  - Ua: ', mask1=umask,               &
-           &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-        !
-      ENDIF     ! end explicit bottom friction
+      ENDIF
+      !
+      IF( l_trddyn ) THEN      ! trends: send trends to trddyn for further diagnostics
+         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
+         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
+         CALL trd_dyn( ztrdu(:,:,:), ztrdv(:,:,:), jpdyn_bfr, kt )
+         DEALLOCATE( ztrdu, ztrdv )
+      ENDIF
+      !                                          ! print mean trends (used for debugging)
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' bfr  - Ua: ', mask1=umask,               &
+         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
       !
       IF( nn_timing == 1 )  CALL timing_stop('dyn_bfr')
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 8143)
@@ -16,4 +16,5 @@
    !!             3.7  ! 2015-11  (J. Chanut) free surface simplification
    !!              -   ! 2016-12  (G. Madec, E. Clementi) update for Stoke-Drift divergence
+   !!             4.0  ! 2017-05  (G. Madec)  drag coef. defined at t-point (zdfdrg.F90)
    !!---------------------------------------------------------------------
 
@@ -27,8 +28,6 @@
    USE dom_oce         ! ocean space and time domain
    USE sbc_oce         ! surface boundary condition: ocean
-   USE zdf_oce         ! Bottom friction coefts
-!!gm new
-   USE zdfdrg         ! vertical physics: top/bottom drag coef.
-!!gm
+   USE zdf_oce         ! vertical physics: variables
+   USE zdfdrg          ! vertical physics: top/bottom drag coef.
    USE sbcisf          ! ice shelf variable (fwfisf)
    USE sbcapr          ! surface boundary condition: atmospheric pressure
@@ -43,4 +42,11 @@
    USE updtide         ! tide potential
    USE sbcwave         ! surface wave
+   USE diatmb          ! Top,middle,bottom output
+#if defined key_agrif
+   USE agrif_opa_interp ! agrif
+#endif
+#if defined key_asminc   
+   USE asminc          ! Assimilation increment
+#endif
    !
    USE in_out_manager  ! I/O manager
@@ -50,14 +56,5 @@
    USE iom             ! IOM library
    USE restart         ! only for lrst_oce
-   USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing    
-   USE diatmb          ! Top,middle,bottom output
-#if defined key_agrif
-   USE agrif_opa_interp ! agrif
-#endif
-#if defined key_asminc   
-   USE asminc          ! Assimilation increment
-#endif
-
 
    IMPLICIT NONE
@@ -69,20 +66,24 @@
    PUBLIC ts_rst            !    "      "     "    "
 
-   INTEGER, SAVE :: icycle  ! Number of barotropic sub-steps for each internal step nn_baro <= 2.5 nn_baro
-   REAL(wp),SAVE :: rdtbt   ! Barotropic time step
-
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) ::   wgtbtp1, wgtbtp2   !: 1st & 2nd weights used in time filtering of barotropic fields
-
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  zwz          !: ff_f/h at F points
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ftnw, ftne   !: triad of coriolis parameter
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ftsw, ftse   !: (only used with een vorticity scheme)
-
    !! Time filtered arrays at baroclinic time step:
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   un_adv , vn_adv     !: Advection vel. at "now" barocl. step
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   un_adv , vn_adv   !: Advection vel. at "now" barocl. step
+
+   INTEGER , SAVE ::   icycle   ! Number of barotropic sub-steps for each internal step nn_baro <= 2.5 nn_baro
+   REAL(wp), SAVE ::   rdtbt    ! Barotropic time step
+   !
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)   ::   wgtbtp1, wgtbtp2   ! 1st & 2nd weights used in time filtering of barotropic fields
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  zwz                 ! ff_f/h at F points
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ftnw, ftne          ! triad of coriolis parameter
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ftsw, ftse          ! (only used with een vorticity scheme)
+
+   REAL(wp) ::   r1_12 = 1._wp / 12._wp   ! local ratios
+   REAL(wp) ::   r1_8  = 0.125_wp         !
+   REAL(wp) ::   r1_4  = 0.25_wp          !
+   REAL(wp) ::   r1_2  = 0.5_wp           !
 
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.5 , NEMO Consortium (2013)
+   !! NEMO/OPA 4.0 , NEMO Consortium (2017)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -140,24 +141,20 @@
       INTEGER, INTENT(in)  ::   kt   ! ocean time-step index
       !
-      LOGICAL  ::   ll_fw_start			! if true, forward integration 
-      LOGICAL  ::   ll_init			    ! if true, special startup of 2d equations
-      LOGICAL  ::   ll_tmp1, ll_tmp2            ! local logical variables used in W/D
       INTEGER  ::   ji, jj, jk, jn   		! dummy loop indices
-      INTEGER  ::   ikbu, ikbv, noffset    	! local integers
-      INTEGER  ::   iktu, iktv             	! local integers
-      REAL(wp) ::   zmdi
-      REAL(wp) ::   zraur, z1_2dt_b, z2dt_bf  	! local scalars
-      REAL(wp) ::   zx1, zy1, zx2, zy2          !   -      -
-      REAL(wp) ::   z1_12, z1_8, z1_4, z1_2 	   !   -      -
-      REAL(wp) ::   zu_spg, zv_spg              !   -      -
-      REAL(wp) ::   zhura, zhvra                !   -      -
-      REAL(wp) ::   za0, za1, za2, za3          !   -      -
-      REAL(wp) ::   zztmp                       !   -      -
-      REAL(wp), DIMENSION(jpi,jpj) :: zsshp2_e
-      REAL(wp), DIMENSION(jpi,jpj) :: zu_trd, zv_trd, zu_frc, zv_frc, zssh_frc
-      REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zhdiv
-      REAL(wp), DIMENSION(jpi,jpj) :: zhup2_e, zhvp2_e, zhust_e, zhvst_e
-      REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zsshv_a
-      REAL(wp), DIMENSION(jpi,jpj) :: zhf
+      LOGICAL  ::   ll_fw_start           ! =T : forward integration 
+      LOGICAL  ::   ll_init			      ! =T : special startup of 2d equations
+      LOGICAL  ::   ll_tmp1, ll_tmp2      ! local logical variables used in W/D
+      INTEGER  ::   ikbu, iktu, noffset   ! local integers
+      INTEGER  ::   ikbv, iktv            !   -      -
+      REAL(wp) ::   z1_2dt_b, z2dt_bf  	   ! local scalars
+      REAL(wp) ::   zx1, zx2, zu_spg, zhura  !   -      -
+      REAL(wp) ::   zy1, zy2, zv_spg, zhvra  !   -      -
+      REAL(wp) ::   za0, za1, za2, za3       !   -      -
+      REAL(wp) ::   zmdi, zztmp              !   -      -
+      REAL(wp), DIMENSION(jpi,jpj) :: zsshp2_e, zhf
+      REAL(wp), DIMENSION(jpi,jpj) :: zwx, zu_trd, zu_frc, zssh_frc
+      REAL(wp), DIMENSION(jpi,jpj) :: zwy, zv_trd, zv_frc, zhdiv
+      REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zhup2_e, zhust_e
+      REAL(wp), DIMENSION(jpi,jpj) :: zsshv_a, zhvp2_e, zhvst_e
       REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v   ! top/bottom stress at u- & v-points
       !
@@ -170,10 +167,5 @@
       !
       zmdi=1.e+20                               !  missing data indicator for masking
-      !                                         !* Local constant initialization
-      z1_12 = 1._wp / 12._wp 
-      z1_8  = 0.125_wp                                   
-      z1_4  = 0.25_wp
-      z1_2  = 0.5_wp     
-      zraur = 1._wp / rau0
+      !
       !  	                                       ! reciprocal of baroclinic time step 
       IF( kt == nit000 .AND. neuler == 0 ) THEN   ;   z2dt_bf =          rdt
@@ -210,29 +202,19 @@
       ENDIF
       !
-!!gm old/new
       IF( ln_isfcav ) THEN    ! top+bottom friction (ocean cavities)
-         zCdU_u(:,:) = bfrua(:,:) + tfrua(:,:)
-         zCdU_v(:,:) = bfrva(:,:) + tfrva(:,:)
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zCdU_u(ji,jj) = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )
+               zCdU_v(ji,jj) = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )
+            END DO
+         END DO
       ELSE                    ! bottom friction only
-         zCdU_u(:,:) = bfrua(:,:)
-         zCdU_v(:,:) = bfrva(:,:)
-      ENDIF
-!!gm new
-!      IF( ln_isfcav ) THEN    ! top+bottom friction (ocean cavities)
-!         DO jj = 2, jpjm1
-!            DO ji = fs_2, fs_jpim1   ! vector opt.
-!               zCdU_u(ji,jj) = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )
-!               zCdU_v(ji,jj) = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )
-!            END DO
-!         END DO
-!      ELSE                    ! bottom friction only
-!         DO jj = 2, jpjm1
-!            DO ji = fs_2, fs_jpim1   ! vector opt.
-!               zCdU_u(ji,jj) = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )
-!               zCdU_v(ji,jj) = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )
-!            END DO
-!         END DO
-!      ENDIF
-!!gm      
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zCdU_u(ji,jj) = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )
+               zCdU_v(ji,jj) = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )
+            END DO
+         END DO
+      ENDIF
       !
       ! Set arrays to remove/compute coriolis trend.
@@ -287,5 +269,5 @@
 !!gm 
 !!            
-              IF ( .not. ln_sco ) THEN
+              IF( .NOT.ln_sco ) THEN
   
    !!gm  agree the JC comment  : this should be done in a much clear way
@@ -338,5 +320,5 @@
       IF (.NOT.ln_bt_fw .AND.( neuler==0 .AND. kt==nit000+1 ) ) THEN
          ll_fw_start=.FALSE.
-         CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2)
+         CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 )
       ENDIF
                           
@@ -387,6 +369,6 @@
                zx2 = ( zwx(ji  ,jj) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj)
                ! energy conserving formulation for planetary vorticity term
-               zu_trd(ji,jj) = z1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
-               zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
+               zu_trd(ji,jj) =   r1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
+               zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
             END DO
          END DO
@@ -395,7 +377,7 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) &
+               zy1 =   r1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) &
                  &            + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) * r1_e1u(ji,jj)
-               zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) &
+               zx1 = - r1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) &
                  &            + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj)
                zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
@@ -407,9 +389,9 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zu_trd(ji,jj) = + z1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) &
+               zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) &
                 &                                         + ftnw(ji+1,jj) * zwy(ji+1,jj  ) &
                 &                                         + ftse(ji,jj  ) * zwy(ji  ,jj-1) &
                 &                                         + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
-               zv_trd(ji,jj) = - z1_12 * r1_e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
+               zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
                 &                                         + ftse(ji,jj+1) * zwx(ji  ,jj+1) &
                 &                                         + ftnw(ji,jj  ) * zwx(ji-1,jj  ) &
@@ -423,66 +405,65 @@
       !                                   ! ----------------------------------------------------
       IF( .NOT.ln_linssh ) THEN                 ! Variable volume : remove surface pressure gradient
-        IF( ln_wd ) THEN                        ! Calculating and applying W/D gravity filters
-           DO jj = 2, jpjm1
-              DO ji = 2, jpim1 
-                ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                &
-                     &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) .AND.            &
-                     &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji+1,jj) + ht_wd(ji+1,jj) ) &
+         IF( ln_wd ) THEN                        ! Calculating and applying W/D gravity filters
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1 
+                  ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                &
+                     &      MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) .AND.            &
+                     &      MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji+1,jj) + ht_wd(ji+1,jj) ) &
                      &                                                         > rn_wdmin1 + rn_wdmin2
-                ll_tmp2 = ( ABS( sshn(ji+1,jj)             -   sshn(ji  ,jj))  > 1.E-12 ).AND.( &
-                     &    MAX(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                &
-                     &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
-   
-                IF(ll_tmp1) THEN
-                  zcpx(ji,jj) = 1.0_wp
-                ELSE IF(ll_tmp2) THEN
-                  ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here
-                  zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) &
-                     &             / (sshn(ji+1,jj) -  sshn(ji  ,jj)) )
-                ELSE
-                  zcpx(ji,jj) = 0._wp
-                END IF
-         
-                ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                &
+                  ll_tmp2 = ( ABS( sshn(ji+1,jj)             -   sshn(ji  ,jj))  > 1.E-12 ).AND.( &
+                     &      MAX(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                &
+                     &      MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
+                     !
+                  IF(ll_tmp1) THEN
+                     zcpx(ji,jj) = 1.0_wp
+                  ELSE IF(ll_tmp2) THEN   ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here
+                     zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) &
+                        &             / (sshn(ji+1,jj) -  sshn(ji  ,jj)) )
+                  ELSE
+                     zcpx(ji,jj) = 0._wp
+                  ENDIF
+                  !
+                  ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                &
                      &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) .AND.            &
                      &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji,jj+1) + ht_wd(ji,jj+1) ) &
                      &                                                         > rn_wdmin1 + rn_wdmin2
-                ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji,jj+1))  > 1.E-12 ).AND.( &
+                  ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji,jj+1))  > 1.E-12 ).AND.( &
                      &    MAX(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                &
                      &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
-   
-                IF(ll_tmp1) THEN
-                  zcpy(ji,jj) = 1.0_wp
-                ELSE IF(ll_tmp2) THEN
-                  ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here
-                  zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) &
-                              &    / (sshn(ji,jj+1) -  sshn(ji,jj  )) )
-                ELSE
-                  zcpy(ji,jj) = 0._wp
-                END IF
-              END DO
-           END DO
- 
-           DO jj = 2, jpjm1
-              DO ji = 2, jpim1
-                 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) )   &
-                        &                        * r1_e1u(ji,jj) * zcpx(ji,jj)
-                 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji  ,jj+1) - sshn(ji  ,jj ) )   &
-                        &                        * r1_e2v(ji,jj) * zcpy(ji,jj)
-              END DO
-           END DO
-
+                     !
+                  IF(ll_tmp1) THEN
+                     zcpy(ji,jj) = 1.0_wp
+                  ELSE IF(ll_tmp2) THEN
+                     ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here
+                     zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) &
+                        &             / (sshn(ji,jj+1) -  sshn(ji,jj  )) )
+                  ELSE
+                     zcpy(ji,jj) = 0._wp
+                  ENDIF
+               END DO
+            END DO
+            !
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1
+                  zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) )   &
+                     &                            * r1_e1u(ji,jj) * zcpx(ji,jj)
+                  zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji  ,jj+1) - sshn(ji  ,jj ) )   &
+                     &                            * r1_e2v(ji,jj) * zcpy(ji,jj)
+               END DO
+            END DO
+            !
          ELSE
-
-           DO jj = 2, jpjm1
-              DO ji = fs_2, fs_jpim1   ! vector opt.
-                 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) * r1_e1u(ji,jj)
-                 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) * r1_e2v(ji,jj) 
-              END DO
-           END DO
-        ENDIF
-
-      ENDIF
-
+            !
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) * r1_e1u(ji,jj)
+                  zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) * r1_e2v(ji,jj) 
+               END DO
+            END DO
+         ENDIF
+         !
+      ENDIF
+      !
       DO jj = 2, jpjm1                          ! Remove coriolis term (and possibly spg) from barotropic trend
          DO ji = fs_2, fs_jpim1
@@ -492,6 +473,6 @@
       END DO 
       !
-      !						! Add bottom stress contribution from baroclinic velocities:      
-      IF (ln_bt_fw) THEN
+      !						! Add BOTTOM stress contribution from baroclinic velocities:      
+      IF( ln_bt_fw ) THEN
          DO jj = 2, jpjm1                          
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -518,11 +499,6 @@
          DO jj = 2, jpjm1                          
             DO ji = fs_2, fs_jpim1   ! vector opt.
-!!gm old
-               zu_frc(ji,jj) = zu_frc(ji,jj) + MAX( r1_hu_n(ji,jj) * bfrua(ji,jj) , zztmp ) * zwx(ji,jj)
-               zv_frc(ji,jj) = zv_frc(ji,jj) + MAX( r1_hv_n(ji,jj) * bfrva(ji,jj) , zztmp ) * zwy(ji,jj)
-!!gm new
-!               zu_frc(ji,jj) = zu_frc(ji,jj) + MAX( r1_hu_n(ji,jj) * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp ) * zwx(ji,jj)
-!               zv_frc(ji,jj) = zv_frc(ji,jj) + MAX( r1_hv_n(ji,jj) * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp ) * zwy(ji,jj)
-!!gm
+               zu_frc(ji,jj) = zu_frc(ji,jj) + MAX( r1_hu_n(ji,jj) * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp ) * zwx(ji,jj)
+               zv_frc(ji,jj) = zv_frc(ji,jj) + MAX( r1_hv_n(ji,jj) * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp ) * zwy(ji,jj)
             END DO
          END DO
@@ -530,62 +506,59 @@
          DO jj = 2, jpjm1                          
             DO ji = fs_2, fs_jpim1   ! vector opt.
-!!gm old
-               zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * bfrua(ji,jj) * zwx(ji,jj)
-               zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * bfrva(ji,jj) * zwy(ji,jj)
-!!gm new
-!               zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zwx(ji,jj)
-!               zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zwy(ji,jj)
-!!gm
+               zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zwx(ji,jj)
+               zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zwy(ji,jj)
             END DO
          END DO
       END IF
       !
-      !                                         ! Add top stress contribution from baroclinic velocities:      
-      IF( ln_bt_fw ) THEN
+      IF( ln_isfcav ) THEN       ! Add TOP stress contribution from baroclinic velocities:      
+         IF( ln_bt_fw ) THEN
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  iktu = miku(ji,jj)
+                  iktv = mikv(ji,jj)
+                  zwx(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj) ! NOW top baroclinic velocities
+                  zwy(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj)
+               END DO
+            END DO
+         ELSE
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  iktu = miku(ji,jj)
+                  iktv = mikv(ji,jj)
+                  zwx(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj) ! BEFORE top baroclinic velocities
+                  zwy(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj)
+               END DO
+            END DO
+         ENDIF
+         !
+         ! Note that the "unclipped" top friction parameter is used even with explicit drag
+         DO jj = 2, jpjm1              
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zwx(ji,jj)
+               zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zwy(ji,jj)
+            END DO
+         END DO
+      ENDIF
+      !       
+      IF( ln_bt_fw ) THEN                        ! Add wind forcing
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               iktu = miku(ji,jj)
-               iktv = mikv(ji,jj)
-               zwx(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj) ! NOW top baroclinic velocities
-               zwy(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj)
+               zu_frc(ji,jj) =  zu_frc(ji,jj) + r1_rau0 * utau(ji,jj) * r1_hu_n(ji,jj)
+               zv_frc(ji,jj) =  zv_frc(ji,jj) + r1_rau0 * vtau(ji,jj) * r1_hv_n(ji,jj)
             END DO
          END DO
       ELSE
+         zztmp = r1_rau0 * r1_2
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               iktu = miku(ji,jj)
-               iktv = mikv(ji,jj)
-               zwx(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj) ! BEFORE top baroclinic velocities
-               zwy(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj)
-            END DO
-         END DO
-      ENDIF
-      !
-      ! Note that the "unclipped" top friction parameter is used even with explicit drag
-      DO jj = 2, jpjm1              
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * tfrua(ji,jj) * zwx(ji,jj)
-            zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * tfrva(ji,jj) * zwy(ji,jj)
-         END DO
-      END DO
-      !       
-      IF (ln_bt_fw) THEN                        ! Add wind forcing
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zu_frc(ji,jj) =  zu_frc(ji,jj) + zraur * utau(ji,jj) * r1_hu_n(ji,jj)
-               zv_frc(ji,jj) =  zv_frc(ji,jj) + zraur * vtau(ji,jj) * r1_hv_n(ji,jj)
-            END DO
-         END DO
-      ELSE
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zu_frc(ji,jj) =  zu_frc(ji,jj) + zraur * z1_2 * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu_n(ji,jj)
-               zv_frc(ji,jj) =  zv_frc(ji,jj) + zraur * z1_2 * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv_n(ji,jj)
+               zu_frc(ji,jj) =  zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu_n(ji,jj)
+               zv_frc(ji,jj) =  zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv_n(ji,jj)
             END DO
          END DO
       ENDIF  
       !
-      IF ( ln_apr_dyn ) THEN                    ! Add atm pressure forcing
-         IF (ln_bt_fw) THEN
+      IF( ln_apr_dyn ) THEN                     ! Add atm pressure forcing
+         IF( ln_bt_fw ) THEN
             DO jj = 2, jpjm1              
                DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -597,10 +570,11 @@
             END DO
          ELSE
+            zztmp = grav * r1_2
             DO jj = 2, jpjm1              
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zu_spg =  grav * z1_2 * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj)    &
-                      &                    + ssh_ibb(ji+1,jj  ) - ssh_ibb(ji,jj)  ) * r1_e1u(ji,jj)
-                  zv_spg =  grav * z1_2 * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj)    &
-                      &                    + ssh_ibb(ji  ,jj+1) - ssh_ibb(ji,jj)  ) * r1_e2v(ji,jj)
+                  zu_spg = zztmp * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj)    &
+                      &             + ssh_ibb(ji+1,jj  ) - ssh_ibb(ji,jj)  ) * r1_e1u(ji,jj)
+                  zv_spg = zztmp * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj)    &
+                      &             + ssh_ibb(ji  ,jj+1) - ssh_ibb(ji,jj)  ) * r1_e2v(ji,jj)
                   zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg
                   zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg
@@ -613,8 +587,9 @@
       !                                         ! Surface net water flux and rivers
       IF (ln_bt_fw) THEN
-         zssh_frc(:,:) = zraur * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) )
+         zssh_frc(:,:) = r1_rau0 * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) )
       ELSE
-         zssh_frc(:,:) = zraur * z1_2 * (  emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:)   &
-                &                        + fwfisf(:,:) + fwfisf_b(:,:)                     )
+         zztmp = r1_rau0 * r1_2
+         zssh_frc(:,:) = zztmp * (  emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:)   &
+                &                 + fwfisf(:,:) + fwfisf_b(:,:)                     )
       ENDIF
       !
@@ -712,8 +687,8 @@
             DO jj = 2, jpjm1                                    ! Sea Surface Height at u- & v-points
                DO ji = 2, fs_jpim1   ! Vector opt.
-                  zwx(ji,jj) = z1_2 * ssumask(ji,jj)  * r1_e1e2u(ji,jj)     &
+                  zwx(ji,jj) = r1_2 * ssumask(ji,jj)  * r1_e1e2u(ji,jj)     &
                      &              * ( e1e2t(ji  ,jj) * zsshp2_e(ji  ,jj)  &
                      &              +   e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) )
-                  zwy(ji,jj) = z1_2 * ssvmask(ji,jj)  * r1_e1e2v(ji,jj)     &
+                  zwy(ji,jj) = r1_2 * ssvmask(ji,jj)  * r1_e1e2v(ji,jj)     &
                      &              * ( e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )  &
                      &              +   e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) )
@@ -789,8 +764,8 @@
             DO jj = 2, jpjm1
                DO ji = 2, jpim1      ! NO Vector Opt.
-                  zsshu_a(ji,jj) = z1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj)    &
+                  zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj)    &
                      &              * ( e1e2t(ji  ,jj  )  * ssha_e(ji  ,jj  ) &
                      &              +   e1e2t(ji+1,jj  )  * ssha_e(ji+1,jj  ) )
-                  zsshv_a(ji,jj) = z1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj)    &
+                  zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj)    &
                      &              * ( e1e2t(ji  ,jj  )  * ssha_e(ji  ,jj  ) &
                      &              +   e1e2t(ji  ,jj+1)  * ssha_e(ji  ,jj+1) )
@@ -868,8 +843,8 @@
             DO jj = 2, jpjm1                            
                DO ji = 2, jpim1
-                  zx1 = z1_2 * ssumask(ji  ,jj) *  r1_e1e2u(ji  ,jj)    &
+                  zx1 = r1_2 * ssumask(ji  ,jj) *  r1_e1e2u(ji  ,jj)    &
                      &      * ( e1e2t(ji  ,jj  ) * zsshp2_e(ji  ,jj)    &
                      &      +   e1e2t(ji+1,jj  ) * zsshp2_e(ji+1,jj  ) )
-                  zy1 = z1_2 * ssvmask(ji  ,jj) *  r1_e1e2v(ji  ,jj  )  &
+                  zy1 = r1_2 * ssvmask(ji  ,jj) *  r1_e1e2v(ji  ,jj  )  &
                      &       * ( e1e2t(ji ,jj  ) * zsshp2_e(ji  ,jj  )  &
                      &       +   e1e2t(ji ,jj+1) * zsshp2_e(ji  ,jj+1) )
@@ -895,6 +870,6 @@
                   zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj)
                   zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj)
-                  zu_trd(ji,jj) = z1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
-                  zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
+                  zu_trd(ji,jj) = r1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
+                  zv_trd(ji,jj) =-r1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
                END DO
             END DO
@@ -903,7 +878,7 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) &
+                  zy1 =   r1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) &
                    &             + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) * r1_e1u(ji,jj)
-                  zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) &
+                  zx1 = - r1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) &
                    &             + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) * r1_e2v(ji,jj)
                   zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
@@ -915,9 +890,9 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zu_trd(ji,jj) = + z1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) &
+                  zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) &
                      &                                       + ftnw(ji+1,jj) * zwy(ji+1,jj  ) &
                      &                                       + ftse(ji,jj  ) * zwy(ji  ,jj-1) & 
                      &                                       + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
-                  zv_trd(ji,jj) = - z1_12 * r1_e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 
+                  zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 
                      &                                       + ftse(ji,jj+1) * zwx(ji  ,jj+1) &
                      &                                       + ftnw(ji,jj  ) * zwx(ji-1,jj  ) & 
@@ -1082,6 +1057,6 @@
          vn_adv(:,:) = zwy(:,:) * r1_hv_n(:,:)
       ELSE
-         un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zwx(:,:) ) * r1_hu_n(:,:)
-         vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zwy(:,:) ) * r1_hv_n(:,:)
+         un_adv(:,:) = r1_2 * ( ub2_b(:,:) + zwx(:,:) ) * r1_hu_n(:,:)
+         vn_adv(:,:) = r1_2 * ( vb2_b(:,:) + zwy(:,:) ) * r1_hv_n(:,:)
       END IF
 
@@ -1101,8 +1076,8 @@
          DO jj = 1, jpjm1
             DO ji = 1, jpim1      ! NO Vector Opt.
-               zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e1e2u(ji,jj) &
+               zsshu_a(ji,jj) = r1_2 * umask(ji,jj,1)  * r1_e1e2u(ji,jj) &
                   &              * ( e1e2t(ji  ,jj) * ssha(ji  ,jj)    &
                   &              +   e1e2t(ji+1,jj) * ssha(ji+1,jj) )
-               zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e1e2v(ji,jj) &
+               zsshv_a(ji,jj) = r1_2 * vmask(ji,jj,1)  * r1_e1e2v(ji,jj) &
                   &              * ( e1e2t(ji,jj  ) * ssha(ji,jj  )    &
                   &              +   e1e2t(ji,jj+1) * ssha(ji,jj+1) )
@@ -1299,10 +1274,8 @@
       INTEGER  ::   ji ,jj              ! dummy loop indices
       REAL(wp) ::   zxr2, zyr2, zcmax   ! local scalar
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zcu
+      REAL(wp), DIMENSION(jpi,jpj) ::   zcu
       !!----------------------------------------------------------------------
       !
       ! Max courant number for ext. grav. waves
-      !
-      CALL wrk_alloc( jpi,jpj,   zcu )
       !
       DO jj = 1, jpj
@@ -1371,6 +1344,4 @@
       ENDIF
       !
-      CALL wrk_dealloc( jpi,jpj,   zcu )
-      !
    END SUBROUTINE dyn_spg_ts_init
 
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 8143)
@@ -9,4 +9,5 @@
    !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
    !!            3.4  !  2012-01  (H. Liu) Semi-implicit bottom friction
+   !!            4.0  !  2017-05  (G. Madec)  drag coef. defined at t-point (zdfdrg.F90)
    !!----------------------------------------------------------------------
 
@@ -22,13 +23,8 @@
    USE dynadv   , ONLY: ln_dynadv_vec ! Momentum advection form
    USE zdf_oce        ! ocean vertical physics
-!!gm new
    USE zdfdrg         ! vertical physics: top/bottom drag coef.
-!!gm old
-   USE zdfbfr         ! Bottom friction setup
-!!gm
    !
    USE in_out_manager ! I/O manager
    USE lib_mpp        ! MPP library
-   USE wrk_nemo       ! Memory Allocation
    USE timing         ! Timing
 
@@ -65,5 +61,5 @@
       !!      with the following surface/top/bottom boundary condition:
       !!      surface: wind stress input (averaged over kt-1/2 & kt+1/2)
-      !!      top & bottom : top stress (iceshelf-ocean) & bottom stress (cf zdfbfr.F)
+      !!      top & bottom : top stress (iceshelf-ocean) & bottom stress (cf zdfdrg.F90)
       !!
       !! ** Action :   (ua,va) after velocity 
@@ -76,10 +72,8 @@
       REAL(wp) ::   zzwi, ze3ua   ! local scalars
       REAL(wp) ::   zzws, ze3va   !   -      -
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwi, zwd, zws
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zwi, zwd, zws
       !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('dyn_zdf_imp')
-      !
-      CALL wrk_alloc( jpi,jpj,jpk, zwi, zwd, zws ) 
       !
       IF( kt == nit000 ) THEN
@@ -115,16 +109,11 @@
       ! column vector of the tri-diagonal matrix equation
       !
-      IF( ln_bfrimp ) THEN
+      IF( ln_drgimp ) THEN
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
                ikbu = mbku(ji,jj)       ! ocean bottom level at u- and v-points 
                ikbv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
-!!gm old
-               avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * e3uw_n(ji,jj,ikbu+1)
-               avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * e3vw_n(ji,jj,ikbv+1)
-!!gm new
-!               avmu(ji,jj,ikbu+1) = -0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * e3uw_n(ji,jj,ikbu+1)
-!               avmv(ji,jj,ikbv+1) = -0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * e3vw_n(ji,jj,ikbv+1)
-!!gm
+               avmu(ji,jj,ikbu+1) = -0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * e3uw_n(ji,jj,ikbu+1)
+               avmv(ji,jj,ikbv+1) = -0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * e3vw_n(ji,jj,ikbv+1)
             END DO
          END DO
@@ -134,12 +123,7 @@
                   ikbu = miku(ji,jj)       ! ocean top level at u- and v-points 
                   ikbv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
-!!gm old
-                  IF( ikbu >= 2 )   avmu(ji,jj,ikbu) = -tfrua(ji,jj) * e3uw_n(ji,jj,ikbu)
-                  IF( ikbv >= 2 )   avmv(ji,jj,ikbv) = -tfrva(ji,jj) * e3vw_n(ji,jj,ikbv)
-!!gm new
                   ! top Cd is masked (=0 outside cavities) no need of test on mik>=2
-!                  avmu(ji,jj,ikbu) = -0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * e3uw_n(ji,jj,ikbu)
-!                  avmv(ji,jj,ikbv) = -0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * e3vw_n(ji,jj,ikbv)
-!!gm
+                  avmu(ji,jj,ikbu) = -0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * e3uw_n(ji,jj,ikbu)
+                  avmv(ji,jj,ikbv) = -0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * e3vw_n(ji,jj,ikbv)
                END DO
             END DO
@@ -152,5 +136,5 @@
       !            not lead to the effective stress seen over the whole barotropic loop. 
       ! G. Madec : in linear free surface, e3u_a = e3u_n = e3u_0, so systematic use of e3u_a
-      IF( ln_bfrimp .AND. ln_dynspg_ts ) THEN
+      IF( ln_drgimp .AND. ln_dynspg_ts ) THEN
          DO jk = 1, jpkm1        ! remove barotropic velocities
             ua(:,:,jk) = ( ua(:,:,jk) - ua_b(:,:) ) * umask(:,:,jk)
@@ -163,11 +147,6 @@
                ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,ikbu) + r_vvl * e3u_a(ji,jj,ikbu)
                ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikbv) + r_vvl * e3v_a(ji,jj,ikbv)
-!!gm old
-               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * bfrua(ji,jj) * ua_b(ji,jj) / ze3ua
-               va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * bfrva(ji,jj) * va_b(ji,jj) / ze3va
-!!gm new
-!               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua
-!               va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va
-!!gm
+               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua
+               va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va
             END DO
          END DO
@@ -179,11 +158,6 @@
                   ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,ikbu) + r_vvl * e3u_a(ji,jj,ikbu)
                   ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikbv) + r_vvl * e3v_a(ji,jj,ikbv)
-!!gm old
-                  ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * tfrua(ji,jj) * ua_b(ji,jj) / ze3ua
-                  va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * tfrva(ji,jj) * va_b(ji,jj) / ze3va
-!!gm new
                   ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua
                   va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va
-!!gm
                END DO
             END DO
@@ -342,30 +316,4 @@
          END DO
       END DO
-      
-      ! J. Chanut: Lines below are useless ?
-      !! restore bottom layer avmu(v) 
-      !!gm  I almost sure it is !!!!
-      IF( ln_bfrimp ) THEN
-        DO jj = 2, jpjm1
-           DO ji = 2, jpim1
-              ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points 
-              ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
-              avmu(ji,jj,ikbu+1) = 0._wp
-              avmv(ji,jj,ikbv+1) = 0._wp
-           END DO
-        END DO
-        IF (ln_isfcav) THEN
-           DO jj = 2, jpjm1
-              DO ji = 2, jpim1
-                 ikbu = miku(ji,jj)         ! ocean top level at u- and v-points 
-                 ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
-                 IF( ikbu > 1 )   avmu(ji,jj,ikbu) = 0._wp
-                 IF( ikbv > 1 )   avmv(ji,jj,ikbv) = 0._wp
-              END DO
-           END DO
-        ENDIF
-      ENDIF
-      !
-      CALL wrk_dealloc( jpi,jpj,jpk,   zwi, zwd, zws) 
       !
       IF( nn_timing == 1 )   CALL timing_stop('dyn_zdf_imp')
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 8143)
@@ -5,26 +5,26 @@
    !!                   shelf
    !!======================================================================
-   !! History :  3.2   !  2011-02  (C.Harris  ) Original code isf cav
-   !!            X.X   !  2006-02  (C. Wang   ) Original code bg03
-   !!            3.4   !  2013-03  (P. Mathiot) Merging + parametrization
+   !! History :  3.2  !  2011-02  (C.Harris  ) Original code isf cav
+   !!            X.X  !  2006-02  (C. Wang   ) Original code bg03
+   !!            3.4  !  2013-03  (P. Mathiot) Merging + parametrization
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
-   !!   sbc_isf        : update sbc under ice shelf
+   !!   sbc_isf       : update sbc under ice shelf
    !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers
-   USE dom_oce         ! ocean space and time domain
-   USE phycst          ! physical constants
-   USE eosbn2          ! equation of state
-   USE sbc_oce         ! surface boundary condition: ocean fields
-   USE zdfbfr          !
+   USE oce            ! ocean dynamics and tracers
+   USE dom_oce        ! ocean space and time domain
+   USE phycst         ! physical constants
+   USE eosbn2         ! equation of state
+   USE sbc_oce        ! surface boundary condition: ocean fields
+   USE zdfdrg         ! vertical physics: top/bottom drag coef.
    !
-   USE in_out_manager  ! I/O manager
-   USE iom             ! I/O manager library
-   USE fldread         ! read input field at current time step
-   USE lbclnk          !
-   USE wrk_nemo        ! Memory allocation
-   USE timing          ! Timing
-   USE lib_fortran     ! glob_sum
+   USE in_out_manager ! I/O manager
+   USE iom            ! I/O manager library
+   USE fldread        ! read input field at current time step
+   USE lbclnk         !
+   USE wrk_nemo       ! Memory allocation
+   USE timing         ! Timing
+   USE lib_fortran    ! glob_sum
 
    IMPLICIT NONE
@@ -77,5 +77,5 @@
 CONTAINS
  
-  SUBROUTINE sbc_isf(kt)
+  SUBROUTINE sbc_isf( kt )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE sbc_isf  ***
@@ -94,5 +94,5 @@
       INTEGER               :: ji, jj, jk           ! loop index
       INTEGER               :: ikt, ikb             ! loop index
-      REAL(wp), DIMENSION (:,:), POINTER :: zt_frz, zdep ! freezing temperature (zt_frz) at depth (zdep) 
+      REAL(wp), DIMENSION(jpi,jpj) ::   zt_frz, zdep   ! freezing temperature (zt_frz) at depth (zdep) 
       REAL(wp), DIMENSION(:,:,:), POINTER :: zfwfisf3d, zqhcisf3d, zqlatisf3d
       REAL(wp), DIMENSION(:,:  ), POINTER :: zqhcisf2d
@@ -100,6 +100,4 @@
       !
       IF( MOD( kt-1, nn_fsbc) == 0 ) THEN
-         ! allocation
-         CALL wrk_alloc( jpi,jpj, zt_frz, zdep  )
 
          ! compute salt and heat flux
@@ -204,6 +202,4 @@
             CALL wrk_dealloc( jpi,jpj,     zqhcisf2d                        )
           END IF
-          ! deallocation
-          CALL wrk_dealloc( jpi,jpj, zt_frz, zdep  )
           !
         END IF
@@ -254,4 +250,5 @@
   END FUNCTION
 
+
   SUBROUTINE sbc_isf_init
       !!---------------------------------------------------------------------
@@ -289,8 +286,6 @@
 
       IF ( lwp ) WRITE(numout,*)
-      IF ( lwp ) WRITE(numout,*) 'sbc_isf: heat flux of the ice shelf'
-      IF ( lwp ) WRITE(numout,*) '~~~~~~~~~'
-      IF ( lwp ) WRITE(numout,*) 'sbcisf :' 
-      IF ( lwp ) WRITE(numout,*) '~~~~~~~~'
+      IF ( lwp ) WRITE(numout,*) 'sbc_isf_init : heat flux of the ice shelf'
+      IF ( lwp ) WRITE(numout,*) '~~~~~~~~~~~'
       IF ( lwp ) WRITE(numout,*) '        nn_isf      = ', nn_isf
       IF ( lwp ) WRITE(numout,*) '        nn_isfblk   = ', nn_isfblk
@@ -299,5 +294,5 @@
       IF ( lwp ) WRITE(numout,*) '        rn_gammat0  = ', rn_gammat0  
       IF ( lwp ) WRITE(numout,*) '        rn_gammas0  = ', rn_gammas0  
-      IF ( lwp ) WRITE(numout,*) '        rn_tfri2    = ', rn_tfri2 
+      IF ( lwp ) WRITE(numout,*) '        rn_Cd0      = ', r_Cdmin_top 
       !
       ! Allocate public variable
@@ -305,6 +300,6 @@
       !
       ! initialisation
-      qisf(:,:)        = 0._wp  ; fwfisf  (:,:) = 0._wp
-      risf_tsc(:,:,:)  = 0._wp  ; fwfisf_b(:,:) = 0._wp
+      qisf    (:,:)    = 0._wp   ;   fwfisf  (:,:) = 0._wp
+      risf_tsc(:,:,:)  = 0._wp   ;   fwfisf_b(:,:) = 0._wp
       !
       ! define isf tbl tickness, top and bottom indice
@@ -312,5 +307,5 @@
       CASE ( 1 ) 
          rhisf_tbl(:,:) = rn_hisf_tbl
-         misfkt(:,:)    = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
+         misfkt   (:,:) = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
 
       CASE ( 2 , 3 )
@@ -346,4 +341,5 @@
             DO jj = 1, jpj
                 ik = 2
+!!gm potential bug: use gdepw_0 not _n
                 DO WHILE ( ik <= mbkt(ji,jj) .AND. gdepw_n(ji,jj,ik) < rzisf_tbl(ji,jj) ) ;  ik = ik + 1 ;  END DO
                 misfkt(ji,jj) = ik-1
@@ -354,5 +350,5 @@
          ! as in nn_isf == 1
          rhisf_tbl(:,:) = rn_hisf_tbl
-         misfkt(:,:)    = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
+         misfkt   (:,:) = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
          
          ! load variable used in fldread (use for temporal interpolation of isf fwf forcing)
@@ -377,5 +373,5 @@
             ! determine the deepest level influenced by the boundary layer
             DO jk = ikt+1, mbkt(ji,jj)
-               IF ( (SUM(e3t_n(ji,jj,ikt:jk-1)) < rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
+               IF( (SUM(e3t_n(ji,jj,ikt:jk-1)) < rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) )   ikb = jk
             END DO
             rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness.
@@ -390,4 +386,5 @@
   END SUBROUTINE sbc_isf_init
 
+
   SUBROUTINE sbc_isf_bg03(kt)
       !!---------------------------------------------------------------------
@@ -402,6 +399,5 @@
       !!         interaction for climate models", Ocean Modelling 5(2003) 157-170.
       !!         (hereafter BG)
-      !! History :
-      !!         06-02  (C. Wang) Original code
+      !! History :  06-02  (C. Wang) Original code
       !!----------------------------------------------------------------------
       INTEGER, INTENT ( in ) :: kt
@@ -415,5 +411,5 @@
       !!----------------------------------------------------------------------
 
-      IF ( nn_timing == 1 ) CALL timing_start('sbc_isf_bg03')
+      IF( nn_timing == 1 )   CALL timing_start('sbc_isf_bg03')
       !
       DO ji = 1, jpi
@@ -441,12 +437,13 @@
                !add to salinity trend
             ELSE
-               qisf(ji,jj) = 0._wp ; fwfisf(ji,jj) = 0._wp
+               qisf(ji,jj) = 0._wp   ;   fwfisf(ji,jj) = 0._wp
             END IF
          END DO
       END DO
       !
-      IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_bg03')
+      IF( nn_timing == 1 )   CALL timing_stop('sbc_isf_bg03')
       !
   END SUBROUTINE sbc_isf_bg03
+
 
   SUBROUTINE sbc_isf_cav( kt )
@@ -463,8 +460,9 @@
       !!                emp, emps  : update freshwater flux below ice shelf
       !!---------------------------------------------------------------------
-      INTEGER, INTENT(in)          ::   kt         ! ocean time step
+      INTEGER, INTENT(in) ::   kt   ! ocean time step
       !
       INTEGER  ::   ji, jj     ! dummy loop indices
       INTEGER  ::   nit
+      LOGICAL  ::   lit
       REAL(wp) ::   zlamb1, zlamb2, zlamb3
       REAL(wp) ::   zeps1,zeps2,zeps3,zeps4,zeps6,zeps7
@@ -472,8 +470,7 @@
       REAL(wp) ::   zeps = 1.e-20_wp        
       REAL(wp) ::   zerr
-      REAL(wp), DIMENSION(:,:), POINTER ::   zfrz
-      REAL(wp), DIMENSION(:,:), POINTER ::   zgammat, zgammas 
-      REAL(wp), DIMENSION(:,:), POINTER ::   zfwflx, zhtflx, zhtflx_b
-      LOGICAL  ::   lit
+      REAL(wp), DIMENSION(jpi,jpj) ::   zfrz
+      REAL(wp), DIMENSION(jpi,jpj) ::   zgammat, zgammas 
+      REAL(wp), DIMENSION(jpi,jpj) ::   zfwflx, zhtflx, zhtflx_b
       !!---------------------------------------------------------------------
       ! coeficient for linearisation of potential tfreez
@@ -484,7 +481,4 @@
       IF( nn_timing == 1 )  CALL timing_start('sbc_isf_cav')
       !
-      CALL wrk_alloc( jpi,jpj, zfrz  , zgammat, zgammas  )
-      CALL wrk_alloc( jpi,jpj, zfwflx, zhtflx , zhtflx_b )
-
       ! initialisation
       zgammat(:,:) = rn_gammat0 ; zgammas (:,:) = rn_gammas0
@@ -578,7 +572,4 @@
       CALL iom_put('isfgammas', zgammas)
       ! 
-      CALL wrk_dealloc( jpi,jpj, zfrz  , zgammat, zgammas  )
-      CALL wrk_dealloc( jpi,jpj, zfwflx, zhtflx , zhtflx_b )
-      !
       IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_cav')
       !
@@ -600,5 +591,4 @@
       INTEGER  :: ikt                        
       INTEGER  :: ji, jj                     ! loop index
-      REAL(wp), DIMENSION(:,:), POINTER :: zustar           ! U, V at T point and friction velocity
       REAL(wp) :: zdku, zdkv                 ! U, V shear 
       REAL(wp) :: zPr, zSc, zRc              ! Prandtl, Scmidth and Richardson number 
@@ -614,6 +604,6 @@
       REAL(wp), PARAMETER :: znu   = 1.95e-6_wp ! kinamatic viscosity of sea water (m2.s-1)
       REAL(wp), DIMENSION(2) :: zts, zab
+      REAL(wp), DIMENSION(jpi,jpj) :: zustar   ! U, V at T point and friction velocity
       !!---------------------------------------------------------------------
-      CALL wrk_alloc( jpi,jpj, zustar )
       !
       SELECT CASE ( nn_gammablk )
@@ -626,7 +616,9 @@
          !! Jenkins et al., 2010, JPO, p2298-2312
          !! Adopted by Asay-Davis et al. (2015)
-
-         !! compute ustar (eq. 24)
-         zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) )
+!!gm  I don't understand the u* expression in those papers... (see for example zdfglf module)
+!!    for me ustar= Cd0 * |U|  not  (Cd0)^1/2 * |U| ....  which is what you can find in Jenkins et al.
+
+         !! compute ustar (eq. 24)           !! NB: here r_Cdmin_top = rn_Cd0 read in namdrg_top namelist)
+         zustar(:,:) = SQRT( r_Cdmin_top * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + r_ke0_top) )
 
          !! Compute gammats
@@ -638,5 +630,5 @@
          !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO)
          !! compute ustar
-         zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) )
+         zustar(:,:) = SQRT( r_Cdmin_top * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + r_ke0_top) )
 
          !! compute Pr and Sc number (can be improved)
@@ -649,13 +641,15 @@
 
          !! compute gamma
-         DO ji=2,jpi
-            DO jj=2,jpj
+         DO ji = 2, jpi
+            DO jj = 2, jpj
                ikt = mikt(ji,jj)
 
-               IF (zustar(ji,jj) == 0._wp) THEN           ! only for kt = 1 I think
+               IF( zustar(ji,jj) == 0._wp ) THEN           ! only for kt = 1 I think
                   pgt = rn_gammat0
                   pgs = rn_gammas0
                ELSE
                   !! compute Rc number (as done in zdfric.F90)
+!!gm better to do it like in the new zdfric.F90   i.e. avm weighted Ri computation
+!!gm moreover, use Max(rn2,0) to take care of static instabilities....
                   zcoef = 0.5_wp / e3w_n(ji,jj,ikt)
                   !                                            ! shear of horizontal velocity
@@ -703,8 +697,8 @@
          CALL lbc_lnk(pgs(:,:),'T',1.)
       END SELECT
-      CALL wrk_dealloc( jpi,jpj, zustar )
       !
    END SUBROUTINE sbc_isf_gammats
 
+
    SUBROUTINE sbc_isf_tbl( pvarin, pvarout, cd_ptin )
       !!----------------------------------------------------------------------
@@ -714,13 +708,12 @@
       !!
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:,:), INTENT( in  ) :: pvarin
-      REAL(wp), DIMENSION(:,:)  , INTENT( out ) :: pvarout
-      CHARACTER(len=1),           INTENT( in  ) :: cd_ptin ! point of variable in/out
-      !
-      REAL(wp) :: ze3, zhk
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) :: pvarin
+      REAL(wp), DIMENSION(:,:)  , INTENT(  out) :: pvarout
+      CHARACTER(len=1),           INTENT(in   ) :: cd_ptin ! point of variable in/out
+      !
+      INTEGER ::   ji, jj, jk                ! loop index
+      INTEGER ::   ikt, ikb                    ! top and bottom index of the tbl
+      REAL(wp) ::   ze3, zhk
       REAL(wp), DIMENSION(:,:), POINTER :: zhisf_tbl ! thickness of the tbl
-
-      INTEGER :: ji, jj, jk                  ! loop index
-      INTEGER :: ikt, ikb                    ! top and bottom index of the tbl
       !!----------------------------------------------------------------------
       ! allocation
@@ -736,5 +729,5 @@
                ikt = miku(ji,jj) ; ikb = miku(ji,jj)
                ! thickness of boundary layer at least the top level thickness
-               zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3u_n(ji,jj,ikt))
+               zhisf_tbl(ji,jj) = MAX( rhisf_tbl_0(ji,jj) , e3u_n(ji,jj,ikt) )
 
                ! determine the deepest level influenced by the boundary layer
@@ -755,6 +748,7 @@
             END DO
          END DO
-         DO jj = 2,jpj
-            DO ji = 2,jpi
+         DO jj = 2, jpj
+            DO ji = 2, jpi
+!!gm a wet-point only average should be used here !!!
                pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji-1,jj))
             END DO
@@ -786,6 +780,7 @@
             END DO
          END DO
-         DO jj = 2,jpj
-            DO ji = 2,jpi
+         DO jj = 2, jpj
+            DO ji = 2, jpi
+!!gm a wet-point only average should be used here !!!
                pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji,jj-1))
             END DO
@@ -882,4 +877,5 @@
       !
    END SUBROUTINE sbc_isf_div
+
    !!======================================================================
 END MODULE sbcisf
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/TRD/trd_oce.F90
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/TRD/trd_oce.F90	(revision 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/TRD/trd_oce.F90	(revision 8143)
@@ -72,5 +72,5 @@
    INTEGER, PUBLIC, PARAMETER ::   jpdyn_atf  = 10     !: Asselin time filter
    INTEGER, PUBLIC, PARAMETER ::   jpdyn_tau  = 11     !: surface stress
-   INTEGER, PUBLIC, PARAMETER ::   jpdyn_bfri = 12     !: implicit bottom friction (ln_bfrimp=.TRUE.)
+   INTEGER, PUBLIC, PARAMETER ::   jpdyn_bfri = 12     !: implicit bottom friction (ln_drgimp=.TRUE.)
    INTEGER, PUBLIC, PARAMETER ::   jpdyn_ken  = 13     !: use for calculation of KE
    !
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/TRD/trddyn.F90
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/TRD/trddyn.F90	(revision 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/TRD/trddyn.F90	(revision 8143)
@@ -15,28 +15,28 @@
    USE oce            ! ocean dynamics and tracers variables
    USE dom_oce        ! ocean space and time domain variables
-   USE zdf_oce        ! ocean vertical physics variables
+   USE phycst         ! physical constants
+   USE sbc_oce        ! surface boundary condition: ocean
+   USE zdf_oce        ! ocean vertical physics: variables
+   USE zdfdrg         ! ocean vertical physics: bottom friction
    USE trd_oce        ! trends: ocean variables
-   USE zdfbfr         ! bottom friction
-   USE sbc_oce        ! surface boundary condition: ocean
-   USE phycst         ! physical constants
    USE trdken         ! trends: Kinetic ENergy 
    USE trdglo         ! trends: global domain averaged
    USE trdvor         ! trends: vertical averaged vorticity 
    USE trdmxl         ! trends: mixed layer averaged 
+   !
    USE in_out_manager ! I/O manager
    USE lbclnk         ! lateral boundary condition 
    USE iom            ! I/O manager library
    USE lib_mpp        ! MPP library
-   USE wrk_nemo       ! Memory allocation
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC trd_dyn        ! called by all dynXX modules
+   PUBLIC trd_dyn        ! called by all dynXXX modules
 
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 4.0 , NEMO Consortium (2017)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -103,6 +103,6 @@
       INTEGER ::   ji, jj, jk   ! dummy loop indices
       INTEGER ::   ikbu, ikbv   ! local integers
-      REAL(wp), POINTER, DIMENSION(:,:)   ::   z2dx, z2dy   ! 2D workspace 
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z3dx, z3dy   ! 3D workspace 
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   z2dx, z2dy   ! 2D workspace 
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   z3dx, z3dy   ! 3D workspace 
       !!----------------------------------------------------------------------
       !
@@ -118,5 +118,5 @@
       CASE( jpdyn_keg )   ;   CALL iom_put( "utrd_keg", putrd )    ! Kinetic Energy gradient (or had)
                               CALL iom_put( "vtrd_keg", pvtrd )
-                              CALL wrk_alloc( jpi, jpj, jpk, z3dx, z3dy )
+                              ALLOCATE( z3dx(jpi,jpj,jpk) , z3dy(jpi,jpj,jpk) )
                               z3dx(:,:,:) = 0._wp                  ! U.dxU & V.dyV (approximation)
                               z3dy(:,:,:) = 0._wp
@@ -133,25 +133,22 @@
                               CALL iom_put( "utrd_udx", z3dx  )
                               CALL iom_put( "vtrd_vdy", z3dy  )
-                              CALL wrk_dealloc( jpi, jpj, jpk, z3dx, z3dy )
-      CASE( jpdyn_zad )   ;   CALL iom_put( "utrd_zad", putrd )    ! vertical   advection
+                              DEALLOCATE( z3dx , z3dy )
+      CASE( jpdyn_zad )   ;   CALL iom_put( "utrd_zad", putrd )    ! vertical advection
                               CALL iom_put( "vtrd_zad", pvtrd )
-      CASE( jpdyn_ldf )   ;   CALL iom_put( "utrd_ldf", putrd )    ! lateral diffusion
+      CASE( jpdyn_ldf )   ;   CALL iom_put( "utrd_ldf", putrd )    ! lateral  diffusion
                               CALL iom_put( "vtrd_ldf", pvtrd )
       CASE( jpdyn_zdf )   ;   CALL iom_put( "utrd_zdf", putrd )    ! vertical diffusion 
                               CALL iom_put( "vtrd_zdf", pvtrd )
+                              !
                               !                                    ! wind stress trends
-                              CALL wrk_alloc( jpi, jpj, z2dx, z2dy )
+                              ALLOCATE( z2dx(jpi,jpj) , z2dy(jpi,jpj) )
                               z2dx(:,:) = ( utau_b(:,:) + utau(:,:) ) / ( e3u_n(:,:,1) * rau0 )
                               z2dy(:,:) = ( vtau_b(:,:) + vtau(:,:) ) / ( e3v_n(:,:,1) * rau0 )
                               CALL iom_put( "utrd_tau", z2dx )
                               CALL iom_put( "vtrd_tau", z2dy )
-                              CALL wrk_dealloc( jpi, jpj, z2dx, z2dy )
-      CASE( jpdyn_bfr )       ! called if ln_bfrimp=T
-                              CALL iom_put( "utrd_bfr", putrd )    ! bottom friction (explicit case)
-                              CALL iom_put( "vtrd_bfr", pvtrd )
-      CASE( jpdyn_atf )   ;   CALL iom_put( "utrd_atf", putrd )        ! asselin filter trends 
-                              CALL iom_put( "vtrd_atf", pvtrd )
-      CASE( jpdyn_bfri )  ;   IF( ln_bfrimp ) THEN                     ! bottom friction (implicit case)
-                              	CALL wrk_alloc( jpi, jpj, jpk, z3dx, z3dy )
+                              DEALLOCATE( z2dx , z2dy )
+                              !                                    ! bottom stress tends (implicit case)
+                              IF( ln_drgimp ) THEN
+                                 ALLOCATE( z3dx(jpi,jpj,jpk) , z3dy(jpi,jpj,jpk) )
                               	z3dx(:,:,:) = 0._wp   ;   z3dy(:,:,:) = 0._wp  ! after velocity known (now filed at this stage)
                               	DO jk = 1, jpkm1
@@ -160,14 +157,19 @@
                               	         ikbu = mbku(ji,jj)          ! deepest ocean u- & v-levels
                                           ikbv = mbkv(ji,jj)
-                                          z3dx(ji,jj,jk) = bfrua(ji,jj) * un(ji,jj,ikbu) / e3u_n(ji,jj,ikbu)
-                                          z3dy(ji,jj,jk) = bfrva(ji,jj) * vn(ji,jj,ikbv) / e3v_n(ji,jj,ikbv)
+                                          z3dx(ji,jj,jk) = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )*un(ji,jj,ikbu)/e3u_n(ji,jj,ikbu)
+                                          z3dy(ji,jj,jk) = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )*vn(ji,jj,ikbv)/e3v_n(ji,jj,ikbv)
                               	      END DO
                               	   END DO
                               	END DO
-                              	CALL lbc_lnk( z3dx, 'U', -1. ) ; CALL lbc_lnk( z3dy, 'V', -1. )
-                              	CALL iom_put( "utrd_bfri", z3dx )
-                              	CALL iom_put( "vtrd_bfri", z3dy )
-                              	CALL wrk_dealloc( jpi, jpj, jpk, z3dx, z3dy )
-         					      ENDIF
+                              	CALL lbc_lnk( z3dx, 'U', -1. )   ;   CALL lbc_lnk( z3dy, 'V', -1. )
+                              	CALL iom_put( "utrd_bfr", z3dx )
+                              	CALL iom_put( "vtrd_bfr", z3dy )
+                                 DEALLOCATE( z3dx , z3dy )
+                              ENDIF
+      CASE( jpdyn_bfr )       ! called if ln_drgimp=F
+                              CALL iom_put( "utrd_bfr", putrd )    ! bottom friction (explicit case)
+                              CALL iom_put( "vtrd_bfr", pvtrd )
+      CASE( jpdyn_atf )   ;   CALL iom_put( "utrd_atf", putrd )        ! asselin filter trends 
+                              CALL iom_put( "vtrd_atf", pvtrd )
       END SELECT
       !
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/TRD/trdglo.F90
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/TRD/trdglo.F90	(revision 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/TRD/trdglo.F90	(revision 8143)
@@ -9,26 +9,25 @@
 
    !!----------------------------------------------------------------------
-   !!   trd_glo      : domain averaged budget of trends (including kinetic energy and T^2 trends)
-   !!   glo_dyn_wri  : print dynamic trends in ocean.output file
-   !!   glo_tra_wri  : print global T & T^2 trends in ocean.output file
-   !!   trd_glo_init : initialization step
+   !!   trd_glo       : domain averaged budget of trends (including kinetic energy and T^2 trends)
+   !!   glo_dyn_wri   : print dynamic trends in ocean.output file
+   !!   glo_tra_wri   : print global T & T^2 trends in ocean.output file
+   !!   trd_glo_init  : initialization step
    !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers variables
-   USE dom_oce         ! ocean space and time domain variables
-   USE sbc_oce         ! surface boundary condition: ocean
-   USE trd_oce         ! trends: ocean variables
-   USE phycst          ! physical constants
-   USE ldftra          ! lateral diffusion: eddy diffusivity & EIV coeff.
-   USE ldfdyn          ! ocean dynamics: lateral physics
-   USE zdf_oce         ! ocean vertical physics
-   USE zdfbfr          ! bottom friction
-   USE zdfddm          ! ocean vertical physics: double diffusion
-   USE eosbn2          ! equation of state
-   USE phycst          ! physical constants
+   USE oce            ! ocean dynamics and tracers variables
+   USE dom_oce        ! ocean space and time domain variables
+   USE sbc_oce        ! surface boundary condition: ocean
+   USE trd_oce        ! trends: ocean variables
+   USE phycst         ! physical constants
+   USE ldftra         ! lateral diffusion: eddy diffusivity & EIV coeff.
+   USE ldfdyn         ! ocean dynamics: lateral physics
+   USE zdf_oce        ! ocean vertical physics
+   USE zdfdrg         ! ocean vertical physics: bottom friction
+   USE zdfddm         ! ocean vertical physics: double diffusion
+   USE eosbn2         ! equation of state
+   USE phycst         ! physical constants
    !
-   USE lib_mpp         ! distibuted memory computing library
-   USE in_out_manager  ! I/O manager
-   USE iom             ! I/O manager library
-   USE wrk_nemo        ! Memory allocation
+   USE lib_mpp        ! distibuted memory computing library
+   USE in_out_manager ! I/O manager
+   USE iom            ! I/O manager library
 
    IMPLICIT NONE
@@ -77,9 +76,7 @@
       INTEGER ::   ikbu, ikbv      ! local integers
       REAL(wp)::   zvm, zvt, zvs, z1_2rau0   ! local scalars
-      REAL(wp), POINTER, DIMENSION(:,:)  :: ztswu, ztswv, z2dx, z2dy   ! 2D workspace 
-      !!----------------------------------------------------------------------
-
-      CALL wrk_alloc( jpi, jpj, ztswu, ztswv, z2dx, z2dy )
-
+      REAL(wp), DIMENSION(jpi,jpj)  :: ztswu, ztswv, z2dx, z2dy   ! 2D workspace 
+      !!----------------------------------------------------------------------
+      !
       IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN
          !
@@ -123,8 +120,8 @@
                DO jj = 1, jpjm1
                   DO ji = 1, jpim1
-                     zvt = ptrdx(ji,jj,jk) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
-                        &                  * e1u    (ji  ,jj  ) * e2u    (ji,jj) * e3u_n(ji,jj,jk)
-                     zvs = ptrdy(ji,jj,jk) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
-                        &                  * e1v    (ji  ,jj  ) * e2v    (ji,jj) * e3u_n(ji,jj,jk)
+                     zvt = ptrdx(ji,jj,jk) * tmask_i(ji+1,jj) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
+                        &                                     * e1e2u  (ji,jj) * e3u_n(ji,jj,jk)
+                     zvs = ptrdy(ji,jj,jk) * tmask_i(ji,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
+                        &                                     * e1e2v  (ji,jj) * e3u_n(ji,jj,jk)
                      umo(ktrd) = umo(ktrd) + zvt
                      vmo(ktrd) = vmo(ktrd) + zvs
@@ -138,8 +135,8 @@
                DO jj = 1, jpjm1
                   DO ji = 1, jpim1
-                     zvt = ( utau_b(ji,jj) + utau(ji,jj) ) * tmask_i(ji+1,jj  ) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
-                        &                       * z1_2rau0 * e1u    (ji  ,jj  ) * e2u    (ji,jj)
-                     zvs = ( vtau_b(ji,jj) + vtau(ji,jj) ) * tmask_i(ji  ,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
-                        &                       * z1_2rau0 * e1v    (ji  ,jj  ) * e2v    (ji,jj) * e3u_n(ji,jj,jk)
+                     zvt = ( utau_b(ji,jj) + utau(ji,jj) ) * tmask_i(ji+1,jj) * tmask_i(ji,jj) * umask(ji,jj,jk)   &
+                        &                                                     * z1_2rau0       * e1e2u(ji,jj)
+                     zvs = ( vtau_b(ji,jj) + vtau(ji,jj) ) * tmask_i(ji,jj+1) * tmask_i(ji,jj) * vmask(ji,jj,jk)   &
+                        &                                                     * z1_2rau0       * e1e2v(ji,jj)
                      umo(jpdyn_tau) = umo(jpdyn_tau) + zvt
                      vmo(jpdyn_tau) = vmo(jpdyn_tau) + zvs
@@ -151,5 +148,5 @@
             IF( ktrd == jpdyn_atf ) THEN     ! last trend (asselin time filter)
                !
-               IF( ln_bfrimp ) THEN                   ! implicit bfr case: compute separately the bottom friction 
+               IF( ln_drgimp ) THEN                   ! implicit drag case: compute separately the bottom friction 
                   z1_2rau0 = 0.5_wp / rau0
                   DO jj = 1, jpjm1
@@ -157,6 +154,6 @@
                         ikbu = mbku(ji,jj)                  ! deepest ocean u- & v-levels
                         ikbv = mbkv(ji,jj)
-                        zvt = bfrua(ji,jj) * un(ji,jj,ikbu) * e1u(ji,jj) * e2v(ji,jj)
-                        zvs = bfrva(ji,jj) * vn(ji,jj,ikbv) * e1v(ji,jj) * e2v(ji,jj)
+                        zvt = 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * un(ji,jj,ikbu) * e1e2u(ji,jj)
+                        zvs = 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vn(ji,jj,ikbv) * e1e2v(ji,jj)
                         umo(jpdyn_bfri) = umo(jpdyn_bfri) + zvt
                         vmo(jpdyn_bfri) = vmo(jpdyn_bfri) + zvs
@@ -165,4 +162,5 @@
                   END DO
                ENDIF
+!!gm top drag case is missing 
                ! 
                CALL glo_dyn_wri( kt )                 ! print the results in ocean.output
@@ -178,6 +176,4 @@
       ENDIF
       !
-      CALL wrk_dealloc( jpi, jpj, ztswu, ztswv, z2dx, z2dy )
-      !
    END SUBROUTINE trd_glo
 
@@ -193,8 +189,6 @@
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       REAL(wp) ::   zcof         ! local scalar
-      REAL(wp), POINTER, DIMENSION(:,:,:)  ::  zkx, zky, zkz, zkepe  
-      !!----------------------------------------------------------------------
-
-      CALL wrk_alloc( jpi, jpj, jpk, zkx, zky, zkz, zkepe )
+      REAL(wp), DIMENSION(jpi,jpj,jpk)  ::  zkx, zky, zkz, zkepe  
+      !!----------------------------------------------------------------------
 
       ! I. Momentum trends
@@ -283,5 +277,5 @@
             &      + vmo(jpdyn_bfr) + vmo(jpdyn_atf) ) / tvolv
             WRITE (numout,9513) umo(jpdyn_tau) / tvolu, vmo(jpdyn_tau) / tvolv
-            IF( ln_bfrimp )   WRITE (numout,9514) umo(jpdyn_bfri) / tvolu, vmo(jpdyn_bfri) / tvolv
+            IF( ln_drgimp )   WRITE (numout,9514) umo(jpdyn_bfri) / tvolu, vmo(jpdyn_bfri) / tvolv
          ENDIF
 
@@ -322,5 +316,5 @@
             &      + hke(jpdyn_bfr) + hke(jpdyn_atf) ) / tvolt
             WRITE (numout,9533) hke(jpdyn_tau) / tvolt
-            IF( ln_bfrimp )   WRITE (numout,9534) hke(jpdyn_bfri) / tvolt
+            IF( ln_drgimp )   WRITE (numout,9534) hke(jpdyn_bfri) / tvolt
          ENDIF
 
@@ -372,6 +366,4 @@
       ENDIF
       !
-      CALL wrk_dealloc( jpi, jpj, jpk, zkx, zky, zkz, zkepe )
-      !
    END SUBROUTINE glo_dyn_wri
 
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/TRD/trdken.F90
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/TRD/trdken.F90	(revision 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/TRD/trdken.F90	(revision 8143)
@@ -13,11 +13,11 @@
    USE oce            ! ocean dynamics and tracers variables
    USE dom_oce        ! ocean space and time domain variables
+   USE phycst         ! physical constants
    USE sbc_oce        ! surface boundary condition: ocean
    USE zdf_oce        ! ocean vertical physics variables
+   USE zdfdrg         ! ocean vertical physics: bottom friction
+!!gm   USE dynhpg          ! hydrostatic pressure gradient   
+   USE ldftra         ! ocean active tracers lateral physics
    USE trd_oce        ! trends: ocean variables
-!!gm   USE dynhpg          ! hydrostatic pressure gradient   
-   USE zdfbfr         ! bottom friction
-   USE ldftra         ! ocean active tracers lateral physics
-   USE phycst         ! physical constants
    USE trdvor         ! ocean vorticity trends 
    USE trdglo         ! trends:global domain averaged
@@ -27,5 +27,4 @@
    USE iom            ! I/O manager library
    USE lib_mpp        ! MPP library
-   USE wrk_nemo       ! Memory allocation
    USE ldfslp         ! Isopycnal slopes
 
@@ -74,5 +73,5 @@
       !!          diagnose separately the KE trend associated with wind stress
       !!              - bottom friction case (jpdyn_bfr):
-      !!          explicit case (ln_bfrimp=F): bottom trend put in the 1st level 
+      !!          explicit case (ln_drgimp=F): bottom trend put in the 1st level 
       !!                                       of putrd, pvtrd
       !
@@ -86,9 +85,7 @@
       INTEGER ::   ikbu  , ikbv     ! local integers
       INTEGER ::   ikbum1, ikbvm1   !   -       -
-      REAL(wp), POINTER, DIMENSION(:,:)   ::   z2dx, z2dy, zke2d   ! 2D workspace 
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zke                 ! 3D workspace 
-      !!----------------------------------------------------------------------
-      !
-      CALL wrk_alloc( jpi, jpj, jpk, zke )
+      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   z2dx, z2dy, zke2d   ! 2D workspace 
+      REAL(wp), DIMENSION(jpi,jpj,jpk)      ::   zke                 ! 3D workspace 
+      !!----------------------------------------------------------------------
       !
       CALL lbc_lnk( putrd, 'U', -1. )   ;   CALL lbc_lnk( pvtrd, 'V', -1. )      ! lateral boundary conditions
@@ -125,5 +122,5 @@
          CASE( jpdyn_zdf )   ;   CALL iom_put( "ketrd_zdf"   , zke )    ! vertical diffusion 
          !                   !                                          ! wind stress trends
-                                 CALL wrk_alloc( jpi, jpj, z2dx, z2dy, zke2d )
+                                 ALLOCATE( z2dx(jpi,jpj) , z2dy(jpi,jpj) , zke2d(jpi,jpj) )
                            z2dx(:,:) = un(:,:,1) * ( utau_b(:,:) + utau(:,:) ) * e1e2u(:,:) * umask(:,:,1)
                            z2dy(:,:) = vn(:,:,1) * ( vtau_b(:,:) + vtau(:,:) ) * e1e2v(:,:) * vmask(:,:,1)
@@ -136,9 +133,9 @@
                            END DO
                                  CALL iom_put( "ketrd_tau"   , zke2d )  ! 
-                                 CALL wrk_dealloc( jpi, jpj     , z2dx, z2dy, zke2d )
+                                 DEALLOCATE( z2dx , z2dy , zke2d )
          CASE( jpdyn_bfr )   ;   CALL iom_put( "ketrd_bfr"   , zke )    ! bottom friction (explicit case) 
 !!gm TO BE DONE properly
-!!gm only valid if ln_bfrimp=F otherwise the bottom stress as to be recomputed at the end of the computation....
-!         IF(.NOT. ln_bfrimp) THEN
+!!gm only valid if ln_drgimp=F otherwise the bottom stress as to be recomputed at the end of the computation....
+!         IF(.NOT. ln_drgimp) THEN
 !            DO jj = 1, jpj    !   
 !               DO ji = 1, jpi
@@ -163,5 +160,5 @@
 !! reflechir a une possible sauvegarde du "vrai" un,vn pour le calcul de atf....
 !
-!         IF( ln_bfrimp ) THEN                                          ! bottom friction (implicit case)
+!         IF( ln_drgimp ) THEN                                          ! bottom friction (implicit case)
 !            DO jj = 1, jpj                                                  ! after velocity known (now filed at this stage)
 !               DO ji = 1, jpi
@@ -192,6 +189,4 @@
       END SELECT
       !
-      CALL wrk_dealloc( jpi, jpj, jpk, zke )
-      !
    END SUBROUTINE trd_ken
 
@@ -207,15 +202,12 @@
       !! ** Work only for full steps and partial steps (ln_hpg_zco or ln_hpg_zps)
       !!---------------------------------------------------------------------- 
-      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
-      !!
-      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   pconv
-      !
-      INTEGER  ::   ji, jj, jk                       ! dummy loop indices
-      INTEGER  ::   iku, ikv                         ! temporary integers
-      REAL(wp) ::   zcoef                            ! temporary scalars
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zconv  ! temporary conv on W-grid
-      !!----------------------------------------------------------------------
-      !
-      CALL wrk_alloc( jpi,jpj,jpk, zconv )
+      INTEGER                   , INTENT(in   ) ::   kt      ! ocean time-step index
+      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pconv   ! 
+      !
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      INTEGER  ::   iku, ikv     ! local integers
+      REAL(wp) ::   zcoef        ! local scalars
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zconv  ! 3D workspace
+      !!----------------------------------------------------------------------
       !
       ! Local constant initialization 
@@ -240,6 +232,4 @@
       END DO
       !
-      CALL wrk_dealloc( jpi,jpj,jpk, zconv )      
-      !
    END SUBROUTINE ken_p2k
 
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 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90	(revision 8143)
@@ -4,6 +4,7 @@
    !! Ocean physics : define vertical mixing variables
    !!=====================================================================
-   !! history :  1.0  !  2002-06  (G. Madec) Original code
-   !!            3.2  !  2009-07  (G.Madec) addition of avm
+   !! history :  1.0  !  2002-06  (G. Madec)  Original code
+   !!            3.2  !  2009-07  (G. Madec)  addition of avm
+   !!            4.0  !  2017-05  (G. Madec)  drag coef. defined at t-point (zdfdrg.F90)
    !!----------------------------------------------------------------------
    USE par_oce        ! ocean parameters
@@ -54,8 +55,4 @@
    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   avtb_2d        !: horizontal shape of background Kz profile
    REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::   avmb , avtb    !: background profile of avm and avt
-!!gm
-   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   bfrua, bfrva   !: bottom friction coefficients
-   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   tfrua, tfrva   !: top    friction coefficients
-!!gm
 
    !!----------------------------------------------------------------------
@@ -71,9 +68,8 @@
       !!----------------------------------------------------------------------
       !
-      ALLOCATE( avm  (jpi,jpj,jpk) , avt  (jpi,jpj,jpk) , avs(jpi,jpj,jpk) ,        &
-         &      avm_k(jpi,jpj,jpk) , avt_k(jpi,jpj,jpk) , en (jpi,jpj,jpk) ,        & 
-         &      avmb(jpk) , bfrua(jpi,jpj) , tfrua(jpi, jpj) ,                      &
-         &      avtb(jpk) , bfrva(jpi,jpj) , tfrva(jpi, jpj) , avtb_2d(jpi,jpj) ,   &
-         &      avmu(jpi,jpj,jpk), avmv(jpi,jpj,jpk) ,  STAT = zdf_oce_alloc )
+      ALLOCATE( avm  (jpi,jpj,jpk) , avt  (jpi,jpj,jpk) , avs(jpi,jpj,jpk) ,   &
+         &      avm_k(jpi,jpj,jpk) , avt_k(jpi,jpj,jpk) , en (jpi,jpj,jpk) ,   & 
+         &      avmb(jpk) , avtb(jpk) ,  avtb_2d(jpi,jpj) ,                    &
+         &      avmu(jpi,jpj,jpk), avmv(jpi,jpj,jpk)    , STAT = zdf_oce_alloc )
          !
       IF( zdf_oce_alloc /= 0 )   CALL ctl_warn('zdf_oce_alloc: failed to allocate arrays')
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90	(revision 8093)
+++ 	(revision )
@@ -1,486 +1,0 @@
-MODULE zdfbfr
-   !!======================================================================
-   !!                       ***  MODULE  zdfbfr  ***
-   !! Ocean physics: Bottom friction
-   !!======================================================================
-   !! History :  OPA  ! 1997-06  (G. Madec, A.-M. Treguier)  Original code
-   !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module
-   !!            3.2  ! 2009-09  (A.C.Coward)  Correction to include barotropic contribution
-   !!            3.3  ! 2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
-   !!            3.4  ! 2011-11  (H. Liu) implementation of semi-implicit bottom friction option
-   !!                 ! 2012-06  (H. Liu) implementation of Log Layer bottom friction option
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   zdf_bfr      : update bottom friction coefficient (non-linear bottom friction only)
-   !!   zdf_bfr_init : read in namelist and control the bottom friction parameters.
-   !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers variables
-   USE dom_oce         ! ocean space and time domain variables
-   USE zdf_oce         ! ocean vertical physics variables
-   USE in_out_manager  ! I/O manager
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE lib_mpp         ! distributed memory computing
-   USE prtctl          ! Print control
-   USE timing          ! Timing
-   USE wrk_nemo        ! Memory Allocation
-   USE phycst, ONLY: vkarmn
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   zdf_bfr         ! called by step.F90
-   PUBLIC   zdf_bfr_init    ! called by nemogcm.F90
-
-   !                                 !!* Namelist nambfr: bottom friction namelist *
-   INTEGER , PUBLIC ::   nn_bfr       ! = 0/1/2/3 type of bottom friction  (PUBLIC for TAM)
-   REAL(wp), PUBLIC ::   rn_bfri1     ! bottom drag coefficient (linear case)  (PUBLIC for TAM)
-   REAL(wp), PUBLIC ::   rn_bfri2     ! bottom drag coefficient (non linear case) (PUBLIC for TAM)
-   REAL(wp), PUBLIC ::   rn_bfri2_max ! Maximum bottom drag coefficient (non linear case and ln_loglayer=T) (PUBLIC for TAM)
-   REAL(wp), PUBLIC ::   rn_bfeb2     ! background bottom turbulent kinetic energy  [m2/s2] (PUBLIC for TAM)
-   REAL(wp), PUBLIC ::   rn_bfrien    ! local factor to enhance coefficient bfri (PUBLIC for TAM)
-   LOGICAL , PUBLIC ::   ln_bfr2d     ! logical switch for 2D enhancement (PUBLIC for TAM)
-   REAL(wp), PUBLIC ::   rn_tfri1     ! top drag coefficient (linear case)  (PUBLIC for TAM)
-   REAL(wp), PUBLIC ::   rn_tfri2     ! top drag coefficient (non linear case) (PUBLIC for TAM)
-   REAL(wp), PUBLIC ::   rn_tfri2_max ! Maximum top drag coefficient (non linear case and ln_loglayer=T) (PUBLIC for TAM)
-   REAL(wp), PUBLIC ::   rn_tfeb2     ! background top turbulent kinetic energy  [m2/s2] (PUBLIC for TAM)
-   REAL(wp), PUBLIC ::   rn_tfrien    ! local factor to enhance coefficient tfri (PUBLIC for TAM)
-   LOGICAL , PUBLIC ::   ln_tfr2d     ! logical switch for 2D enhancement (PUBLIC for TAM)
-
-   LOGICAL , PUBLIC ::   ln_loglayer  ! switch for log layer bfr coeff. (PUBLIC for TAM)
-   REAL(wp), PUBLIC ::   rn_bfrz0     ! bottom roughness for loglayer bfr coeff (PUBLIC for TAM)
-   REAL(wp), PUBLIC ::   rn_tfrz0     ! bottom roughness for loglayer bfr coeff (PUBLIC for TAM)
-   LOGICAL , PUBLIC ::   ln_bfrimp    ! logical switch for implicit bottom friction
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::  bfrcoef2d, tfrcoef2d   ! 2D bottom/top drag coefficient (PUBLIC for TAM)
-
-   !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
-   !! $Id$
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   INTEGER FUNCTION zdf_bfr_alloc()
-      !!----------------------------------------------------------------------
-      !!                ***  FUNCTION zdf_bfr_alloc  ***
-      !!----------------------------------------------------------------------
-      ALLOCATE( bfrcoef2d(jpi,jpj), tfrcoef2d(jpi,jpj), STAT=zdf_bfr_alloc )
-      !
-      IF( lk_mpp             )   CALL mpp_sum ( zdf_bfr_alloc )
-      IF( zdf_bfr_alloc /= 0 )   CALL ctl_warn('zdf_bfr_alloc: failed to allocate arrays.')
-   END FUNCTION zdf_bfr_alloc
-
-
-   SUBROUTINE zdf_bfr( kt )
-      !!----------------------------------------------------------------------
-      !!                   ***  ROUTINE zdf_bfr  ***
-      !!
-      !! ** Purpose :   compute the bottom friction coefficient.
-      !!
-      !! ** Method  :   Calculate and store part of the momentum trend due
-      !!              to bottom friction following the chosen friction type
-      !!              (free-slip, linear, or quadratic). The component
-      !!              calculated here is multiplied by the bottom velocity in
-      !!              dyn_bfr to provide the trend term.
-      !!                The coefficients are updated at each time step only
-      !!              in the quadratic case.
-      !!
-      !! ** Action  :   bfrua , bfrva   bottom friction coefficients
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
-      !!
-      INTEGER  ::   ji, jj                       ! dummy loop indices
-      INTEGER  ::   ikbt, ikbu, ikbv             ! local integers
-      REAL(wp) ::   zvu, zuv, zecu, zecv, ztmp   ! temporary scalars
-      REAL(wp), POINTER, DIMENSION(:,:) ::  zbfrt, ztfrt
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('zdf_bfr')
-      !
-      IF( nn_bfr == 2 ) THEN                 ! quadratic bottom friction only
-         !
-         CALL wrk_alloc( jpi, jpj, zbfrt, ztfrt )
-
-         IF ( ln_loglayer.AND. .NOT.ln_linssh ) THEN ! "log layer" bottom friction coefficient
-
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ikbt = mbkt(ji,jj)
-!! JC: possible WAD implementation should modify line below if layers vanish
-                  ztmp = tmask(ji,jj,ikbt) * ( vkarmn / LOG( 0.5_wp * e3t_n(ji,jj,ikbt) / rn_bfrz0 ))**2._wp
-                  zbfrt(ji,jj) = MAX(bfrcoef2d(ji,jj), ztmp)
-                  zbfrt(ji,jj) = MIN(zbfrt(ji,jj), rn_bfri2_max)
-               END DO
-            END DO
-! (ISF)
-            IF ( ln_isfcav ) THEN
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     ikbt = mikt(ji,jj)
-! JC: possible WAD implementation should modify line below if layers vanish
-                     ztmp = (1.-tmask(ji,jj,1)) * ( vkarmn / LOG( 0.5_wp * e3t_n(ji,jj,ikbt) / rn_bfrz0 ))**2._wp
-                     ztfrt(ji,jj) = MAX(tfrcoef2d(ji,jj), ztmp)
-                     ztfrt(ji,jj) = MIN(ztfrt(ji,jj), rn_tfri2_max)
-                  END DO
-               END DO
-            END IF
-            !   
-         ELSE
-            zbfrt(:,:) = bfrcoef2d(:,:)
-            ztfrt(:,:) = tfrcoef2d(:,:)
-         ENDIF
-
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
-               ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points
-               ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
-               !
-               zvu  = 0.25 * (  vn(ji,jj  ,ikbu) + vn(ji+1,jj  ,ikbu)     &
-                  &           + vn(ji,jj-1,ikbu) + vn(ji+1,jj-1,ikbu)  )
-               zuv  = 0.25 * (  un(ji,jj  ,ikbv) + un(ji-1,jj  ,ikbv)     &
-                  &           + un(ji,jj+1,ikbv) + un(ji-1,jj+1,ikbv)  )
-               !
-               zecu = SQRT(  un(ji,jj,ikbu) * un(ji,jj,ikbu) + zvu*zvu + rn_bfeb2  )
-               zecv = SQRT(  vn(ji,jj,ikbv) * vn(ji,jj,ikbv) + zuv*zuv + rn_bfeb2  )
-               !
-               bfrua(ji,jj) = - 0.5_wp * ( zbfrt(ji,jj) + zbfrt(ji+1,jj  ) ) * zecu
-               bfrva(ji,jj) = - 0.5_wp * ( zbfrt(ji,jj) + zbfrt(ji  ,jj+1) ) * zecv
-               !
-               ! in case of 2 cell water column, we assume each cell feels the top and bottom friction
-               IF ( ln_isfcav ) THEN
-                  IF ( miku(ji,jj) + 1  >=  mbku(ji,jj) ) THEN
-                     bfrua(ji,jj) = - 0.5_wp * ( ( zbfrt(ji,jj) + zbfrt(ji+1,jj  ) )   &
-                                  &            + ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) ) ) &
-                                  &          * zecu * (1._wp - umask(ji,jj,1))
-                  ENDIF
-                  IF( mikv(ji,jj) + 1  >=  mbkv(ji,jj) ) THEN
-                     bfrva(ji,jj) = - 0.5_wp * ( ( zbfrt(ji,jj) + zbfrt(ji  ,jj+1) )   &
-                                  &            + ( ztfrt(ji,jj) + ztfrt(ji  ,jj+1) ) ) &
-                                  &          * zecv * (1._wp - vmask(ji,jj,1))
-                  ENDIF
-               ENDIF
-            END DO
-         END DO
-         CALL lbc_lnk( bfrua, 'U', 1. )   ;   CALL lbc_lnk( bfrva, 'V', 1. )      ! Lateral boundary condition
-
-         IF( ln_isfcav ) THEN
-            DO jj = 2, jpjm1
-               DO ji = 2, jpim1
-                  ! (ISF) ========================================================================
-                  ikbu = miku(ji,jj)         ! ocean top level at u- and v-points 
-                  ikbv = mikv(ji,jj)         ! (1st wet ocean u- and v-points)
-                  !
-                  zvu  = 0.25 * (  vn(ji,jj  ,ikbu) + vn(ji+1,jj  ,ikbu)     &
-                     &           + vn(ji,jj-1,ikbu) + vn(ji+1,jj-1,ikbu)  )
-                  zuv  = 0.25 * (  un(ji,jj  ,ikbv) + un(ji-1,jj  ,ikbv)     &
-                     &           + un(ji,jj+1,ikbv) + un(ji-1,jj+1,ikbv)  )
-              !
-                  zecu = SQRT(  un(ji,jj,ikbu) * un(ji,jj,ikbu) + zvu*zvu + rn_tfeb2 )
-                  zecv = SQRT(  vn(ji,jj,ikbv) * vn(ji,jj,ikbv) + zuv*zuv + rn_tfeb2 )
-              !
-                  tfrua(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) ) * zecu * (1._wp - umask(ji,jj,1))
-                  tfrva(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji  ,jj+1) ) * zecv * (1._wp - vmask(ji,jj,1))
-              ! (ISF) END ====================================================================
-              ! in case of 2 cell water column, we assume each cell feels the top and bottom friction
-                  IF ( miku(ji,jj) + 1 .GE. mbku(ji,jj) ) THEN
-                     tfrua(ji,jj) = - 0.5_wp * ( ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) )   &
-                                  &            + ( zbfrt(ji,jj) + zbfrt(ji+1,jj  ) ) ) &
-                                  &          * zecu * (1._wp - umask(ji,jj,1))
-                  END IF
-                  IF ( mikv(ji,jj) + 1 .GE. mbkv(ji,jj) ) THEN
-                     tfrva(ji,jj) = - 0.5_wp * ( ( ztfrt(ji,jj) + ztfrt(ji  ,jj+1) )   &
-                                  &            + ( zbfrt(ji,jj) + zbfrt(ji  ,jj+1) ) ) &
-                                  &          * zecv * (1._wp - vmask(ji,jj,1))
-                  END IF
-               END DO
-            END DO
-            CALL lbc_lnk( tfrua, 'U', 1. )   ;   CALL lbc_lnk( tfrva, 'V', 1. )      ! Lateral boundary condition
-         END IF
-         !
-         !
-         IF(ln_ctl)   CALL prt_ctl( tab2d_1=bfrua, clinfo1=' bfr  - u: ', mask1=umask,        &
-            &                       tab2d_2=bfrva, clinfo2=       ' v: ', mask2=vmask,ovlap=1 )
-         CALL wrk_dealloc( jpi,jpj,zbfrt, ztfrt )
-      ENDIF
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('zdf_bfr')
-      !
-   END SUBROUTINE zdf_bfr
-
-
-   SUBROUTINE zdf_bfr_init
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE zdf_bfr_init  ***
-      !!
-      !! ** Purpose :   Initialization of the bottom friction
-      !!
-      !! ** Method  :   Read the nambfr namelist and check their consistency
-      !!                called at the first timestep (nit000)
-      !!----------------------------------------------------------------------
-      USE iom   ! I/O module for ehanced bottom friction file
-      !!
-      INTEGER   ::   inum         ! logical unit for enhanced bottom friction file
-      INTEGER   ::   ji, jj       ! dummy loop indexes
-      INTEGER   ::   ikbt, ikbu, ikbv   ! temporary integers
-      INTEGER   ::   ictu, ictv         !    -          -
-      INTEGER   ::   ios
-      REAL(wp)  ::   zminbfr, zmaxbfr   ! temporary scalars
-      REAL(wp)  ::   zmintfr, zmaxtfr   ! temporary scalars
-      REAL(wp)  ::   ztmp, zfru, zfrv   !    -         -
-      !!
-      NAMELIST/nambfr/ nn_bfr, rn_bfri1, rn_bfri2, rn_bfri2_max, rn_bfeb2, rn_bfrz0, ln_bfr2d, &
-                    &          rn_tfri1, rn_tfri2, rn_tfri2_max, rn_tfeb2, rn_tfrz0, ln_tfr2d, &
-                    &  rn_bfrien, rn_tfrien, ln_bfrimp, ln_loglayer
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('zdf_bfr_init')
-      !
-      !                              !* Allocate zdfbfr arrays
-      IF( zdf_bfr_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_bfr_init : unable to allocate arrays' )
-      !
-      !                              !* Parameter control and print
-      !
-      REWIND( numnam_ref )              ! Namelist nambfr in reference namelist : Bottom momentum boundary condition
-      READ  ( numnam_ref, nambfr, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambfr in reference namelist', lwp )
-
-      REWIND( numnam_cfg )              ! Namelist nambfr in configuration namelist : Bottom momentum boundary condition
-      READ  ( numnam_cfg, nambfr, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambfr in configuration namelist', lwp )
-      IF(lwm) WRITE ( numond, nambfr )
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'zdf_bfr_init : momentum bottom friction'
-      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
-      IF(lwp) WRITE(numout,*) '   Namelist nam_bfr : set bottom friction parameters'
-      !
-      SELECT CASE (nn_bfr)
-      !
-      CASE( 0 )
-         IF(lwp) WRITE(numout,*) '      free-slip '
-         bfrua(:,:) = 0._wp
-         bfrva(:,:) = 0._wp
-         tfrua(:,:) = 0._wp
-         tfrva(:,:) = 0._wp
-         !
-      CASE( 1 )
-         IF(lwp) WRITE(numout,*) '      linear botton friction'
-         IF(lwp) WRITE(numout,*) '      bottom friction coef.   rn_bfri1  = ', rn_bfri1
-         IF( ln_bfr2d ) THEN
-            IF(lwp) WRITE(numout,*) '      read coef. enhancement distribution from file   ln_bfr2d  = ', ln_bfr2d
-            IF(lwp) WRITE(numout,*) '      coef rn_bfri2 enhancement factor                rn_bfrien  = ',rn_bfrien
-         ENDIF
-         IF ( ln_isfcav ) THEN
-            IF(lwp) WRITE(numout,*) '      top    friction coef.   rn_bfri1  = ', rn_tfri1
-            IF( ln_tfr2d ) THEN
-               IF(lwp) WRITE(numout,*) '      read coef. enhancement distribution from file   ln_tfr2d  = ', ln_tfr2d
-               IF(lwp) WRITE(numout,*) '      coef rn_tfri2 enhancement factor                rn_tfrien  = ',rn_tfrien
-            ENDIF
-         END IF
-         !
-         IF(ln_bfr2d) THEN
-            ! bfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement
-            CALL iom_open('bfr_coef.nc',inum)
-            CALL iom_get (inum, jpdom_data, 'bfr_coef',bfrcoef2d,1) ! bfrcoef2d is used as tmp array
-            CALL iom_close(inum)
-            bfrcoef2d(:,:) = rn_bfri1 * ( 1 + rn_bfrien * bfrcoef2d(:,:) )
-         ELSE
-            bfrcoef2d(:,:) = rn_bfri1  ! initialize bfrcoef2d to the namelist variable
-         ENDIF
-         !
-         bfrua(:,:) = - bfrcoef2d(:,:)
-         bfrva(:,:) = - bfrcoef2d(:,:)
-         !
-         IF ( ln_isfcav ) THEN
-            IF(ln_tfr2d) THEN
-               ! tfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement
-               CALL iom_open('tfr_coef.nc',inum)
-               CALL iom_get (inum, jpdom_data, 'tfr_coef',tfrcoef2d,1) ! tfrcoef2d is used as tmp array
-               CALL iom_close(inum)
-               tfrcoef2d(:,:) = rn_tfri1 * ( 1 + rn_tfrien * tfrcoef2d(:,:) )
-            ELSE
-               tfrcoef2d(:,:) = rn_tfri1  ! initialize tfrcoef2d to the namelist variable
-            ENDIF
-            !
-            tfrua(:,:) = - tfrcoef2d(:,:)
-            tfrva(:,:) = - tfrcoef2d(:,:)
-         END IF
-         !
-      CASE( 2 )
-         IF(lwp) WRITE(numout,*) '      quadratic bottom friction'
-         IF(lwp) WRITE(numout,*) '      friction coef.   rn_bfri2  = ', rn_bfri2
-         IF(lwp) WRITE(numout,*) '      Max. coef. (log case)   rn_bfri2_max  = ', rn_bfri2_max
-         IF(lwp) WRITE(numout,*) '      background tke   rn_bfeb2  = ', rn_bfeb2
-         IF(lwp) WRITE(numout,*) '      log formulation   ln_bfr2d = ', ln_loglayer
-         IF(lwp) WRITE(numout,*) '      bottom roughness  rn_bfrz0 [m] = ', rn_bfrz0
-         IF( rn_bfrz0 <= 0._wp ) THEN
-            WRITE(ctmp1,*) '      bottom roughness must be strictly positive'
-            CALL ctl_stop( ctmp1 )
-         ENDIF
-         IF( ln_bfr2d ) THEN
-            IF(lwp) WRITE(numout,*) '      read coef. enhancement distribution from file   ln_bfr2d  = ', ln_bfr2d
-            IF(lwp) WRITE(numout,*) '      coef rn_bfri2 enhancement factor                rn_bfrien  = ',rn_bfrien
-         ENDIF
-         IF ( ln_isfcav ) THEN
-            IF(lwp) WRITE(numout,*) '      quadratic top    friction'
-            IF(lwp) WRITE(numout,*) '      friction coef.    rn_tfri2     = ', rn_tfri2
-            IF(lwp) WRITE(numout,*) '      Max. coef. (log case)   rn_tfri2_max  = ', rn_tfri2_max
-            IF(lwp) WRITE(numout,*) '      background tke    rn_tfeb2     = ', rn_tfeb2
-            IF(lwp) WRITE(numout,*) '      log formulation   ln_tfr2d     = ', ln_loglayer
-            IF(lwp) WRITE(numout,*) '      top roughness     rn_tfrz0 [m] = ', rn_tfrz0
-            IF( rn_tfrz0 <= 0._wp ) THEN
-               WRITE(ctmp1,*) '      top roughness must be strictly positive'
-               CALL ctl_stop( ctmp1 )
-            ENDIF
-            IF( ln_tfr2d ) THEN
-               IF(lwp) WRITE(numout,*) '      read coef. enhancement distribution from file   ln_tfr2d  = ', ln_tfr2d
-               IF(lwp) WRITE(numout,*) '      coef rn_tfri2 enhancement factor                rn_tfrien  = ',rn_tfrien
-            ENDIF
-         END IF
-         !
-         IF(ln_bfr2d) THEN
-            ! bfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement
-            CALL iom_open('bfr_coef.nc',inum)
-            CALL iom_get (inum, jpdom_data, 'bfr_coef',bfrcoef2d,1) ! bfrcoef2d is used as tmp array
-            CALL iom_close(inum)
-            !
-            bfrcoef2d(:,:) = rn_bfri2 * ( 1 + rn_bfrien * bfrcoef2d(:,:) )
-         ELSE
-            bfrcoef2d(:,:) = rn_bfri2  ! initialize bfrcoef2d to the namelist variable
-         ENDIF
-         
-         IF ( ln_isfcav ) THEN
-            IF(ln_tfr2d) THEN
-               ! tfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement
-               CALL iom_open('tfr_coef.nc',inum)
-               CALL iom_get (inum, jpdom_data, 'tfr_coef',tfrcoef2d,1) ! tfrcoef2d is used as tmp array
-               CALL iom_close(inum)
-               !
-               tfrcoef2d(:,:) = rn_tfri2 * ( 1 + rn_tfrien * tfrcoef2d(:,:) )
-            ELSE
-               tfrcoef2d(:,:) = rn_tfri2  ! initialize tfrcoef2d to the namelist variable
-            ENDIF
-         END IF
-         !
-         IF( ln_loglayer.AND. ln_linssh ) THEN ! set "log layer" bottom friction once for all
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ikbt = mbkt(ji,jj)
-                  ztmp = tmask(ji,jj,ikbt) * ( vkarmn / LOG( 0.5_wp * e3t_n(ji,jj,ikbt) / rn_bfrz0 ))**2._wp
-                  bfrcoef2d(ji,jj) = MAX(bfrcoef2d(ji,jj), ztmp)
-                  bfrcoef2d(ji,jj) = MIN(bfrcoef2d(ji,jj), rn_bfri2_max)
-               END DO
-            END DO
-            IF ( ln_isfcav ) THEN
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     ikbt = mikt(ji,jj)
-                     ztmp = tmask(ji,jj,ikbt) * ( vkarmn / LOG( 0.5_wp * e3t_n(ji,jj,ikbt) / rn_tfrz0 ))**2._wp
-                     tfrcoef2d(ji,jj) = MAX(tfrcoef2d(ji,jj), ztmp)
-                     tfrcoef2d(ji,jj) = MIN(tfrcoef2d(ji,jj), rn_tfri2_max)
-                  END DO
-               END DO
-            END IF
-         ENDIF
-         !
-      CASE DEFAULT
-         IF(lwp) WRITE(ctmp1,*) '         bad flag value for nn_bfr = ', nn_bfr
-         CALL ctl_stop( ctmp1 )
-         !
-      END SELECT
-      !
-      IF(lwp) WRITE(numout,*) '      implicit bottom friction switch                ln_bfrimp  = ', ln_bfrimp
-      !
-      !                              ! Make sure ln_zdfexp=.false. when use implicit bfr
-      IF( ln_bfrimp .AND. ln_zdfexp ) THEN
-         IF(lwp) THEN
-            WRITE(numout,*)
-            WRITE(numout,*) 'Implicit bottom friction can only be used when ln_zdfexp=.false.'
-            WRITE(numout,*) '         but you set: ln_bfrimp=.true. and ln_zdfexp=.true.!!!!'
-            WRITE(ctmp1,*)  '         set either ln_zdfexp = .false or ln_bfrimp = .false.'
-            CALL ctl_stop( ctmp1 )
-         END IF
-      END IF
-      !
-      ! Basic stability check on bottom friction coefficient
-      !
-      ictu = 0               ! counter for stability criterion breaches at U-pts
-      ictv = 0               ! counter for stability criterion breaches at V-pts
-      zminbfr =  1.e10_wp    ! initialise tracker for minimum of bottom friction coefficient
-      zmaxbfr = -1.e10_wp    ! initialise tracker for maximum of bottom friction coefficient
-      zmintfr =  1.e10_wp    ! initialise tracker for minimum of bottom friction coefficient
-      zmaxtfr = -1.e10_wp    ! initialise tracker for maximum of bottom friction coefficient
-      !
-      DO jj = 2, jpjm1
-         DO ji = 2, jpim1
-             ikbu = mbku(ji,jj)       ! deepest ocean level at u- and v-points
-             ikbv = mbkv(ji,jj)
-             zfru = 0.5 * e3u_n(ji,jj,ikbu) / rdt
-             zfrv = 0.5 * e3v_n(ji,jj,ikbv) / rdt
-             IF( ABS( bfrcoef2d(ji,jj) ) > zfru ) THEN
-                IF( ln_ctl ) THEN
-                   WRITE(numout,*) 'BFR ', narea, nimpp+ji, njmpp+jj, ikbu
-                   WRITE(numout,*) 'BFR ', ABS( bfrcoef2d(ji,jj) ), zfru
-                ENDIF
-                ictu = ictu + 1
-             ENDIF
-             IF( ABS( bfrcoef2d(ji,jj) ) > zfrv ) THEN
-                 IF( ln_ctl ) THEN
-                     WRITE(numout,*) 'BFR ', narea, nimpp+ji, njmpp+jj, ikbv
-                     WRITE(numout,*) 'BFR ', bfrcoef2d(ji,jj), zfrv
-                 ENDIF
-                 ictv = ictv + 1
-             ENDIF
-             zminbfr = MIN(  zminbfr, MIN( zfru, ABS( bfrcoef2d(ji,jj) ) )  )
-             zmaxbfr = MAX(  zmaxbfr, MIN( zfrv, ABS( bfrcoef2d(ji,jj) ) )  )
-! (ISF)
-             IF ( ln_isfcav ) THEN
-                ikbu = miku(ji,jj)       ! 1st wet ocean level at u- and v-points
-                ikbv = mikv(ji,jj)
-                zfru = 0.5 * e3u_n(ji,jj,ikbu) / rdt
-                zfrv = 0.5 * e3v_n(ji,jj,ikbv) / rdt
-                IF( ABS( tfrcoef2d(ji,jj) ) > zfru ) THEN
-                   IF( ln_ctl ) THEN
-                      WRITE(numout,*) 'TFR ', narea, nimpp+ji, njmpp+jj, ikbu
-                      WRITE(numout,*) 'TFR ', ABS( tfrcoef2d(ji,jj) ), zfru
-                   ENDIF
-                   ictu = ictu + 1
-                ENDIF
-                IF( ABS( tfrcoef2d(ji,jj) ) > zfrv ) THEN
-                   IF( ln_ctl ) THEN
-                      WRITE(numout,*) 'TFR ', narea, nimpp+ji, njmpp+jj, ikbv
-                      WRITE(numout,*) 'TFR ', tfrcoef2d(ji,jj), zfrv
-                   ENDIF
-                   ictv = ictv + 1
-                ENDIF
-                zmintfr = MIN(  zmintfr, MIN( zfru, ABS( tfrcoef2d(ji,jj) ) )  )
-                zmaxtfr = MAX(  zmaxtfr, MIN( zfrv, ABS( tfrcoef2d(ji,jj) ) )  )
-             END IF
-! END ISF
-         END DO
-      END DO
-      IF( lk_mpp ) THEN
-         CALL mpp_sum( ictu )
-         CALL mpp_sum( ictv )
-         CALL mpp_min( zminbfr )
-         CALL mpp_max( zmaxbfr )
-         IF ( ln_isfcav) CALL mpp_min( zmintfr )
-         IF ( ln_isfcav) CALL mpp_max( zmaxtfr )
-      ENDIF
-      IF( .NOT.ln_bfrimp) THEN
-      IF( lwp .AND. ictu + ictv > 0 ) THEN
-         WRITE(numout,*) ' Bottom/Top friction stability check failed at ', ictu, ' U-points '
-         WRITE(numout,*) ' Bottom/Top friction stability check failed at ', ictv, ' V-points '
-         WRITE(numout,*) ' Bottom friction coefficient now ranges from: ', zminbfr, ' to ', zmaxbfr
-         IF ( ln_isfcav ) WRITE(numout,*) ' Top friction coefficient now ranges from: ', zmintfr, ' to ', zmaxtfr
-         WRITE(numout,*) ' Bottom/Top friction coefficient will be reduced where necessary'
-      ENDIF
-      ENDIF
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('zdf_bfr_init')
-      !
-   END SUBROUTINE zdf_bfr_init
-
-   !!======================================================================
-END MODULE zdfbfr
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfdrg.F90
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfdrg.F90	(revision 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfdrg.F90	(revision 8143)
@@ -10,4 +10,6 @@
    !!            3.4  ! 2011-11  (H. Liu) implementation of semi-implicit bottom friction option
    !!                 ! 2012-06  (H. Liu) implementation of Log Layer bottom friction option
+   !!            4.0  ! 2017-05  (G. Madec) zdfbfr becomes zdfdrg + variable names change
+   !!                                     + drag defined at t-point + new user interface + top drag (ocean cavities)
    !!----------------------------------------------------------------------
 
@@ -28,5 +30,4 @@
    USE prtctl         ! Print control
    USE timing         ! Timing
-   USE wrk_nemo       ! Memory Allocation
 
    IMPLICIT NONE
@@ -137,8 +138,4 @@
       ENDIF
       !
-!!gm to be moved at the end of zdfphy
-      CALL lbc_lnk( pCdU, 'T', 1. )    ! Lateral boundary condition
-!!gm end
-      !
       IF(ln_ctl)   CALL prt_ctl( tab2d_1=pCdU, clinfo1=' Cd*U ')
       !
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 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 8143)
@@ -8,4 +8,5 @@
    !!            3.3  !  2010-10  (C. Bricaud)  Add in the reference
    !!            4.0  !  2017-04  (G. Madec)  remove CPP keys & avm at t-point only 
+   !!             -   !  2017-05  (G. Madec)  add top friction as boundary condition
    !!----------------------------------------------------------------------
 
@@ -19,10 +20,6 @@
    USE domvvl         ! ocean space and time domain : variable volume layer
    USE zdf_oce        ! ocean vertical physics
-!!gm old
-   USE zdfbfr  , ONLY : rn_tfrz0, rn_bfrz0   ! top/bottom roughness
-!!gm new
    USE zdfdrg  , ONLY : r_z0_top , r_z0_bot   ! top/bottom roughness
    USE zdfdrg  , ONLY : rCdU_top , rCdU_bot   ! top/bottom friction
-!!gm
    USE sbc_oce        ! surface boundary condition: ocean
    USE phycst         ! physical constants
@@ -32,5 +29,4 @@
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
    USE lib_mpp        ! MPP manager
-   USE wrk_nemo       ! work arrays
    USE prtctl         ! Print control
    USE in_out_manager ! I/O manager
@@ -149,4 +145,5 @@
       REAL(wp) ::   prod, buoy, diss, zdiss, sm         !   -      -
       REAL(wp) ::   gh, gm, shr, dif, zsqen, zavt, zavm !   -      -
+      REAL(wp) ::   zmsku, zmskv                        !   -      -
       REAL(wp), DIMENSION(jpi,jpj)     ::   zdep
       REAL(wp), DIMENSION(jpi,jpj)     ::   zkar
@@ -158,7 +155,5 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwall_psi   ! Wall function use in the wb case (ln_sigpsi)
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   psi         ! psi at time now
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z_elem_a    ! element of the first  matrix diagonal
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z_elem_b    ! element of the second matrix diagonal
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z_elem_c    ! element of the third  matrix diagonal
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zd_lw, zd_up, zdiag   ! lower, upper  and diagonal of the matrix
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zstt, zstm  ! stability function on tracer and momentum
       !!--------------------------------------------------------------------
@@ -171,6 +166,4 @@
       ustar2_top (:,:) = 0._wp   ;   zwall_psi(:,:,:) = 0._wp
       ustar2_bot (:,:) = 0._wp
-      
-
 
       ! Compute surface, top and bottom friction at T-points
@@ -181,34 +174,23 @@
             ustar2_surf(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1)
             !   
-!!gm old
-            ! bottom friction (explicit before friction)        
-            ! Note that we chose here not to bound the friction as in dynbfr)   
-            ztx2 = (  bfrua(ji,jj)  * ub(ji,jj,mbku(ji,jj)) + bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj))  )   &         
-               & * ( 1._wp - 0.5_wp * umask(ji,jj,1) * umask(ji-1,jj,1)  )      
-            zty2 = (  bfrva(ji,jj)  * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1))  )   &         
-               & * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1)  )      
-            ustar2_bot(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)         
-         END DO         
-      END DO    
-!!gm new
 !!gm Rq we may add here r_ke0(_top/_bot) ?  ==>> think about that...
-!          ! bottom friction (explicit before friction)
-!          zmsku = ( 2._wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) )
-!          zmskv = ( 2._wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )     ! (CAUTION: CdU<0)
-!          ustar2_bot(ji,jj) = - rCdU_bot(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mbkt(ji,jj))+ub(ji-1,jj,mbkt(ji,jj)) ) )**2
-!             &                                         + ( zmskv*( vb(ji,jj,mbkt(ji,jj))+vb(ji,jj-1,mbkt(ji,jj)) ) )**2  )
-!       END DO
-!    END DO
-!    IF( ln_isfcav ) THEN       !top friction
-!       DO jj = 2, jpjm1
-!          DO ji = fs_2, fs_jpim1   ! vector opt.
-!             zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) )
-!             zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) )     ! (CAUTION: CdU<0)
-!             ustar2_top(ji,jj) = - rCdU_top(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mikt(ji,jj))+ub(ji-1,jj,mikt(ji,jj)) ) )**2
-!                &                                         + ( zmskv*( vb(ji,jj,mikt(ji,jj))+vb(ji,jj-1,mikt(ji,jj)) ) )**2  )
-!          END DO
-!       END DO
-!    ENDIF
-!!gm
+          ! bottom friction (explicit before friction)
+          zmsku = ( 2._wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) )
+          zmskv = ( 2._wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )     ! (CAUTION: CdU<0)
+          ustar2_bot(ji,jj) = - rCdU_bot(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mbkt(ji,jj))+ub(ji-1,jj,mbkt(ji,jj)) ) )**2  &
+             &                                         + ( zmskv*( vb(ji,jj,mbkt(ji,jj))+vb(ji,jj-1,mbkt(ji,jj)) ) )**2  )
+         END DO
+      END DO
+      IF( ln_isfcav ) THEN       !top friction
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) )
+               zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) )     ! (CAUTION: CdU<0)
+               ustar2_top(ji,jj) = - rCdU_top(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mikt(ji,jj))+ub(ji-1,jj,mikt(ji,jj)) ) )**2  &
+                  &                                         + ( zmskv*( vb(ji,jj,mikt(ji,jj))+vb(ji,jj-1,mikt(ji,jj)) ) )**2  )
+            END DO
+         END DO
+      ENDIF
+   
       SELECT CASE ( nn_z0_met )      !==  Set surface roughness length  ==!
       CASE ( 0 )                          ! Constant roughness          
@@ -261,5 +243,5 @@
       ! The surface boundary condition are set after
       ! The bottom boundary condition are also set after. In standard e(bottom)=0.
-      ! z_elem_b : diagonal z_elem_c : upper diagonal z_elem_a : lower diagonal
+      ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
       ! Warning : after this step, en : right hand side of the matrix
 
@@ -296,9 +278,9 @@
                zcof = rfact_tke * tmask(ji,jj,jk)
                !                                               ! lower diagonal
-               z_elem_a(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) )
+               zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) )
                !                                               ! upper diagonal
-               z_elem_c(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) ) / ( e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk) )
+               zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) ) / ( e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk) )
                !                                               ! diagonal
-               z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk)  + rdt * zdiss * wmask(ji,jj,jk) 
+               zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk)  + rdt * zdiss * wmask(ji,jj,jk) 
                !                                               ! right hand side in en
                en(ji,jj,jk) = en(ji,jj,jk) + rdt * zesh2 * wmask(ji,jj,jk)
@@ -307,5 +289,5 @@
       END DO
       !
-      z_elem_b(:,:,jpk) = 1._wp
+      zdiag(:,:,jpk) = 1._wp
       !
       ! Set surface condition on zwall_psi (1 at the bottom)
@@ -318,39 +300,36 @@
       SELECT CASE ( nn_bc_surf )
       !
-      CASE ( 0 )             ! Dirichlet case
+      CASE ( 0 )             ! Dirichlet boundary condition (set e at k=1 & 2) 
       ! First level
-      en(:,:,1) = rc02r * ustar2_surf(:,:) * (1._wp + rsbc_tke1)**r2_3
-      en(:,:,1) = MAX(en(:,:,1), rn_emin) 
-      z_elem_a(:,:,1) = en(:,:,1)
-      z_elem_c(:,:,1) = 0._wp
-      z_elem_b(:,:,1) = 1._wp
+      en   (:,:,1) = MAX(  rn_emin , rc02r * ustar2_surf(:,:) * (1._wp + rsbc_tke1)**r2_3  )
+      zd_lw(:,:,1) = en(:,:,1)
+      zd_up(:,:,1) = 0._wp
+      zdiag(:,:,1) = 1._wp
       ! 
       ! One level below
-      en(:,:,2) = rc02r * ustar2_surf(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+gdepw_n(:,:,2)) &
-         &               / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp)
-      en(:,:,2) = MAX(en(:,:,2), rn_emin )
-      z_elem_a(:,:,2) = 0._wp 
-      z_elem_c(:,:,2) = 0._wp
-      z_elem_b(:,:,2) = 1._wp
-      !
-      !
-      CASE ( 1 )             ! Neumann boundary condition on d(e)/dz
+      en   (:,:,2) =  MAX(  rc02r * ustar2_surf(:,:) * (  1._wp + rsbc_tke1 * ((zhsro(:,:)+gdepw_n(:,:,2))   &
+         &                 / zhsro(:,:) )**(1.5_wp*ra_sf)  )**(2._wp/3._wp)                      , rn_emin   )
+      zd_lw(:,:,2) = 0._wp 
+      zd_up(:,:,2) = 0._wp
+      zdiag(:,:,2) = 1._wp
+      !
+      !
+      CASE ( 1 )             ! Neumann boundary condition (set d(e)/dz)
       !
       ! Dirichlet conditions at k=1
-      en(:,:,1)       = rc02r * ustar2_surf(:,:) * (1._wp + rsbc_tke1)**r2_3
-      en(:,:,1)       = MAX(en(:,:,1), rn_emin)      
-      z_elem_a(:,:,1) = en(:,:,1)
-      z_elem_c(:,:,1) = 0._wp
-      z_elem_b(:,:,1) = 1._wp
+      en   (:,:,1) = MAX(  rc02r * ustar2_surf(:,:) * (1._wp + rsbc_tke1)**r2_3 , rn_emin  )
+      zd_lw(:,:,1) = en(:,:,1)
+      zd_up(:,:,1) = 0._wp
+      zdiag(:,:,1) = 1._wp
       !
       ! at k=2, set de/dz=Fw
       !cbr
-      z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
-      z_elem_a(:,:,2) = 0._wp
-      zkar(:,:)       = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*gdept_n(:,:,1)/zhsro(:,:)) ))
-      zflxs(:,:)      = rsbc_tke2 * ustar2_surf(:,:)**1.5_wp * zkar(:,:) &
-          &                       * ((zhsro(:,:)+gdept_n(:,:,1)) / zhsro(:,:) )**(1.5_wp*ra_sf)
-
-      en(:,:,2) = en(:,:,2) + zflxs(:,:)/e3w_n(:,:,2)
+      zdiag(:,:,2) = zdiag(:,:,2) +  zd_lw(:,:,2) ! Remove zd_lw from zdiag
+      zd_lw(:,:,2) = 0._wp
+      zkar (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*gdept_n(:,:,1)/zhsro(:,:)) ))
+      zflxs(:,:)   = rsbc_tke2 * ustar2_surf(:,:)**1.5_wp * zkar(:,:) &
+          &                    * (  ( zhsro(:,:)+gdept_n(:,:,1) ) / zhsro(:,:)  )**(1.5_wp*ra_sf)
+!!gm why not   :                        * ( 1._wp + gdept_n(:,:,1) / zhsro(:,:) )**(1.5_wp*ra_sf)
+      en(:,:,2) = en(:,:,2) + zflxs(:,:) / e3w_n(:,:,2)
       !
       !
@@ -377,29 +356,29 @@
                ! Dirichlet condition applied at: 
                !     Bottom level (ibot)      &      Just above it (ibotm1)   
-               z_elem_a(ji,jj,ibot) = 0._wp   ;   z_elem_a(ji,jj,ibotm1) = 0._wp
-               z_elem_c(ji,jj,ibot) = 0._wp   ;   z_elem_c(ji,jj,ibotm1) = 0._wp
-               z_elem_b(ji,jj,ibot) = 1._wp   ;   z_elem_b(ji,jj,ibotm1) = 1._wp
-               en      (ji,jj,ibot) = z_en    ;   en      (ji,jj,ibotm1) = z_en
-            END DO
-         END DO
-!!gm new
-!         IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
-!            DO jj = 2, jpjm1
-!               DO ji = fs_2, fs_jpim1   ! vector opt.
-!                  itop   = mikt(ji,jj)       ! k   top w-point
-!                  itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one
-!                  !                                                ! mask at the ocean surface points
-!                  z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) )
-!                  !
-!                  ! Dirichlet condition applied at: 
-!                  !     top level (itop)         &      Just below it (itopp1)   
-!                  z_elem_a(ji,jj,itop) = 0._wp   ;   z_elem_a(ji,jj,ipotm1) = 0._wp
-!                  z_elem_c(ji,jj,itop) = 0._wp   ;   z_elem_c(ji,jj,itopp1) = 0._wp
-!                  z_elem_b(ji,jj,itop) = 1._wp   ;   z_elem_b(ji,jj,itopp1) = 1._wp
-!                  en      (ji,jj,itop) = zen     ;   en      (ji,jj,itopp1) = z_en
-!               END DO
-!            END DO
-!         ENDIF
-!!gm
+               zd_lw(ji,jj,ibot) = 0._wp   ;   zd_lw(ji,jj,ibotm1) = 0._wp
+               zd_up(ji,jj,ibot) = 0._wp   ;   zd_up(ji,jj,ibotm1) = 0._wp
+               zdiag(ji,jj,ibot) = 1._wp   ;   zdiag(ji,jj,ibotm1) = 1._wp
+               en   (ji,jj,ibot) = z_en    ;   en   (ji,jj,ibotm1) = z_en
+            END DO
+         END DO
+         !
+         IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  itop   = mikt(ji,jj)       ! k   top w-point
+                  itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one
+                  !                                                ! mask at the ocean surface points
+                  z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) )
+                  !
+ !!gm TO BE VERIFIED !!!
+                  ! Dirichlet condition applied at: 
+                  !     top level (itop)         &      Just below it (itopp1)   
+                  zd_lw(ji,jj,itop) = 0._wp   ;   zd_lw(ji,jj,itopp1) = 0._wp
+                  zd_up(ji,jj,itop) = 0._wp   ;   zd_up(ji,jj,itopp1) = 0._wp
+                  zdiag(ji,jj,itop) = 1._wp   ;   zdiag(ji,jj,itopp1) = 1._wp
+                  en   (ji,jj,itop) = z_en    ;   en   (ji,jj,itopp1) = z_en
+               END DO
+            END DO
+         ENDIF
          !
       CASE ( 1 )             ! Neumman boundary condition
@@ -415,28 +394,26 @@
                !     Bottom level (ibot)      &      Just above it (ibotm1)   
                !         Dirichlet            !         Neumann
-               z_elem_a(ji,jj,ibot) = 0._wp   !   ! Remove z_elem_c from z_elem_b
-               z_elem_b(ji,jj,ibot) = 1._wp   ;   z_elem_b(ji,jj,ibotm1) = z_elem_b(ji,jj,ibotm1) + z_elem_c(ji,jj,ibotm1)
-               z_elem_c(ji,jj,ibot) = 0._wp   ;   z_elem_c(ji,jj,ibotm1) = 0._wp
-            END DO
-         END DO
-!!gm new
-!         IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
-!            DO jj = 2, jpjm1
-!               DO ji = fs_2, fs_jpim1   ! vector opt.
-!                  itop   = mikt(ji,jj)       ! k   top w-point
-!                  itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one
-!                  !                                                ! mask at the ocean surface points
-!                  z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) )
-!                  !
-!                  ! Bottom level Dirichlet condition:
-!                  !     Bottom level (ibot)      &      Just above it (ibotm1)   
-!                  !         Dirichlet            !         Neumann
-!                  z_elem_a(ji,jj,itop) = 0._wp   !   ! Remove z_elem_c from z_elem_b
-!                  z_elem_b(ji,jj,itop) = 1._wp   ;   z_elem_b(ji,jj,itopp1) = z_elem_b(ji,jj,itopp1) + z_elem_c(ji,jj,itopp1)
-!                  z_elem_c(ji,jj,itop) = 0._wp   ;   z_elem_c(ji,jj,itopp1) = 0._wp
-!               END DO
-!            END DO
-!         ENDIF
-!!gm
+               zd_lw(ji,jj,ibot) = 0._wp   !   ! Remove zd_up from zdiag
+               zdiag(ji,jj,ibot) = 1._wp   ;   zdiag(ji,jj,ibotm1) = zdiag(ji,jj,ibotm1) + zd_up(ji,jj,ibotm1)
+               zd_up(ji,jj,ibot) = 0._wp   ;   zd_up(ji,jj,ibotm1) = 0._wp
+            END DO
+         END DO
+         IF( ln_isfcav) THEN     ! top boundary   (ocean cavity)
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  itop   = mikt(ji,jj)       ! k   top w-point
+                  itopp1 = mikt(ji,jj) + 1   ! k+1 1st w-point below the top one
+                  !                                                ! mask at the ocean surface points
+                  z_en = MAX( rc02r * ustar2_top(ji,jj), rn_emin ) * ( 1._wp - tmask(ji,jj,1) )
+                  !
+                  ! Bottom level Dirichlet condition:
+                  !     Bottom level (ibot)      &      Just above it (ibotm1)   
+                  !         Dirichlet            !         Neumann
+                  zd_lw(ji,jj,itop) = 0._wp   !   ! Remove zd_up from zdiag
+                  zdiag(ji,jj,itop) = 1._wp   ;   zdiag(ji,jj,itopp1) = zdiag(ji,jj,itopp1) + zd_up(ji,jj,itopp1)
+                  zd_up(ji,jj,itop) = 0._wp   ;   zd_up(ji,jj,itopp1) = 0._wp
+               END DO
+            END DO
+         ENDIF
          !
       END SELECT
@@ -448,5 +425,5 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1    ! vector opt.
-               z_elem_b(ji,jj,jk) = z_elem_b(ji,jj,jk) - z_elem_a(ji,jj,jk) * z_elem_c(ji,jj,jk-1) / z_elem_b(ji,jj,jk-1)
+               zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
             END DO
          END DO
@@ -455,5 +432,5 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1    ! vector opt.
-               z_elem_a(ji,jj,jk) = en(ji,jj,jk) - z_elem_a(ji,jj,jk) / z_elem_b(ji,jj,jk-1) * z_elem_a(ji,jj,jk-1)
+               zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) * zd_lw(ji,jj,jk-1)
             END DO
          END DO
@@ -462,5 +439,5 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1    ! vector opt.
-               en(ji,jj,jk) = ( z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) * en(ji,jj,jk+1) ) / z_elem_b(ji,jj,jk)
+               en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
             END DO
          END DO
@@ -519,5 +496,5 @@
       ! Resolution of a tridiagonal linear system by a "methode de chasse"
       ! computation from level 2 to jpkm1  (e(1) already computed and e(jpk)=0 ).
-      ! z_elem_b : diagonal z_elem_c : upper diagonal z_elem_a : lower diagonal
+      ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
       ! Warning : after this step, en : right hand side of the matrix
 
@@ -543,5 +520,5 @@
                diss = rpsi2 * zratio * zwall(ji,jj,jk) * eps(ji,jj,jk)
                !
-               zdir = 0.5_wp + SIGN( 0.5_wp, prod + buoy )   ! zdir =1(=0) if shear(ji,jj,jk)+buoy >0(<0)
+               zdir = 0.5_wp + SIGN( 0.5_wp, prod + buoy )     ! zdir =1(=0) if shear(ji,jj,jk)+buoy >0(<0)
                !
                zesh2 = zdir * ( prod + buoy )          + (1._wp - zdir ) * prod                        ! production term
@@ -551,9 +528,9 @@
                zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk)
                !                                               ! lower diagonal
-               z_elem_a(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) )
+               zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) ) / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk) )
                !                                               ! upper diagonal
-               z_elem_c(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) ) / ( e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk) )
+               zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) ) / ( e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk) )
                !                                               ! diagonal
-               z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk)
+               zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk)
                !                                               ! right hand side in psi
                psi(ji,jj,jk) = psi(ji,jj,jk) + rdt * zesh2 * wmask(ji,jj,jk)
@@ -562,5 +539,5 @@
       END DO
       !
-      z_elem_b(:,:,jpk) = 1._wp
+      zdiag(:,:,jpk) = 1._wp
 
       ! Surface boundary condition on psi
@@ -574,7 +551,7 @@
          zdep    (:,:)   = zhsro(:,:) * rl_sf ! Cosmetic
          psi     (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
-         z_elem_a(:,:,1) = psi(:,:,1)
-         z_elem_c(:,:,1) = 0._wp
-         z_elem_b(:,:,1) = 1._wp
+         zd_lw(:,:,1) = psi(:,:,1)
+         zd_up(:,:,1) = 0._wp
+         zdiag(:,:,1) = 1._wp
          !
          ! One level below
@@ -582,7 +559,7 @@
          zdep    (:,:)   = (zhsro(:,:) + gdepw_n(:,:,2)) * zkar(:,:)
          psi     (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
-         z_elem_a(:,:,2) = 0._wp
-         z_elem_c(:,:,2) = 0._wp
-         z_elem_b(:,:,2) = 1._wp
+         zd_lw(:,:,2) = 0._wp
+         zd_up(:,:,2) = 0._wp
+         zdiag(:,:,2) = 1._wp
          ! 
       CASE ( 1 )             ! Neumann boundary condition on d(psi)/dz
@@ -591,11 +568,11 @@
          zdep    (:,:)   = zhsro(:,:) * rl_sf
          psi     (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
-         z_elem_a(:,:,1) = psi(:,:,1)
-         z_elem_c(:,:,1) = 0._wp
-         z_elem_b(:,:,1) = 1._wp
+         zd_lw(:,:,1) = psi(:,:,1)
+         zd_up(:,:,1) = 0._wp
+         zdiag(:,:,1) = 1._wp
          !
          ! Neumann condition at k=2
-         z_elem_b(:,:,2) = z_elem_b(:,:,2) +  z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b
-         z_elem_a(:,:,2) = 0._wp
+         zdiag(:,:,2) = zdiag(:,:,2) +  zd_lw(:,:,2) ! Remove zd_lw from zdiag
+         zd_lw(:,:,2) = 0._wp
          !
          ! Set psi vertical flux at the surface:
@@ -613,8 +590,11 @@
       ! --------------------------------
       !
-      SELECT CASE ( nn_bc_bot )
+!!gm should be done for ISF (top boundary cond.)
+!!gm so, totally new staff needed      ===>>> think about that !
+!
+      SELECT CASE ( nn_bc_bot )     ! bottom boundary
       !
       CASE ( 0 )             ! Dirichlet 
-         !                      ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = vkarmn * rn_bfrz0
+         !                      ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = vkarmn * r_z0_bot
          !                      ! Balance between the production and the dissipation terms
          DO jj = 2, jpjm1
@@ -622,16 +602,16 @@
                ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point
                ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1
-               zdep(ji,jj) = vkarmn * rn_bfrz0
+               zdep(ji,jj) = vkarmn * r_z0_bot
                psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn
-               z_elem_a(ji,jj,ibot) = 0._wp
-               z_elem_c(ji,jj,ibot) = 0._wp
-               z_elem_b(ji,jj,ibot) = 1._wp
+               zd_lw(ji,jj,ibot) = 0._wp
+               zd_up(ji,jj,ibot) = 0._wp
+               zdiag(ji,jj,ibot) = 1._wp
                !
                ! Just above last level, Dirichlet condition again (GOTM like)
-               zdep(ji,jj) = vkarmn * ( rn_bfrz0 + e3t_n(ji,jj,ibotm1) )
+               zdep(ji,jj) = vkarmn * ( r_z0_bot + e3t_n(ji,jj,ibotm1) )
                psi (ji,jj,ibotm1) = rc0**rpp * en(ji,jj,ibot  )**rmm * zdep(ji,jj)**rnn
-               z_elem_a(ji,jj,ibotm1) = 0._wp
-               z_elem_c(ji,jj,ibotm1) = 0._wp
-               z_elem_b(ji,jj,ibotm1) = 1._wp
+               zd_lw(ji,jj,ibotm1) = 0._wp
+               zd_up(ji,jj,ibotm1) = 0._wp
+               zdiag(ji,jj,ibotm1) = 1._wp
             END DO
          END DO
@@ -645,17 +625,17 @@
                !
                ! Bottom level Dirichlet condition:
-               zdep(ji,jj) = vkarmn * rn_bfrz0
+               zdep(ji,jj) = vkarmn * r_z0_bot
                psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn
                !
-               z_elem_a(ji,jj,ibot) = 0._wp
-               z_elem_c(ji,jj,ibot) = 0._wp
-               z_elem_b(ji,jj,ibot) = 1._wp
+               zd_lw(ji,jj,ibot) = 0._wp
+               zd_up(ji,jj,ibot) = 0._wp
+               zdiag(ji,jj,ibot) = 1._wp
                !
                ! Just above last level: Neumann condition with flux injection
-               z_elem_b(ji,jj,ibotm1) = z_elem_b(ji,jj,ibotm1) + z_elem_c(ji,jj,ibotm1) ! Remove z_elem_c from z_elem_b
-               z_elem_c(ji,jj,ibotm1) = 0.
+               zdiag(ji,jj,ibotm1) = zdiag(ji,jj,ibotm1) + zd_up(ji,jj,ibotm1) ! Remove zd_up from zdiag
+               zd_up(ji,jj,ibotm1) = 0.
                !
                ! Set psi vertical flux at the bottom:
-               zdep(ji,jj) = rn_bfrz0 + 0.5_wp*e3t_n(ji,jj,ibotm1)
+               zdep(ji,jj) = r_z0_bot + 0.5_wp*e3t_n(ji,jj,ibotm1)
                zflxb = rsbc_psi2 * ( p_avm(ji,jj,ibot) + p_avm(ji,jj,ibotm1) )   &
                   &  * (0.5_wp*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1._wp)
@@ -672,5 +652,5 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1    ! vector opt.
-               z_elem_b(ji,jj,jk) = z_elem_b(ji,jj,jk) - z_elem_a(ji,jj,jk) * z_elem_c(ji,jj,jk-1) / z_elem_b(ji,jj,jk-1)
+               zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
             END DO
          END DO
@@ -679,5 +659,5 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1    ! vector opt.
-               z_elem_a(ji,jj,jk) = psi(ji,jj,jk) - z_elem_a(ji,jj,jk) / z_elem_b(ji,jj,jk-1) * z_elem_a(ji,jj,jk-1)
+               zd_lw(ji,jj,jk) = psi(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) * zd_lw(ji,jj,jk-1)
             END DO
          END DO
@@ -686,5 +666,5 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1    ! vector opt.
-               psi(ji,jj,jk) = ( z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) * psi(ji,jj,jk+1) ) / z_elem_b(ji,jj,jk)
+               psi(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * psi(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
             END DO
          END DO
@@ -814,5 +794,4 @@
       zstm(:,:,1) = zstm(:,:,2)
 
-!!gm should be done for ISF (top boundary cond.)
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -820,4 +799,6 @@
          END DO
       END DO
+!!gm should be done for ISF (top boundary cond.)
+!!gm so, totally new staff needed!!gm
 
       ! Compute diffusivities/viscosities
@@ -904,12 +885,8 @@
          WRITE(numout,*) '      Type of closure                               nn_clos        = ', nn_clos
          WRITE(numout,*) '      Surface roughness (m)                         rn_hsro        = ', rn_hsro
-!!gm old
-         WRITE(numout,*) '      top    roughness (m) (nambfr namelist)        rn_tfrz0       = ', rn_tfrz0
-         WRITE(numout,*) '      Bottom roughness (m) (nambfr namelist)        rn_bfrz0       = ', rn_bfrz0
-!!gm new
-!         WRITE(numout,*) '   Namelist namdrg_top/_bot used values:'
-!         WRITE(numout,*) '      top    ocean cavity roughness (m)             rn_z0(_top)   = ', r_z0_top
-!         WRITE(numout,*) '      Bottom seafloor     roughness (m)             rn_z0(_bot)   = ', r_z0_bot
-!!gm
+         WRITE(numout,*)
+         WRITE(numout,*) '   Namelist namdrg_top/_bot:   used values:'
+         WRITE(numout,*) '      top    ocean cavity roughness (m)             rn_z0(_top)   = ', r_z0_top
+         WRITE(numout,*) '      Bottom seafloor     roughness (m)             rn_z0(_bot)   = ', r_z0_bot
          WRITE(numout,*)
       ENDIF
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90	(revision 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90	(revision 8143)
@@ -11,25 +11,24 @@
    !!   zdf_mxl      : Compute the turbocline and mixed layer depths.
    !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers variables
-   USE dom_oce         ! ocean space and time domain variables
-   USE trc_oce, ONLY: l_offline         ! ocean space and time domain variables
-   USE zdf_oce         ! ocean vertical physics
-   USE in_out_manager  ! I/O manager
-   USE prtctl          ! Print control
-   USE phycst          ! physical constants
-   USE iom             ! I/O library
-   USE lib_mpp         ! MPP library
-   USE wrk_nemo        ! work arrays
-   USE timing          ! Timing
+   USE oce            ! ocean dynamics and tracers variables
+   USE dom_oce        ! ocean space and time domain variables
+   USE trc_oce  , ONLY: l_offline         ! ocean space and time domain variables
+   USE zdf_oce        ! ocean vertical physics
+   USE in_out_manager ! I/O manager
+   USE prtctl         ! Print control
+   USE phycst         ! physical constants
+   USE iom            ! I/O library
+   USE lib_mpp        ! MPP library
+   USE timing         ! Timing
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   zdf_mxl       ! called by step.F90
+   PUBLIC   zdf_mxl   ! called by zdfphy.F90
 
-   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   nmln    !: number of level in the mixed layer (used by TOP)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmld    !: mixing layer depth (turbocline)      [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlp    !: mixed layer depth  (rho=rho0+zdcrit) [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlpt   !: depth of the last T-point inside the mixed layer [m]
+   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   nmln    !: number of level in the mixed layer (used by LDF, ZDF, TRD, TOP)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmld    !: mixing layer depth (turbocline)      [m]   (used by TOP)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlp    !: mixed layer depth  (rho=rho0+zdcrit) [m]   (used by LDF)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlpt   !: depth of the last T-point inside the mixed layer [m] (used by LDF)
 
    REAL(wp), PUBLIC ::   rho_c = 0.01_wp    !: density criterion for mixed layer depth
@@ -37,5 +36,5 @@
 
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
+   !! NEMO/OPA 4.0 , NEMO Consortium (2017)
    !! $Id$ 
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -80,11 +79,9 @@
       INTEGER  ::   iikn, iiki, ikt ! local integer
       REAL(wp) ::   zN2_c           ! local scalar
-      INTEGER, POINTER, DIMENSION(:,:) ::   imld   ! 2D workspace
+      INTEGER, DIMENSION(jpi,jpj) ::   imld   ! 2D workspace
       !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('zdf_mxl')
       !
-      CALL wrk_alloc( jpi,jpj, imld )
-
       IF( kt == nit000 ) THEN
          IF(lwp) WRITE(numout,*)
@@ -144,6 +141,4 @@
       IF(ln_ctl)   CALL prt_ctl( tab2d_1=REAL(nmln,wp), clinfo1=' nmln : ', tab2d_2=hmlp, clinfo2=' hmlp : ', ovlap=1 )
       !
-      CALL wrk_dealloc( jpi,jpj, imld )
-      !
       IF( nn_timing == 1 )  CALL timing_stop('zdf_mxl')
       !
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 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfphy.F90	(revision 8143)
@@ -14,7 +14,4 @@
    USE zdf_oce        ! vertical physics: shared variables         
    USE zdfdrg         ! vertical physics: top/bottom drag coef.
-!!gm old
-   USE zdfbfr         ! vertical physics: bottom friction
-!!gm
    USE zdfsh2         ! vertical physics: shear production term of TKE
    USE zdfric         ! vertical physics: RIChardson dependent vertical mixing   
@@ -195,7 +192,4 @@
       !                          !== top/bottom friction  ==!
       CALL zdf_drg_init
-!!gm old
-      CALL zdf_bfr_init
-!!gm
       !
       !                          !== time-stepping  ==!
@@ -224,20 +218,16 @@
       !! ---------------------------------------------------------------------
       !
-!!gm old
-      CALL zdf_bfr( kt )                        !* bottom friction (if quadratic)
-!!gm
-      !
-!      IF( l_zdfdrg ) THEN     !==  update top/bottom drag  ==!   (non-linear cases)
-!         !
-!         !                       !* bottom drag
-!         CALL zdf_drg( kt, mbkt    , r_Cdmin_bot, r_Cdmax_bot,   &   ! <<== in 
-!            &              r_z0_bot,   r_ke0_bot,    rCd0_bot,   &
-!            &                                        rCdU_bot  )     ! ==>> out : bottom drag [m/s]
-!         IF( ln_isfcav ) THEN    !* top drag   (ocean cavities)
-!            CALL zdf_drg( kt, mikt    , r_Cdmin_top, r_Cdmax_top,   &   ! <<== in 
-!               &              r_z0_top,   r_ke0_top,    rCd0_top,   &
-!               &                                        rCdU_top  )     ! ==>> out : bottom drag [m/s]
-!         ENDIF
-!      ENDIF
+      IF( l_zdfdrg ) THEN     !==  update top/bottom drag  ==!   (non-linear cases)
+         !
+         !                       !* bottom drag
+         CALL zdf_drg( kt, mbkt    , r_Cdmin_bot, r_Cdmax_bot,   &   ! <<== in 
+            &              r_z0_bot,   r_ke0_bot,    rCd0_bot,   &
+            &                                        rCdU_bot  )     ! ==>> out : bottom drag [m/s]
+         IF( ln_isfcav ) THEN    !* top drag   (ocean cavities)
+            CALL zdf_drg( kt, mikt    , r_Cdmin_top, r_Cdmax_top,   &   ! <<== in 
+               &              r_z0_top,   r_ke0_top,    rCd0_top,   &
+               &                                        rCdU_top  )     ! ==>> out : bottom drag [m/s]
+         ENDIF
+      ENDIF
       !
       !                       !==  Kz from chosen turbulent closure  ==!   (avm_k, avt_k)
@@ -290,6 +280,9 @@
       CALL lbc_lnk( avs  , 'W', 1. )                  !!gm  To be tested
       !
-
-
+      IF( l_zdfdrg ) THEN     ! drag  have been updated (non-linear cases)
+         IF( ln_isfcav )   CALL lbc_lnk( rCdU_top, 'T', 1. )   ! top    drag
+                           CALL lbc_lnk( rCdU_bot, 'T', 1. )   ! bottom drag
+      ENDIF
+      !
       CALL zdf_mxl( kt )                        !* mixed layer depth, and level
 
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 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90	(revision 8143)
@@ -61,5 +61,5 @@
    SUBROUTINE zdf_ric_init
       !!----------------------------------------------------------------------
-      !!                 ***  ROUTINE zdfbfr_init  ***
+      !!                 ***  ROUTINE zdf_ric_init  ***
       !!                    
       !! ** Purpose :   Initialization of the vertical eddy diffusivity and
@@ -90,6 +90,6 @@
       IF(lwp) THEN                   ! Control print
          WRITE(numout,*)
-         WRITE(numout,*) 'zdf_ric : Ri depend vertical mixing scheme'
-         WRITE(numout,*) '~~~~~~~'
+         WRITE(numout,*) 'zdf_ric_init : Ri depend vertical mixing scheme'
+         WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namzdf_ric : set Kz=F(Ri) parameters'
          WRITE(numout,*) '      maximum vertical viscosity        rn_avmri  = ', rn_avmri
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 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 8143)
@@ -28,4 +28,5 @@
    !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability
    !!            4.0  !  2017-04  (G. Madec)  remove CPP ddm key & avm at t-point only 
+   !!             -   !  2017-05  (G. Madec)  add top/bottom friction as boundary condition (ln_drg)
    !!----------------------------------------------------------------------
 
@@ -43,7 +44,5 @@
    USE sbc_oce        ! surface boundary condition: ocean
    USE zdf_oce        ! vertical physics: ocean variables
-!!gm new
    USE zdfdrg         ! vertical physics: top/bottom drag coef.
-!!gm
    USE zdfmxl         ! vertical physics: mixed layer
 #if defined key_agrif
@@ -57,5 +56,4 @@
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
    USE prtctl         ! Print control
-   USE wrk_nemo       ! work arrays
    USE timing         ! Timing
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
@@ -79,9 +77,10 @@
    REAL(wp) ::   rn_emin0  ! surface minimum value of tke   [m2/s2]
    REAL(wp) ::   rn_bshear ! background shear (>0) currently a numerical threshold (do not change it)
+   LOGICAL  ::   ln_drg    ! top/bottom friction forcing flag 
    INTEGER  ::   nn_etau   ! type of depth penetration of surface tke (=0/1/2/3)
-   INTEGER  ::   nn_htau   ! type of tke profile of penetration (=0/1)
-   REAL(wp) ::   rn_efr    ! fraction of TKE surface value which penetrates in the ocean
+   INTEGER  ::      nn_htau   ! type of tke profile of penetration (=0/1)
+   REAL(wp) ::      rn_efr    ! fraction of TKE surface value which penetrates in the ocean
    LOGICAL  ::   ln_lc     ! Langmuir cells (LC) as a source term of TKE or not
-   REAL(wp) ::   rn_lc     ! coef to compute vertical velocity of Langmuir cells
+   REAL(wp) ::      rn_lc     ! coef to compute vertical velocity of Langmuir cells
 
    REAL(wp) ::   ri_cri    ! critic Richardson number (deduced from rn_ediff and rn_ediss values)
@@ -204,13 +203,13 @@
       !
       INTEGER ::   ji, jj, jk              ! dummy loop arguments
-!!bfr      REAL(wp) ::   zebot, zmshu, zmskv      ! local scalars
-      REAL(wp) ::   zrhoa  = 1.22            ! Air density kg/m3
-      REAL(wp) ::   zcdrag = 1.5e-3          ! drag coefficient
-      REAL(wp) ::   zbbrau, zri              ! local scalars
-      REAL(wp) ::   zfact1, zfact2, zfact3   !   -         -
-      REAL(wp) ::   ztx2  , zty2  , zcof     !   -         -
-      REAL(wp) ::   ztau  , zdif             !   -         -
-      REAL(wp) ::   zus   , zwlc  , zind     !   -         -
-      REAL(wp) ::   zzd_up, zzd_lw           !   -         -
+      REAL(wp) ::   zetop, zebot, zmsku, zmskv ! local scalars
+      REAL(wp) ::   zrhoa  = 1.22              ! Air density kg/m3
+      REAL(wp) ::   zcdrag = 1.5e-3            ! drag coefficient
+      REAL(wp) ::   zbbrau, zri                ! local scalars
+      REAL(wp) ::   zfact1, zfact2, zfact3     !   -         -
+      REAL(wp) ::   ztx2  , zty2  , zcof       !   -         -
+      REAL(wp) ::   ztau  , zdif               !   -         -
+      REAL(wp) ::   zus   , zwlc  , zind       !   -         -
+      REAL(wp) ::   zzd_up, zzd_lw             !   -         -
       INTEGER , DIMENSION(jpi,jpj)     ::   imlc
       REAL(wp), DIMENSION(jpi,jpj)     ::   zhlc
@@ -227,6 +226,12 @@
       !
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      !                     !  Surface boundary condition on tke
-      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                     !  Surface/top/bottom boundary condition on tke
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      
+      DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0)
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
+         END DO
+      END DO
       IF ( ln_isfcav ) THEN
          DO jj = 2, jpjm1            ! en(mikt(ji,jj))   = rn_emin
@@ -235,64 +240,40 @@
             END DO
          END DO
-      END IF
-      DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0)
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
-         END DO
-      END DO
+      ENDIF
       
-!!bfr   - start commented area
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !                     !  Bottom boundary condition on tke
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      ! Tests to date have found the bottom boundary condition on tke to have very little effect.
-      ! The condition is coded here for completion but commented out until there is proof that the
-      ! computational cost is justified
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      !                     en(bot)   = (rn_ebb0/rau0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin)
-!!gm old
-!!    DO jj = 2, jpjm1
-!!       DO ji = fs_2, fs_jpim1   ! vector opt.
-!!          ztx2 = bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) + &
-!!                 bfrua(ji  ,jj) * ub(ji  ,jj,mbku(ji  ,jj) )
-!!          zty2 = bfrva(ji,jj  ) * vb(ji,jj  ,mbkv(ji,jj  )) + &
-!!                 bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1) )
-!!          zebot = 0.001875_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 )   !  where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000.
-!!          en(ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * ssmask(ji,jj)
-!!       END DO
-!!    END DO
-!!gm new
-!!
-!!    ====>>>>  add below an wet-only calculation of u and v at t-point like in zdfsh2:
-!!                   zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) )
-!!                   zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )
-!!
-!!
-!!    DO jj = 2, jpjm1
-!!       DO ji = fs_2, fs_jpim1   ! vector opt.
-!!          zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) )
-!!          zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )
-!!          !                             ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0)
-!!          zebot = - 0.001875_wp * rCdU_bot(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mbkt(ji,jj))+ub(ji-1,jj,mbkt(ji,jj)) ) )**2
-!!             &                                           + ( zmskv*( vb(ji,jj,mbkt(ji,jj))+vb(ji,jj-1,mbkt(ji,jj)) ) )**2  )
-!!          en(ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * ssmask(ji,jj)
-!!       END DO
-!!    END DO
-!!    IF( ln_isfcav ) THEN       !top friction
-!!       DO jj = 2, jpjm1
-!!          DO ji = fs_2, fs_jpim1   ! vector opt.
-!!             zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) )
-!!             zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) )
-!!             !                             ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000.  (CAUTION CdU<0)
-!!             zebot = - 0.001875_wp * rCdU_top(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mikt(ji,jj))+ub(ji-1,jj,mikt(ji,jj)) ) )**2
-!!                &                                           + ( zmskv*( vb(ji,jj,mikt(ji,jj))+vb(ji,jj-1,mikt(ji,jj)) ) )**2  )
-!!             en(ji,jj,mikt(ji,jj)+1) = MAX( zebot, rn_emin ) * (1._wp - tmask(ji,jj,1))   ! masked at ocean surface
-!!          END DO
-!!       END DO
-!!    ENDIF
-!!
-!!bfr   - end commented area
+      !   en(bot)   = (ebb0/rau0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin)
+      ! where ebb0 does not includes surface wave enhancement (i.e. ebb0=3.75)
+      ! Note that stress averaged is done using an wet-only calculation of u and v at t-point like in zdfsh2
+      !
+      IF( ln_drg ) THEN       !== friction used as top/bottom boundary condition on TKE
+         !
+         DO jj = 2, jpjm1           ! bottom friction
+            DO ji = fs_2, fs_jpim1     ! vector opt.
+               zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) )
+               zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )
+               !                       ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0)
+               zebot = - 0.001875_wp * rCdU_bot(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mbkt(ji,jj))+ub(ji-1,jj,mbkt(ji,jj)) ) )**2  &
+                  &                                           + ( zmskv*( vb(ji,jj,mbkt(ji,jj))+vb(ji,jj-1,mbkt(ji,jj)) ) )**2  )
+               en(ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * ssmask(ji,jj)
+            END DO
+         END DO
+         IF( ln_isfcav ) THEN       ! top friction
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) )
+                  zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) )
+                  !                             ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000.  (CAUTION CdU<0)
+                  zetop = - 0.001875_wp * rCdU_top(ji,jj) * SQRT(  ( zmsku*( ub(ji,jj,mikt(ji,jj))+ub(ji-1,jj,mikt(ji,jj)) ) )**2  &
+                     &                                           + ( zmskv*( vb(ji,jj,mikt(ji,jj))+vb(ji,jj-1,mikt(ji,jj)) ) )**2  )
+                  en(ji,jj,mikt(ji,jj)) = MAX( zetop, rn_emin ) * (1._wp - tmask(ji,jj,1))   ! masked at ocean surface
+               END DO
+            END DO
+         ENDIF
+         !
+      ENDIF
       !
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
@@ -426,4 +407,5 @@
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 !!gm BUG : in the exp  remove the depth of ssh !!!
+!!gm       i.e. use gde3w in argument (pdepw)
       
       
@@ -678,5 +660,5 @@
       NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin  ,   &
          &                 rn_emin0, rn_bshear, nn_mxl , ln_mxl0  ,   &
-         &                 rn_mxl0 , nn_pdl   , ln_lc  , rn_lc    ,   &
+         &                 rn_mxl0 , nn_pdl   , ln_drg , ln_lc    , rn_lc    ,   &
          &                 nn_etau , nn_htau  , rn_efr   
       !!----------------------------------------------------------------------
@@ -703,16 +685,24 @@
          WRITE(numout,*) '      minimum value of tke                        rn_emin   = ', rn_emin
          WRITE(numout,*) '      surface minimum value of tke                rn_emin0  = ', rn_emin0
+         WRITE(numout,*) '      prandl number flag                          nn_pdl    = ', nn_pdl
          WRITE(numout,*) '      background shear (>0)                       rn_bshear = ', rn_bshear
          WRITE(numout,*) '      mixing length type                          nn_mxl    = ', nn_mxl
-         WRITE(numout,*) '      prandl number flag                          nn_pdl    = ', nn_pdl
-         WRITE(numout,*) '      surface mixing length = F(stress) or not    ln_mxl0   = ', ln_mxl0
-         WRITE(numout,*) '      surface  mixing length minimum value        rn_mxl0   = ', rn_mxl0
-         WRITE(numout,*) '      flag to take into acc.  Langmuir circ.      ln_lc     = ', ln_lc
-         WRITE(numout,*) '      coef to compute verticla velocity of LC     rn_lc     = ', rn_lc
+         WRITE(numout,*) '         surface mixing length = F(stress) or not    ln_mxl0   = ', ln_mxl0
+         WRITE(numout,*) '         surface  mixing length minimum value        rn_mxl0   = ', rn_mxl0
+         WRITE(numout,*) '      top/bottom friction forcing flag            ln_drg    = ', ln_drg
+         WRITE(numout,*) '      Langmuir cells parametrization              ln_lc     = ', ln_lc
+         WRITE(numout,*) '         coef to compute vertical velocity of LC     rn_lc  = ', rn_lc
          WRITE(numout,*) '      test param. to add tke induced by wind      nn_etau   = ', nn_etau
-         WRITE(numout,*) '      flag for computation of exp. tke profile    nn_htau   = ', nn_htau
-         WRITE(numout,*) '      fraction of en which pene. the thermocline  rn_efr    = ', rn_efr
+         WRITE(numout,*) '          type of tke penetration profile            nn_htau   = ', nn_htau
+         WRITE(numout,*) '          fraction of TKE that penetrates            rn_efr    = ', rn_efr
          WRITE(numout,*)
-         WRITE(numout,*) '      critical Richardson nb with your parameters  ri_cri = ', ri_cri
+         IF( ln_drg ) THEN
+            WRITE(numout,*) '   Namelist namdrg_top/_bot:   used values:'
+            WRITE(numout,*) '      top    ocean cavity roughness (m)          rn_z0(_top)= ', r_z0_top
+            WRITE(numout,*) '      Bottom seafloor     roughness (m)          rn_z0(_bot)= ', r_z0_bot
+         ENDIF
+         WRITE(numout,*)
+         WRITE(numout,*)
+         WRITE(numout,*) '   ==>> critical Richardson nb with your parameters  ri_cri = ', ri_cri
          WRITE(numout,*)
       ENDIF
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 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 8143)
@@ -129,5 +129,5 @@
 
       !  VERTICAL PHYSICS
-                         CALL zdf_phy( kstp )         ! vertical physics update (bfr, avt, avs, avm + MLD)
+                         CALL zdf_phy( kstp )         ! vertical physics update (top/bot drag, avt, avs, avm + MLD)
 
 
@@ -211,5 +211,6 @@
       ENDIF
 
-                         CALL dyn_bfr       ( kstp )  ! bottom friction
+      IF( .NOT.ln_drgimp)   CALL dyn_bfr       ( kstp )  ! bottom friction
+      
                          CALL dyn_zdf       ( kstp )  ! vertical diffusion
 
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 8093)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 8143)
@@ -7,61 +7,62 @@
    !!             3.7  !  2014-01  (G. Madec) LDF simplication 
    !!----------------------------------------------------------------------
-   USE oce              ! ocean dynamics and tracers variables
-   USE dom_oce          ! ocean space and time domain variables
-   USE zdf_oce          ! ocean vertical physics variables
+   USE oce             ! ocean dynamics and tracers variables
+   USE dom_oce         ! ocean space and time domain variables
+   USE zdf_oce         ! ocean vertical physics variables
+   USE zdfdrg  ,  ONLY : ln_drgimp   ! implicit top/bottom friction
 
-   USE daymod           ! calendar                         (day     routine)
+   USE daymod          ! calendar                         (day     routine)
 
-   USE sbc_oce          ! surface boundary condition: ocean
-   USE sbcmod           ! surface boundary condition       (sbc     routine)
-   USE sbcrnf           ! surface boundary condition: runoff variables
-   USE sbccpl           ! surface boundary condition: coupled formulation (call send at end of step)
-   USE sbcapr           ! surface boundary condition: atmospheric pressure
-   USE sbctide          ! Tide initialisation
-   USE sbcwave          ! Wave intialisation
+   USE sbc_oce         ! surface boundary condition: ocean
+   USE sbcmod          ! surface boundary condition       (sbc     routine)
+   USE sbcrnf          ! surface boundary condition: runoff variables
+   USE sbccpl          ! surface boundary condition: coupled formulation (call send at end of step)
+   USE sbcapr          ! surface boundary condition: atmospheric pressure
+   USE sbctide         ! Tide initialisation
+   USE sbcwave         ! Wave intialisation
 
-   USE traqsr           ! solar radiation penetration      (tra_qsr routine)
-   USE trasbc           ! surface boundary condition       (tra_sbc routine)
-   USE trabbc           ! bottom boundary condition        (tra_bbc routine)
-   USE trabbl           ! bottom boundary layer            (tra_bbl routine)
-   USE tradmp           ! internal damping                 (tra_dmp routine)
-   USE traadv           ! advection scheme control     (tra_adv_ctl routine)
-   USE traldf           ! lateral mixing                   (tra_ldf routine)
-   USE trazdf           ! vertical mixing                  (tra_zdf routine)
-   USE tranxt           ! time-stepping                    (tra_nxt routine)
-   USE tranpc           ! non-penetrative convection       (tra_npc routine)
+   USE traqsr          ! solar radiation penetration      (tra_qsr routine)
+   USE trasbc          ! surface boundary condition       (tra_sbc routine)
+   USE trabbc          ! bottom boundary condition        (tra_bbc routine)
+   USE trabbl          ! bottom boundary layer            (tra_bbl routine)
+   USE tradmp          ! internal damping                 (tra_dmp routine)
+   USE traadv          ! advection scheme control     (tra_adv_ctl routine)
+   USE traldf          ! lateral mixing                   (tra_ldf routine)
+   USE trazdf          ! vertical mixing                  (tra_zdf routine)
+   USE tranxt          ! time-stepping                    (tra_nxt routine)
+   USE tranpc          ! non-penetrative convection       (tra_npc routine)
 
-   USE eosbn2           ! equation of state                (eos_bn2 routine)
+   USE eosbn2          ! equation of state                (eos_bn2 routine)
 
-   USE divhor           ! horizontal divergence            (div_hor routine)
-   USE dynadv           ! advection                        (dyn_adv routine)
-   USE dynbfr           ! Bottom friction terms            (dyn_bfr routine)
-   USE dynvor           ! vorticity term                   (dyn_vor routine)
-   USE dynhpg           ! hydrostatic pressure grad.       (dyn_hpg routine)
-   USE dynldf           ! lateral momentum diffusion       (dyn_ldf routine)
-   USE dynzdf           ! vertical diffusion               (dyn_zdf routine)
-   USE dynspg           ! surface pressure gradient        (dyn_spg routine)
+   USE divhor          ! horizontal divergence            (div_hor routine)
+   USE dynadv          ! advection                        (dyn_adv routine)
+   USE dynbfr          ! Bottom friction terms            (dyn_bfr routine)
+   USE dynvor          ! vorticity term                   (dyn_vor routine)
+   USE dynhpg          ! hydrostatic pressure grad.       (dyn_hpg routine)
+   USE dynldf          ! lateral momentum diffusion       (dyn_ldf routine)
+   USE dynzdf          ! vertical diffusion               (dyn_zdf routine)
+   USE dynspg          ! surface pressure gradient        (dyn_spg routine)
 
-   USE dynnxt           ! time-stepping                    (dyn_nxt routine)
+   USE dynnxt          ! time-stepping                    (dyn_nxt routine)
 
-   USE stopar           ! Stochastic parametrization       (sto_par routine)
+   USE stopar          ! Stochastic parametrization       (sto_par routine)
    USE stopts 
 
-   USE bdy_oce    , ONLY: ln_bdy
-   USE bdydta           ! open boundary condition data     (bdy_dta routine)
-   USE bdytra           ! bdy cond. for tracers            (bdy_tra routine)
-   USE bdydyn3d         ! bdy cond. for baroclinic vel.  (bdy_dyn3d routine)
+   USE bdy_oce  , ONLY : ln_bdy
+   USE bdydta          ! open boundary condition data     (bdy_dta routine)
+   USE bdytra          ! bdy cond. for tracers            (bdy_tra routine)
+   USE bdydyn3d        ! bdy cond. for baroclinic vel.  (bdy_dyn3d routine)
 
-   USE sshwzv           ! vertical velocity and ssh        (ssh_nxt routine)
+   USE sshwzv          ! vertical velocity and ssh        (ssh_nxt routine)
    !                                                       (ssh_swp routine)
    !                                                       (wzv     routine)
-   USE domvvl           ! variable vertical scale factors  (dom_vvl_sf_nxt routine)
+   USE domvvl          ! variable vertical scale factors  (dom_vvl_sf_nxt routine)
    !                                                       (dom_vvl_sf_swp routine)
 
-   USE ldfslp           ! iso-neutral slopes               (ldf_slp routine)
-   USE ldfdyn           ! lateral eddy viscosity coef.     (ldf_dyn routine)
-   USE ldftra           ! lateral eddy diffusive coef.     (ldf_tra routine)
+   USE ldfslp          ! iso-neutral slopes               (ldf_slp routine)
+   USE ldfdyn          ! lateral eddy viscosity coef.     (ldf_dyn routine)
+   USE ldftra          ! lateral eddy diffusive coef.     (ldf_tra routine)
 
-   USE zdfphy         ! vertical physics manager      (zdf_phy_init routine)
+   USE zdfphy          ! vertical physics manager      (zdf_phy_init routine)
 
    USE step_diu        ! Time stepping for diurnal sst
@@ -70,35 +71,35 @@
    USE sbc_oce         ! surface fluxes  
    
-   USE zpshde           ! partial step: hor. derivative     (zps_hde routine)
+   USE zpshde          ! partial step: hor. derivative     (zps_hde routine)
 
-   USE diawri           ! Standard run outputs             (dia_wri routine)
-   USE diaptr           ! poleward transports              (dia_ptr routine)
-   USE diadct           ! sections transports              (dia_dct routine)
-   USE diaar5           ! AR5 diagnosics                   (dia_ar5 routine)
-   USE diahth           ! thermocline depth                (dia_hth routine)
-   USE diahsb           ! heat, salt and volume budgets    (dia_hsb routine)
+   USE diawri          ! Standard run outputs             (dia_wri routine)
+   USE diaptr          ! poleward transports              (dia_ptr routine)
+   USE diadct          ! sections transports              (dia_dct routine)
+   USE diaar5          ! AR5 diagnosics                   (dia_ar5 routine)
+   USE diahth          ! thermocline depth                (dia_hth routine)
+   USE diahsb          ! heat, salt and volume budgets    (dia_hsb routine)
    USE diaharm
    USE diacfl
-   USE flo_oce          ! floats variables
-   USE floats           ! floats computation               (flo_stp routine)
+   USE flo_oce         ! floats variables
+   USE floats          ! floats computation               (flo_stp routine)
 
-   USE crsfld           ! Standard output on coarse grid   (crs_fld routine)
+   USE crsfld          ! Standard output on coarse grid   (crs_fld routine)
 
-   USE asminc           ! assimilation increments      (tra_asm_inc routine)
+   USE asminc          ! assimilation increments      (tra_asm_inc routine)
    !                                                   (dyn_asm_inc routine)
    USE asmbkg
-   USE stpctl           ! time stepping control            (stp_ctl routine)
-   USE restart          ! ocean restart                    (rst_wri routine)
-   USE prtctl           ! Print control                    (prt_ctl routine)
+   USE stpctl          ! time stepping control            (stp_ctl routine)
+   USE restart         ! ocean restart                    (rst_wri routine)
+   USE prtctl          ! Print control                    (prt_ctl routine)
 
-   USE diaobs           ! Observation operator
+   USE diaobs          ! Observation operator
 
-   USE in_out_manager   ! I/O manager
-   USE iom              !
+   USE in_out_manager  ! I/O manager
+   USE iom             !
    USE lbclnk
-   USE timing           ! Timing
+   USE timing          ! Timing
 
 #if defined key_iomput
-   USE xios
+   USE xios            ! I/O server
 #endif
 #if defined key_agrif
