Index: NEMO/branches/2020/ticket_2459/cfgs/WED025/EXPREF/file_def_nemo-ice.xml
===================================================================
--- NEMO/branches/2020/ticket_2459/cfgs/WED025/EXPREF/file_def_nemo-ice.xml	(revision 12887)
+++ NEMO/branches/2020/ticket_2459/cfgs/WED025/EXPREF/file_def_nemo-ice.xml	(revision 12889)
@@ -78,29 +78,4 @@
      </file>
      
-     <file id="file22" name_suffix="_SBC_scalar" description="scalar variables" enabled=".true." >
-       <!-- global contents -->
-       <field field_ref="ibgvol_tot"     grid_ref="grid_1point"   name="ibgvol_tot"   />
-       <field field_ref="sbgvol_tot"     grid_ref="grid_1point"   name="sbgvol_tot"   />
-       <field field_ref="ibgarea_tot"    grid_ref="grid_1point"   name="ibgarea_tot"  />
-       <field field_ref="ibgsalt_tot"    grid_ref="grid_1point"   name="ibgsalt_tot"  />
-       <field field_ref="ibgheat_tot"    grid_ref="grid_1point"   name="ibgheat_tot"  />
-       <field field_ref="sbgheat_tot"    grid_ref="grid_1point"   name="sbgheat_tot"  />
-       
-       <!-- global drifts (conservation checks) -->
-       <field field_ref="ibgvolume"      grid_ref="grid_1point"   name="ibgvolume"    />
-       <field field_ref="ibgsaltco"      grid_ref="grid_1point"   name="ibgsaltco"    />
-       <field field_ref="ibgheatco"      grid_ref="grid_1point"   name="ibgheatco"    />
-       <field field_ref="ibgheatfx"      grid_ref="grid_1point"   name="ibgheatfx"    />
-       
-       <!-- global forcings  -->
-       <field field_ref="ibgfrcvoltop"   grid_ref="grid_1point"   name="ibgfrcvoltop" />
-       <field field_ref="ibgfrcvolbot"   grid_ref="grid_1point"   name="ibgfrcvolbot" />
-       <field field_ref="ibgfrctemtop"   grid_ref="grid_1point"   name="ibgfrctemtop" />
-       <field field_ref="ibgfrctembot"   grid_ref="grid_1point"   name="ibgfrctembot" />
-       <field field_ref="ibgfrcsal"      grid_ref="grid_1point"   name="ibgfrcsal"    />
-       <field field_ref="ibgfrchfxtop"   grid_ref="grid_1point"   name="ibgfrchfxtop" />
-       <field field_ref="ibgfrchfxbot"   grid_ref="grid_1point"   name="ibgfrchfxbot" />
-     </file>
-     
    </file_group>
    
Index: NEMO/branches/2020/ticket_2459/cfgs/WED025/EXPREF/namelist_cfg
===================================================================
--- NEMO/branches/2020/ticket_2459/cfgs/WED025/EXPREF/namelist_cfg	(revision 12887)
+++ NEMO/branches/2020/ticket_2459/cfgs/WED025/EXPREF/namelist_cfg	(revision 12889)
@@ -5,5 +5,5 @@
 !! namelists    2 - Surface boundary (namsbc, namsbc_flx, namsbc_blk, namsbc_cpl,
 !!                                    namsbc_sas, namtra_qsr, namsbc_rnf,
-!!                                    namsbc_isf, namsbc_iscpl, namsbc_apr, 
+!!                                    namisf, namsbc_apr, 
 !!                                    namsbc_ssr, namsbc_wave, namberg)
 !!              3 - lateral boundary (namlbc, namagrif, nambdy, nambdy_tide)
@@ -38,5 +38,5 @@
    nn_it000    =   1       !  first time step
    nn_itend    =  26280    !  last  time step (std 5475)
-   nn_date0    = 19760301  !  date at nit_0000 (format yyyymmdd) used if ln_rstart=F or (ln_rstart=T and nn_rstctl=0 or 1)
+   nn_date0    = 20000101  !  date at nit_0000 (format yyyymmdd) used if ln_rstart=F or (ln_rstart=T and nn_rstctl=0 or 1)
    ln_rstart   = .false.   !  start from rest (F) or from a restart file (T)
       nn_rstctl   =    2      !  restart control ==> activated only if ln_rstart=T
@@ -61,11 +61,11 @@
    ln_tsd_init = .true.          !  ocean initialisation
    ln_tsd_dmp  = .false.         !  T-S restoring   (see namtra_dmp)
-   
+
    cn_dir      = './'      !  root directory for the T-S data location
-   !___________!_________________________!___________________!___________!_____________!________!___________!__________________!__________!_______________!
-   !           !  file name              ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ ! weights filename ! rotation ! land/sea mask !
-   !           !                         !  (if <0  months)  !   name    !   (logical) !  (T/F) ! 'monthly' !                  ! pairing  !    filename   !
-   sn_tem = 'dta_temp_WED025'            ,         -12       , 'votemper',   .true.    , .true. , 'yearly'  ,    ''            ,    ''    ,    ''
-   sn_sal = 'dta_sal_WED025'             ,         -12       , 'vosaline',   .true.    , .true. , 'yearly'  ,    ''            ,    ''    ,    ''
+   !___________!_____________________!___________________!___________!_____________!________!___________!__________________!__________!_______________!
+   !           !  file name          ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ ! weights filename ! rotation ! land/sea mask !
+   !           !                     !  (if <0  months)  !   name    !   (logical) !  (T/F) ! 'monthly' !                  ! pairing  !    filename   !
+   sn_tem = 'WED025_init_JRA_200001.nc',       -12       , 'votemper',   .false.   , .true. , 'yearly'  ,    ''            ,    ''    ,    ''
+   sn_sal = 'WED025_init_JRA_200001.nc',       -12       , 'vosaline',   .false.   , .true. , 'yearly'  ,    ''            ,    ''    ,    ''
 /
 !-----------------------------------------------------------------------
@@ -124,5 +124,5 @@
                      ! Misc. options of sbc : 
    ln_traqsr   = .true.    !  Light penetration in the ocean            (T => fill namtra_qsr)
-   ln_dm2dc    = .true.    !  daily mean to diurnal cycle on short wave
+   ln_dm2dc    = .false.   !  daily mean to diurnal cycle on short wave
    ln_ssr      = .false.   !  Sea Surface Restoring on T and/or S       (T => fill namsbc_ssr)
    nn_fwb      = 0         !  FreshWater Budget: =0 unchecked
@@ -141,6 +141,6 @@
    ln_NCAR     = .true.   ! "NCAR"      algorithm   (Large and Yeager 2008)
    ln_COARE_3p0 = .false.   ! "COARE 3.0" algorithm   (Fairall et al. 2003)
-   ln_COARE_3p5 = .false.   ! "COARE 3.5" algorithm   (Edson et al. 2013)
-   ln_ECMWF    = .false.   ! "ECMWF"     algorithm   (IFS cycle 31)
+   ln_COARE_3p6 = .false.   ! "COARE 3.6" algorithm   (Edson et al. 2013)
+   ln_ECMWF     = .false.   ! "ECMWF"     algorithm   (IFS cycle 45r1)
 
    cn_dir      = './'      !  root directory for the bulk data location
@@ -148,13 +148,13 @@
    !           !  file name              ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ !       weights filename               ! rotation ! land/sea mask !
    !           !                         !  (if <0  months)  !   name    !   (logical) !  (T/F) ! 'monthly' !                                      ! pairing  !    filename   !
-   sn_wndi     = 'u10_core'              ,         6         , 'U_10_MOD',   .true.    , .false. , 'yearly'  , 'weights_bicubic_core.nc'           , 'Uwnd'   , ''
-   sn_wndj     = 'v10_core'              ,         6         , 'V_10_MOD',   .true.    , .false. , 'yearly'  , 'weights_bicubic_core.nc'           , 'Vwnd'   , ''
-   sn_qsr      = 'qsw_core'              ,        24         , 'SWDN_MOD',   .false.   , .false. , 'yearly'  , 'weights_bilin_core.nc'             , ''       , ''
-   sn_qlw      = 'qlw_core'              ,        24         , 'LWDN_MOD',   .false.   , .false. , 'yearly'  , 'weights_bilin_core.nc'             , ''       , ''
-   sn_tair     = 't10_core'              ,         6         , 'T_10_MOD',   .true.    , .false. , 'yearly'  , 'weights_bilin_core.nc'             , ''       , ''
-   sn_humi     = 'q10_core'              ,         6         , 'Q_10_MOD',   .true.    , .false. , 'yearly'  , 'weights_bilin_core.nc'             , ''       , ''
-   sn_prec     = 'precip_core'           ,        -1         , 'TPRECIP',    .true.    , .false. , 'yearly'  , 'weights_bilin_core.nc'             , ''       , ''
-   sn_snow     = 'snow_core'             ,        -1         , 'SNOW'    ,   .true.    , .false. , 'yearly'  , 'weights_bilin_core.nc'             , ''       , ''
-   sn_slp      = 'slp_core'              ,         6         , 'SLP'     ,   .true.    , .false. , 'yearly'  , 'weights_bilin_core.nc'             , ''       , ''
+   sn_wndi     = 'u10_JRA'              ,         3         , 'uas_10m' ,   .true.    , .false. , 'yearly'  , 'weights_bicubic_JRA.nc'           , 'Uwnd'   , ''
+   sn_wndj     = 'v10_JRA'              ,         3         , 'vas_10m' ,   .true.    , .false. , 'yearly'  , 'weights_bicubic_JRA.nc'           , 'Vwnd'   , ''
+   sn_qsr      = 'rsds_JRA'             ,         3         , 'rsds'    ,   .true.    , .false. , 'yearly'  , 'weights_bilin_JRA.nc'             , ''       , ''
+   sn_qlw      = 'rlds_JRA'             ,         3         , 'rlds'    ,   .true.    , .false. , 'yearly'  , 'weights_bilin_JRA.nc'             , ''       , ''
+   sn_tair     = 't10_JRA'              ,         3         , 'tas_10m' ,   .true.    , .false. , 'yearly'  , 'weights_bilin_JRA.nc'             , ''       , ''
+   sn_humi     = 'q10_JRA'              ,         3         , 'huss_10m',   .true.    , .false. , 'yearly'  , 'weights_bilin_JRA.nc'             , ''       , ''
+   sn_prec     = 'precip_JRA'           ,         3         , 'prto'    ,   .true.    , .false. , 'yearly'  , 'weights_bilin_JRA.nc'             , ''       , ''
+   sn_snow     = 'snow_JRA'             ,         3         , 'prsn'    ,   .true.    , .false. , 'yearly'  , 'weights_bilin_JRA.nc'             , ''       , ''
+   sn_slp      = 'slp_JRA'              ,         3         , 'psl'     ,   .true.    , .false. , 'yearly'  , 'weights_bilin_JRA.nc'             , ''       , ''
 /
 !-----------------------------------------------------------------------
@@ -201,5 +201,5 @@
    !           !  file name              ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ ! weights filename ! rotation ! land/sea mask !
    !           !                         !  (if <0  months)  !   name    !   (logical) !  (T/F) ! 'monthly' !                  ! pairing  !    filename   !
-   sn_rnf      = 'runoff_WED025'         ,  -1               , 'runoff'  ,   .true.    , .false., 'yearly'  , ''               , ''       , ''
+   sn_rnf      = 'WED025_icb'         ,  -1               , 'runoff'  ,   .true.    , .false., 'yearly'  , ''               , ''       , ''
 /
 !-----------------------------------------------------------------------
@@ -221,13 +221,14 @@
          cn_isfcav_mlt = '3eq'   ! ice shelf melting formulation (spe/2eq/3eq/oasis)
          !                       ! spe = fwfisf is read from a forcing field
-         !                       ! 2eq = ISOMIP  like: 2 equations formulation (Hunter et al., 2006)
-         !                       ! 3eq = ISOMIP+ like: 3 equations formulation (Asay-Davis et al., 2015)
+         !                       ! 2eq = ISOMIP  like: 2 equations formulation (Hunter et al., 2006 for a short description)
+         !                       ! 3eq = ISOMIP+ like: 3 equations formulation (Asay-Davis et al., 2016 for a short description)
          !                       ! oasis = fwfisf is given by oasis and pattern by file sn_isfcav_fwf
          !              !  cn_isfcav_mlt = 2eq or 3eq cases:
          cn_gammablk = 'vel'     ! scheme to compute gammat/s (spe,ad15,hj99)
-         !                       ! ad15 = velocity dependend Gamma (u* * gammat/s)  (Jenkins et al. 2010)
-         !                       ! hj99 = velocity and stability dependent Gamma    (Holland et al. 1999)
-         rn_gammat0  = 1.4e-2    ! gammat coefficient used in blk formula
-         rn_gammas0  = 4.e-4    ! gammas coefficient used in blk formula
+         !                       ! spe      = constant transfert velocity (rn_gammat0, rn_gammas0)
+         !                       ! vel      = velocity dependent transfert velocity (u* * gammat/s) (Asay-Davis et al. 2016 for a short description)
+         !                       ! vel_stab = velocity and stability dependent transfert coeficient (Holland et al. 1999 for a complete description)
+         rn_gammat0  = 1.4e-2    ! gammat coefficient used in spe, vel and vel_stab gamma computation method
+         rn_gammas0  = 4.0e-4    ! gammas coefficient used in spe, vel and vel_stab gamma computation method
          !
          rn_htbl     =  30.      ! thickness of the top boundary layer    (Losh et al. 2008)
@@ -255,7 +256,7 @@
          sn_isfpar_zmin = 'isfmlt_par',      -12.      , 'sozisfmin' ,  .false.    , .true.  , 'yearly'  ,    ''    ,   ''     ,    ''
          !* 'spe' and 'oasis' case
-         sn_isfpar_fwf = 'isfmlt_par' ,      -12.      , 'sofwfisf' ,  .false.    , .true.  , 'yearly'   ,    ''    ,   ''     ,    ''
+         sn_isfpar_fwf = 'isfmlt_par' ,      -12.      ,'sofwfisf' ,  .false.    , .true.  , 'yearly'   ,    ''    ,   ''     ,    ''
          !* 'bg03' case
-         sn_isfpar_Leff = 'isfmlt_par',       0.       , 'Leff'     ,  .false.    , .true.  , 'yearly'   ,    ''    ,   ''     ,    ''
+         sn_isfpar_Leff = 'isfmlt_par',       0.       ,'Leff'     ,  .false.    , .true.  , 'yearly'   ,    ''    ,   ''     ,    ''
       !
       ! ---------------- ice sheet coupling -------------------------------
@@ -300,8 +301,8 @@
    ln_tide     = .true.       ! Activate tides
       ln_tide_pot   = .false.               !  use tidal potential forcing
-      clname(1) = 'M2'  !  name of constituent - all tidal components must be set in namelist_cfg
-      clname(2) = 'S2'
-      clname(3) = 'K1'
-      clname(4) = 'O1'
+      sn_tide_cnames(1) = 'M2'  !  name of constituent - all tidal components must be set in namelist_cfg
+      sn_tide_cnames(2) = 'S2'
+      sn_tide_cnames(3) = 'K1'
+      sn_tide_cnames(4) = 'O1'
 /
 !-----------------------------------------------------------------------
@@ -340,20 +341,20 @@
    !           !  file name              ! frequency (hours) ! variable  ! time interp.!  clim  ! 'yearly'/ ! weights filename ! rotation ! land/sea mask !
    !           !                         !  (if <0  months)  !   name    !   (logical) !  (T/F) ! 'monthly' !                  ! pairing  !    filename   !
-   bn_ssh      =    'bdyT_ssh_WED025'    ,         -1        , 'sossheig' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     ''
-   bn_u2d      =    'bdyU_u2d_WED025'    ,         -1        , 'vobtcrtx' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     ''
-   bn_v2d      =    'bdyV_u2d_WED025'    ,         -1        , 'vobtcrty' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     ''
-   bn_u3d      =    'bdyU_u3d_WED025'    ,         -1        , 'vozocrtx' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     ''
-   bn_v3d      =    'bdyV_u3d_WED025'    ,         -1        , 'vomecrty' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     ''
-   bn_tem      =    'bdyT_tra_WED025'    ,         -1        , 'votemper' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     ''
-   bn_sal      =    'bdyT_tra_WED025'    ,         -1        , 'vosaline' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     ''
+   bn_ssh      =    'WED025_bdyT_ssh'    ,         -1        , 'sossheig' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     ''
+   bn_u2d      =    'WED025_bdyU_u2d'    ,         -1        , 'vobtcrtx' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     ''
+   bn_v2d      =    'WED025_bdyV_u2d'    ,         -1        , 'vobtcrty' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     ''
+   bn_u3d      =    'WED025_bdyU_u3d'    ,         -1        , 'vozocrtx' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     ''
+   bn_v3d      =    'WED025_bdyV_u3d'    ,         -1        , 'vomecrty' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     ''
+   bn_tem      =    'WED025_bdyT_tra'    ,         -1        , 'votemper' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     ''
+   bn_sal      =    'WED025_bdyT_tra'    ,         -1        , 'vosaline' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     ''
 !* for si3
-   bn_a_i      =    'bdyT_ice_WED025'    ,         -1        , 'ileadfra' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     ''
-   bn_h_i      =    'bdyT_ice_WED025'    ,         -1        , 'iicethic' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     ''
-   bn_h_s      =    'bdyT_ice_WED025'    ,         -1        , 'isnowthi' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     ''
+   bn_a_i      =    'WED025_bdyT_ice'    ,         -1        , 'ileadfra' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     ''
+   bn_h_i      =    'WED025_bdyT_ice'    ,         -1        , 'iicethic' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     ''
+   bn_h_s      =    'WED025_bdyT_ice'    ,         -1        , 'isnowthi' ,     .true. , .false., 'yearly'  ,    ''            ,   ''     ,     ''
 /
 !-----------------------------------------------------------------------
 &nambdy_tide   !  tidal forcing at open boundaries                      (default: OFF)
 !-----------------------------------------------------------------------
-   filtide          = 'bdytide_WED025_'         !  file name root of tidal forcing files
+   filtide          = 'WED025_bdytide_'         !  file name root of tidal forcing files
 /
 
@@ -658,12 +659,14 @@
 &namctl        !   Control prints                                       (default: OFF)
 !-----------------------------------------------------------------------
-   ln_ctl = .FALSE.                 ! Toggle all report printing on/off (T/F); Ignored if sn_cfctl%l_config is T
+   sn_cfctl%l_glochk = .FALSE.    ! Range sanity checks are local (F) or global (T). Set T for debugging only
+   sn_cfctl%l_allon  = .FALSE.    ! IF T activate all options. If F deactivate all unless l_config is T
      sn_cfctl%l_config = .TRUE.     ! IF .true. then control which reports are written with the following
-       sn_cfctl%l_runstat = .FALSE. ! switches and which areas produce reports with the proc integer settings.
+       sn_cfctl%l_runstat = .TRUE.  ! switches and which areas produce reports with the proc integer settings.
        sn_cfctl%l_trcstat = .FALSE. ! The default settings for the proc integers should ensure
        sn_cfctl%l_oceout  = .FALSE. ! that  all areas report.
        sn_cfctl%l_layout  = .FALSE. !
-       sn_cfctl%l_mppout  = .FALSE. !
-       sn_cfctl%l_mpptop  = .FALSE. !
+       sn_cfctl%l_prtctl  = .FALSE. !
+       sn_cfctl%l_prttrc  = .FALSE. !
+       sn_cfctl%l_oasout  = .FALSE. !
        sn_cfctl%procmin   = 0       ! Minimum area number for reporting [default:0]
        sn_cfctl%procmax   = 1000000 ! Maximum area number for reporting [default:1000000]
Index: NEMO/branches/2020/ticket_2459/cfgs/WED025/EXPREF/namelist_ice_cfg
===================================================================
--- NEMO/branches/2020/ticket_2459/cfgs/WED025/EXPREF/namelist_ice_cfg	(revision 12887)
+++ NEMO/branches/2020/ticket_2459/cfgs/WED025/EXPREF/namelist_ice_cfg	(revision 12889)
@@ -74,4 +74,18 @@
 &namini         !   Ice initialization
 !------------------------------------------------------------------------------
+   ln_iceini        = .true.          !  activate ice initialization (T) or not (F)
+   ln_iceini_file   = .true.          !  netcdf file provided for initialization (T) or not (F)
+   ! -- for ln_iceini_file = T
+   sn_hti = 'WED025_init_JRA_200001.nc', -12 ,'icethic_cea',  .false.  , .true., 'yearly'  , '' , '', ''
+   sn_hts = 'WED025_init_JRA_200001.nc', -12 ,'icesnow_cea',  .false.  , .true., 'yearly'  , '' , '', ''
+   sn_ati = 'WED025_init_JRA_200001.nc', -12 ,'ice_cover'  ,  .false.  , .true., 'yearly'  , '' , '', ''
+   sn_smi = 'NOT USED'              , -12 ,'smi'   ,  .false.  , .true., 'yearly'  , '' , '', ''
+   sn_tmi = 'NOT USED'              , -12 ,'tmi'   ,  .false.  , .true., 'yearly'  , '' , '', ''
+   sn_tsu = 'NOT USED'              , -12 ,'tsu'   ,  .false.  , .true., 'yearly'  , '' , '', ''
+   sn_tms = 'NOT USED'              , -12 ,'tms'   ,  .false.  , .true., 'yearly'  , '' , '', ''
+   !      melt ponds (be careful, sn_apd is the pond concentration (not fraction), so it differs from rn_apd)
+   sn_apd = 'NOT USED'              , -12 ,'apd'   ,  .false.  , .true., 'yearly'  , '' , '', ''
+   sn_hpd = 'NOT USED'              , -12 ,'hpd'   ,  .false.  , .true., 'yearly'  , '' , '', ''
+   cn_dir='./'
 /
 !------------------------------------------------------------------------------
Index: NEMO/branches/2020/ticket_2459/src/OCE/ISF/isfdiags.F90
===================================================================
--- NEMO/branches/2020/ticket_2459/src/OCE/ISF/isfdiags.F90	(revision 12887)
+++ NEMO/branches/2020/ticket_2459/src/OCE/ISF/isfdiags.F90	(revision 12889)
@@ -88,5 +88,5 @@
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: phtbl, pfrac  ! thickness of the tbl and fraction of last cell affected by the tbl
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pvar2d        ! 2d var to map in 3d
-      CHARACTER(LEN=256), INTENT(in) :: cdvar
+      CHARACTER(LEN=*), INTENT(in) :: cdvar
       !!---------------------------------------------------------------------
       INTEGER  :: ji, jj, jk                       ! loop indices
Index: NEMO/branches/2020/ticket_2459/tests/ISOMIP+/EXPREF/file_def_nemo-oce.xml
===================================================================
--- NEMO/branches/2020/ticket_2459/tests/ISOMIP+/EXPREF/file_def_nemo-oce.xml	(revision 12887)
+++ NEMO/branches/2020/ticket_2459/tests/ISOMIP+/EXPREF/file_def_nemo-oce.xml	(revision 12889)
@@ -21,25 +21,26 @@
       <file_group id="5d" output_freq="5d"  output_level="10" enabled=".TRUE.">  <!-- 5d files -->  
 
-      <file_group id="1m" output_freq="1mo" output_level="10" enabled=".TRUE."/> <!-- real monthly files -->
-	<file id="file1" output_freq="1mo" name_suffix="_grid_T" description="ocean T grid variables" >
-	  <field field_ref="toce"         name="votemper"  />
-	  <field field_ref="soce"         name="vosaline"  />
-	  <field field_ref="ssh"          name="sossheig"  />
+	<file id="file1" output_freq="5d" name_suffix="_grid_T" description="ocean T grid variables" >
+	  <field field_ref="toce"         name="votemper"  operation="average" freq_op="5d" > @toce_e3t / @e3t </field>
+	  <field field_ref="soce"         name="vosaline"  operation="average" freq_op="5d" > @soce_e3t / @e3t </field>
+	  <field field_ref="ssh"          name="sossheig" />
           <!-- variable for ice shelf -->
-          <field field_ref="fwfisf_cav"       name="sowflisf"  />
-          <field field_ref="isfgammat"    name="sogammat"  />
-          <field field_ref="isfgammas"    name="sogammas"  />
+          <field field_ref="fwfisf_cav"  name="sowflisf"  />
+          <field field_ref="isfgammat"   name="sogammat"  />
+          <field field_ref="isfgammas"   name="sogammas"  />
           <field field_ref="ttbl_cav"    name="ttbl"  />
-          <field field_ref="stbl"    name="stbl"  />
-          <field field_ref="utbl"    name="utbl"  />
-          <field field_ref="vtbl"    name="vtbl"  />
+          <field field_ref="stbl"        name="stbl"  />
+          <field field_ref="utbl"        name="utbl"  />
+          <field field_ref="vtbl"        name="vtbl"  />
         </file>
-	<file id="file2" output_freq="1mo" name_suffix="_grid_U" description="ocean U grid variables" >
-          <field field_ref="uoce"         name="vozocrtx"  />
+	<file id="file2" output_freq="5d" name_suffix="_grid_U" description="ocean U grid variables" >
+          <field field_ref="uoce"         name="vozocrtx" operation="average" freq_op="5d" > @uoce_e3u / @e3u </field> />
         </file>
-	<file id="file3" output_freq="1mo" name_suffix="_grid_V" description="ocean V grid variables" >
-          <field field_ref="voce"         name="vomecrty"  /> 
+	<file id="file3" output_freq="5d" name_suffix="_grid_V" description="ocean V grid variables" >
+          <field field_ref="voce"         name="vomecrty" operation="average" freq_op="5d" > @voce_e3v / @e3v </field> /> 
         </file>
       </file_group>
+
+      <file_group id="1m" output_freq="1mo" output_level="10" enabled=".TRUE."/> <!-- real monthly files -->
       <file_group id="2m" output_freq="2mo" output_level="10" enabled=".TRUE."/> <!-- real 2m files -->
       <file_group id="3m" output_freq="3mo" output_level="10" enabled=".TRUE."/> <!-- real 3m files -->
Index: NEMO/branches/2020/ticket_2459/tests/ISOMIP+/EXPREF/namelist_cfg
===================================================================
--- NEMO/branches/2020/ticket_2459/tests/ISOMIP+/EXPREF/namelist_cfg	(revision 12887)
+++ NEMO/branches/2020/ticket_2459/tests/ISOMIP+/EXPREF/namelist_cfg	(revision 12889)
@@ -308,7 +308,5 @@
 &nameos        !   ocean Equation Of Seawater                           (default: NO selection)
 !-----------------------------------------------------------------------
-   ln_teos10   = .false.         !  = Use TEOS-10
-   ln_eos80    = .false.         !  = Use EOS80
-   ln_leos     = .true.          !  = Use S-EOS (simplified Eq.)
+   ln_leos     = .true.          !  = Use L-EOS (linear Eq.)
                                  !
    !                     ! S-EOS coefficients (ln_seos=T):
Index: NEMO/branches/2020/ticket_2459/tests/ISOMIP+/MY_SRC/dtatsd.F90
===================================================================
--- NEMO/branches/2020/ticket_2459/tests/ISOMIP+/MY_SRC/dtatsd.F90	(revision 12887)
+++ NEMO/branches/2020/ticket_2459/tests/ISOMIP+/MY_SRC/dtatsd.F90	(revision 12889)
@@ -36,4 +36,6 @@
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_tsddmp ! structure of input SST (file informations, fields read)
 
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -67,8 +69,6 @@
       ierr0 = 0  ;  ierr1 = 0  ;  ierr2 = 0  ;  ierr3 = 0
       !
-      REWIND( numnam_ref )              ! Namelist namtsd in reference namelist : 
       READ  ( numnam_ref, namtsd, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtsd in reference namelist' )
-      REWIND( numnam_cfg )              ! Namelist namtsd in configuration namelist : Parameters of the run
       READ  ( numnam_cfg, namtsd, IOSTAT = ios, ERR = 902 )
 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtsd in configuration namelist' )
@@ -191,32 +191,30 @@
          ENDIF
          !
-         DO jj = 1, jpj                         ! vertical interpolation of T & S
-            DO ji = 1, jpi
-               DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points
-                  zl = gdept_0(ji,jj,jk)
-                  IF(     zl < gdept_1d(1  ) ) THEN          ! above the first level of data
-                     ztp(jk) =  ptsd(ji,jj,1    ,jp_tem)
-                     zsp(jk) =  ptsd(ji,jj,1    ,jp_sal)
-                  ELSEIF( zl > gdept_1d(jpk) ) THEN          ! below the last level of data
-                     ztp(jk) =  ptsd(ji,jj,jpkm1,jp_tem)
-                     zsp(jk) =  ptsd(ji,jj,jpkm1,jp_sal)
-                  ELSE                                      ! inbetween : vertical interpolation between jkk & jkk+1
-                     DO jkk = 1, jpkm1                                  ! when  gdept(jkk) < zl < gdept(jkk+1)
-                        IF( (zl-gdept_1d(jkk)) * (zl-gdept_1d(jkk+1)) <= 0._wp ) THEN
-                           zi = ( zl - gdept_1d(jkk) ) / (gdept_1d(jkk+1)-gdept_1d(jkk))
-                           ztp(jk) = ptsd(ji,jj,jkk,jp_tem) + ( ptsd(ji,jj,jkk+1,jp_tem) - ptsd(ji,jj,jkk,jp_tem) ) * zi 
-                           zsp(jk) = ptsd(ji,jj,jkk,jp_sal) + ( ptsd(ji,jj,jkk+1,jp_sal) - ptsd(ji,jj,jkk,jp_sal) ) * zi
-                        ENDIF
-                     END DO
-                  ENDIF
-               END DO
-               DO jk = 1, jpkm1
-                  ptsd(ji,jj,jk,jp_tem) = ztp(jk) * tmask(ji,jj,jk)     ! mask required for mixed zps-s-coord
-                  ptsd(ji,jj,jk,jp_sal) = zsp(jk) * tmask(ji,jj,jk)
-               END DO
-               ptsd(ji,jj,jpk,jp_tem) = 0._wp
-               ptsd(ji,jj,jpk,jp_sal) = 0._wp
+         DO_2D_11_11
+            DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points
+               zl = gdept_0(ji,jj,jk)
+               IF(     zl < gdept_1d(1  ) ) THEN          ! above the first level of data
+                  ztp(jk) =  ptsd(ji,jj,1    ,jp_tem)
+                  zsp(jk) =  ptsd(ji,jj,1    ,jp_sal)
+               ELSEIF( zl > gdept_1d(jpk) ) THEN          ! below the last level of data
+                  ztp(jk) =  ptsd(ji,jj,jpkm1,jp_tem)
+                  zsp(jk) =  ptsd(ji,jj,jpkm1,jp_sal)
+               ELSE                                      ! inbetween : vertical interpolation between jkk & jkk+1
+                  DO jkk = 1, jpkm1                                  ! when  gdept(jkk) < zl < gdept(jkk+1)
+                     IF( (zl-gdept_1d(jkk)) * (zl-gdept_1d(jkk+1)) <= 0._wp ) THEN
+                        zi = ( zl - gdept_1d(jkk) ) / (gdept_1d(jkk+1)-gdept_1d(jkk))
+                        ztp(jk) = ptsd(ji,jj,jkk,jp_tem) + ( ptsd(ji,jj,jkk+1,jp_tem) - ptsd(ji,jj,jkk,jp_tem) ) * zi 
+                        zsp(jk) = ptsd(ji,jj,jkk,jp_sal) + ( ptsd(ji,jj,jkk+1,jp_sal) - ptsd(ji,jj,jkk,jp_sal) ) * zi
+                     ENDIF
+                  END DO
+               ENDIF
             END DO
-         END DO
+            DO jk = 1, jpkm1
+               ptsd(ji,jj,jk,jp_tem) = ztp(jk) * tmask(ji,jj,jk)     ! mask required for mixed zps-s-coord
+               ptsd(ji,jj,jk,jp_sal) = zsp(jk) * tmask(ji,jj,jk)
+            END DO
+            ptsd(ji,jj,jpk,jp_tem) = 0._wp
+            ptsd(ji,jj,jpk,jp_sal) = 0._wp
+         END_2D
          ! 
       ELSE                                !==   z- or zps- coordinate   ==!
@@ -226,20 +224,18 @@
          !
          IF( ln_zps ) THEN                      ! zps-coordinate (partial steps) interpolation at the last ocean level
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ik = mbkt(ji,jj) 
-                  IF( ik > 1 ) THEN
-                     zl = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
-                     ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik-1,jp_tem)
-                     ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik-1,jp_sal)
-                  ENDIF
-                  ik = mikt(ji,jj)
-                  IF( ik > 1 ) THEN
-                     zl = ( gdept_0(ji,jj,ik) - gdept_1d(ik) ) / ( gdept_1d(ik+1) - gdept_1d(ik) ) 
-                     ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik+1,jp_tem)
-                     ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik+1,jp_sal)
-                  END IF
-               END DO
-            END DO
+            DO_2D_11_11
+               ik = mbkt(ji,jj) 
+               IF( ik > 1 ) THEN
+                  zl = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
+                  ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik-1,jp_tem)
+                  ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik-1,jp_sal)
+               ENDIF
+               ik = mikt(ji,jj)
+               IF( ik > 1 ) THEN
+                  zl = ( gdept_0(ji,jj,ik) - gdept_1d(ik) ) / ( gdept_1d(ik+1) - gdept_1d(ik) ) 
+                  ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik+1,jp_tem)
+                  ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik+1,jp_sal)
+               END IF
+            END_2D
          ENDIF
          !
Index: NEMO/branches/2020/ticket_2459/tests/ISOMIP+/MY_SRC/eosbn2.F90
===================================================================
--- NEMO/branches/2020/ticket_2459/tests/ISOMIP+/MY_SRC/eosbn2.F90	(revision 12887)
+++ NEMO/branches/2020/ticket_2459/tests/ISOMIP+/MY_SRC/eosbn2.F90	(revision 12889)
@@ -180,4 +180,6 @@
    REAL(wp) ::   BPE002
 
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -241,12 +243,129 @@
       CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
          !
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  !
-                  zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth
-                  zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature
-                  zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
-                  ztm = tmask(ji,jj,jk)                                         ! tmask
+         DO_3D_11_11( 1, jpkm1 )
+            !
+            zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth
+            zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature
+            zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
+            ztm = tmask(ji,jj,jk)                                         ! tmask
+            !
+            zn3 = EOS013*zt   &
+               &   + EOS103*zs+EOS003
+               !
+            zn2 = (EOS022*zt   &
+               &   + EOS112*zs+EOS012)*zt   &
+               &   + (EOS202*zs+EOS102)*zs+EOS002
+               !
+            zn1 = (((EOS041*zt   &
+               &   + EOS131*zs+EOS031)*zt   &
+               &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   &
+               &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   &
+               &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
+               !
+            zn0 = (((((EOS060*zt   &
+               &   + EOS150*zs+EOS050)*zt   &
+               &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   &
+               &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   &
+               &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   &
+               &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   &
+               &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
+               !
+            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+            !
+            prd(ji,jj,jk) = (  zn * r1_rho0 - 1._wp  ) * ztm  ! density anomaly (masked)
+            !
+         END_3D
+         !
+      CASE( np_seos )                !==  simplified EOS  ==!
+         !
+         DO_3D_11_11( 1, jpkm1 )
+            zt  = pts  (ji,jj,jk,jp_tem) - 10._wp
+            zs  = pts  (ji,jj,jk,jp_sal) - 35._wp
+            zh  = pdep (ji,jj,jk)
+            ztm = tmask(ji,jj,jk)
+            !
+            zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt   &
+               &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   &
+               &  - rn_nu * zt * zs
+               !                                 
+            prd(ji,jj,jk) = zn * r1_rho0 * ztm                ! density anomaly (masked)
+         END_3D
+         !
+      CASE( np_leos )                !==  linear ISOMIP EOS  ==!
+         !
+         DO_3D_11_11( 1, jpkm1 )
+            zt  = pts  (ji,jj,jk,jp_tem) - (-1._wp)
+            zs  = pts  (ji,jj,jk,jp_sal) - 34.2_wp
+            zh  = pdep (ji,jj,jk)
+            ztm = tmask(ji,jj,jk)
+            !
+            zn =  rho0 * ( - rn_a0 * zt + rn_b0 * zs )
+            !                                 
+            prd(ji,jj,jk) = zn * r1_rho0 * ztm                ! density anomaly (masked)
+         END_3D
+         !
+      END SELECT
+      !
+      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-insitu  : ', kdim=jpk )
+      !
+      IF( ln_timing )   CALL timing_stop('eos-insitu')
+      !
+   END SUBROUTINE eos_insitu
+
+
+   SUBROUTINE eos_insitu_pot( pts, prd, prhop, pdep )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE eos_insitu_pot  ***
+      !!
+      !! ** Purpose :   Compute the in situ density (ratio rho/rho0) and the
+      !!      potential volumic mass (Kg/m3) from potential temperature and
+      !!      salinity fields using an equation of state selected in the
+      !!     namelist.
+      !!
+      !! ** Action  : - prd  , the in situ density (no units)
+      !!              - prhop, the potential volumic mass (Kg/m3)
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celsius]
+      !                                                                ! 2 : salinity               [psu]
+      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prd    ! in situ density            [-]
+      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prhop  ! potential density (surface referenced)
+      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pdep   ! depth                      [m]
+      !
+      INTEGER  ::   ji, jj, jk, jsmp             ! dummy loop indices
+      INTEGER  ::   jdof
+      REAL(wp) ::   zt , zh , zstemp, zs , ztm   ! local scalars
+      REAL(wp) ::   zn , zn0, zn1, zn2, zn3      !   -      -
+      REAL(wp), DIMENSION(:), ALLOCATABLE :: zn0_sto, zn_sto, zsign    ! local vectors
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('eos-pot')
+      !
+      SELECT CASE ( neos )
+      !
+      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
+         !
+         ! Stochastic equation of state
+         IF ( ln_sto_eos ) THEN
+            ALLOCATE(zn0_sto(1:2*nn_sto_eos))
+            ALLOCATE(zn_sto(1:2*nn_sto_eos))
+            ALLOCATE(zsign(1:2*nn_sto_eos))
+            DO jsmp = 1, 2*nn_sto_eos, 2
+              zsign(jsmp)   = 1._wp
+              zsign(jsmp+1) = -1._wp
+            END DO
+            !
+            DO_3D_11_11( 1, jpkm1 )
+               !
+               ! compute density (2*nn_sto_eos) times:
+               ! (1) for t+dt, s+ds (with the random TS fluctutation computed in sto_pts)
+               ! (2) for t-dt, s-ds (with the opposite fluctuation)
+               DO jsmp = 1, nn_sto_eos*2
+                  jdof   = (jsmp + 1) / 2
+                  zh     = pdep(ji,jj,jk) * r1_Z0                                  ! depth
+                  zt     = (pts (ji,jj,jk,jp_tem) + pts_ran(ji,jj,jk,jp_tem,jdof) * zsign(jsmp)) * r1_T0    ! temperature
+                  zstemp = pts  (ji,jj,jk,jp_sal) + pts_ran(ji,jj,jk,jp_sal,jdof) * zsign(jsmp)
+                  zs     = SQRT( ABS( zstemp + rdeltaS ) * r1_S0 )   ! square root salinity
+                  ztm    = tmask(ji,jj,jk)                                         ! tmask
                   !
                   zn3 = EOS013*zt   &
@@ -263,5 +382,5 @@
                      &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
                      !
-                  zn0 = (((((EOS060*zt   &
+                  zn0_sto(jsmp) = (((((EOS060*zt   &
                      &   + EOS150*zs+EOS050)*zt   &
                      &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   &
@@ -271,278 +390,26 @@
                      &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
                      !
-                  zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+                  zn_sto(jsmp)  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0_sto(jsmp)
+               END DO
+               !
+               ! compute stochastic density as the mean of the (2*nn_sto_eos) densities
+               prhop(ji,jj,jk) = 0._wp ; prd(ji,jj,jk) = 0._wp
+               DO jsmp = 1, nn_sto_eos*2
+                  prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp)                      ! potential density referenced at the surface
                   !
-                  prd(ji,jj,jk) = (  zn * r1_rho0 - 1._wp  ) * ztm  ! density anomaly (masked)
-                  !
+                  prd(ji,jj,jk) = prd(ji,jj,jk) + (  zn_sto(jsmp) * r1_rho0 - 1._wp  )   ! density anomaly (masked)
                END DO
-            END DO
-         END DO
-         !
-      CASE( np_seos )                !==  simplified EOS  ==!
-         !
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zt  = pts  (ji,jj,jk,jp_tem) - 10._wp
-                  zs  = pts  (ji,jj,jk,jp_sal) - 35._wp
-                  zh  = pdep (ji,jj,jk)
-                  ztm = tmask(ji,jj,jk)
-                  !
-                  zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt   &
-                     &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   &
-                     &  - rn_nu * zt * zs
-                     !                                 
-                  prd(ji,jj,jk) = zn * r1_rho0 * ztm                ! density anomaly (masked)
-               END DO
-            END DO
-         END DO
-         !
-      CASE( np_leos )                !==  linear ISOMIP EOS  ==!
-         !
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zt  = pts  (ji,jj,jk,jp_tem) - (-1._wp)
-                  zs  = pts  (ji,jj,jk,jp_sal) - 34.2_wp
-                  zh  = pdep (ji,jj,jk)
-                  ztm = tmask(ji,jj,jk)
-                  !
-                  zn =  rho0 * ( - rn_a0 * zt + rn_b0 * zs )
-                  !                                 
-                  prd(ji,jj,jk) = zn * r1_rho0 * ztm                ! density anomaly (masked)
-               END DO
-            END DO
-         END DO
-         !
-      END SELECT
-      !
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-insitu  : ', kdim=jpk )
-      !
-      IF( ln_timing )   CALL timing_stop('eos-insitu')
-      !
-   END SUBROUTINE eos_insitu
-
-
-   SUBROUTINE eos_insitu_pot( pts, prd, prhop, pdep )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE eos_insitu_pot  ***
-      !!
-      !! ** Purpose :   Compute the in situ density (ratio rho/rho0) and the
-      !!      potential volumic mass (Kg/m3) from potential temperature and
-      !!      salinity fields using an equation of state selected in the
-      !!     namelist.
-      !!
-      !! ** Action  : - prd  , the in situ density (no units)
-      !!              - prhop, the potential volumic mass (Kg/m3)
-      !!
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celsius]
-      !                                                                ! 2 : salinity               [psu]
-      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prd    ! in situ density            [-]
-      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prhop  ! potential density (surface referenced)
-      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pdep   ! depth                      [m]
-      !
-      INTEGER  ::   ji, jj, jk, jsmp             ! dummy loop indices
-      INTEGER  ::   jdof
-      REAL(wp) ::   zt , zh , zstemp, zs , ztm   ! local scalars
-      REAL(wp) ::   zn , zn0, zn1, zn2, zn3      !   -      -
-      REAL(wp), DIMENSION(:), ALLOCATABLE :: zn0_sto, zn_sto, zsign    ! local vectors
-      !!----------------------------------------------------------------------
-      !
-      IF( ln_timing )   CALL timing_start('eos-pot')
-      !
-      SELECT CASE ( neos )
-      !
-      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
-         !
-         ! Stochastic equation of state
-         IF ( ln_sto_eos ) THEN
-            ALLOCATE(zn0_sto(1:2*nn_sto_eos))
-            ALLOCATE(zn_sto(1:2*nn_sto_eos))
-            ALLOCATE(zsign(1:2*nn_sto_eos))
-            DO jsmp = 1, 2*nn_sto_eos, 2
-              zsign(jsmp)   = 1._wp
-              zsign(jsmp+1) = -1._wp
-            END DO
-            !
-            DO jk = 1, jpkm1
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     !
-                     ! compute density (2*nn_sto_eos) times:
-                     ! (1) for t+dt, s+ds (with the random TS fluctutation computed in sto_pts)
-                     ! (2) for t-dt, s-ds (with the opposite fluctuation)
-                     DO jsmp = 1, nn_sto_eos*2
-                        jdof   = (jsmp + 1) / 2
-                        zh     = pdep(ji,jj,jk) * r1_Z0                                  ! depth
-                        zt     = (pts (ji,jj,jk,jp_tem) + pts_ran(ji,jj,jk,jp_tem,jdof) * zsign(jsmp)) * r1_T0    ! temperature
-                        zstemp = pts  (ji,jj,jk,jp_sal) + pts_ran(ji,jj,jk,jp_sal,jdof) * zsign(jsmp)
-                        zs     = SQRT( ABS( zstemp + rdeltaS ) * r1_S0 )   ! square root salinity
-                        ztm    = tmask(ji,jj,jk)                                         ! tmask
-                        !
-                        zn3 = EOS013*zt   &
-                           &   + EOS103*zs+EOS003
-                           !
-                        zn2 = (EOS022*zt   &
-                           &   + EOS112*zs+EOS012)*zt   &
-                           &   + (EOS202*zs+EOS102)*zs+EOS002
-                           !
-                        zn1 = (((EOS041*zt   &
-                           &   + EOS131*zs+EOS031)*zt   &
-                           &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   &
-                           &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   &
-                           &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
-                           !
-                        zn0_sto(jsmp) = (((((EOS060*zt   &
-                           &   + EOS150*zs+EOS050)*zt   &
-                           &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   &
-                           &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   &
-                           &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   &
-                           &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   &
-                           &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
-                           !
-                        zn_sto(jsmp)  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0_sto(jsmp)
-                     END DO
-                     !
-                     ! compute stochastic density as the mean of the (2*nn_sto_eos) densities
-                     prhop(ji,jj,jk) = 0._wp ; prd(ji,jj,jk) = 0._wp
-                     DO jsmp = 1, nn_sto_eos*2
-                        prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp)                      ! potential density referenced at the surface
-                        !
-                        prd(ji,jj,jk) = prd(ji,jj,jk) + (  zn_sto(jsmp) * r1_rho0 - 1._wp  )   ! density anomaly (masked)
-                     END DO
-                     prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos
-                     prd  (ji,jj,jk) = 0.5_wp * prd  (ji,jj,jk) * ztm / nn_sto_eos
-                  END DO
-               END DO
-            END DO
+               prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos
+               prd  (ji,jj,jk) = 0.5_wp * prd  (ji,jj,jk) * ztm / nn_sto_eos
+            END_3D
             DEALLOCATE(zn0_sto,zn_sto,zsign)
          ! Non-stochastic equation of state
          ELSE
-            DO jk = 1, jpkm1
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     !
-                     zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth
-                     zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature
-                     zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
-                     ztm = tmask(ji,jj,jk)                                         ! tmask
-                     !
-                     zn3 = EOS013*zt   &
-                        &   + EOS103*zs+EOS003
-                        !
-                     zn2 = (EOS022*zt   &
-                        &   + EOS112*zs+EOS012)*zt   &
-                        &   + (EOS202*zs+EOS102)*zs+EOS002
-                        !
-                     zn1 = (((EOS041*zt   &
-                        &   + EOS131*zs+EOS031)*zt   &
-                        &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   &
-                        &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   &
-                        &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
-                        !
-                     zn0 = (((((EOS060*zt   &
-                        &   + EOS150*zs+EOS050)*zt   &
-                        &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   &
-                        &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   &
-                        &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   &
-                        &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   &
-                        &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
-                        !
-                     zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
-                     !
-                     prhop(ji,jj,jk) = zn0 * ztm                           ! potential density referenced at the surface
-                     !
-                     prd(ji,jj,jk) = (  zn * r1_rho0 - 1._wp  ) * ztm      ! density anomaly (masked)
-                  END DO
-               END DO
-            END DO
-         ENDIF
-         
-      CASE( np_seos )                !==  simplified EOS  ==!
-         !
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zt  = pts  (ji,jj,jk,jp_tem) - 10._wp
-                  zs  = pts  (ji,jj,jk,jp_sal) - 35._wp
-                  zh  = pdep (ji,jj,jk)
-                  ztm = tmask(ji,jj,jk)
-                  !                                                     ! potential density referenced at the surface
-                  zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt ) * zt   &
-                     &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs ) * zs   &
-                     &  - rn_nu * zt * zs
-                  prhop(ji,jj,jk) = ( rho0 + zn ) * ztm
-                  !                                                     ! density anomaly (masked)
-                  zn = zn - ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zh
-                  prd(ji,jj,jk) = zn * r1_rho0 * ztm
-                  !
-               END DO
-            END DO
-         END DO
-         !
-      CASE( np_leos )                !==  linear ISOMIP EOS  ==!
-         !
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zt  = pts  (ji,jj,jk,jp_tem) - (-1._wp)
-                  zs  = pts  (ji,jj,jk,jp_sal) - 34.2_wp
-                  zh  = pdep (ji,jj,jk)
-                  ztm = tmask(ji,jj,jk)
-                  !                                                     ! potential density referenced at the surface
-                  zn =  rho0 * ( - rn_a0 * zt + rn_b0 * zs )
-                  prhop(ji,jj,jk) = ( rho0 + zn ) * ztm
-                  !                                                     ! density anomaly (masked)
-                  prd(ji,jj,jk) = zn * r1_rho0 * ztm
-                  !
-               END DO
-            END DO
-         END DO
-         !
-      END SELECT
-      !
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-pot: ', tab3d_2=prhop, clinfo2=' pot : ', kdim=jpk )
-      !
-      IF( ln_timing )   CALL timing_stop('eos-pot')
-      !
-   END SUBROUTINE eos_insitu_pot
-
-
-   SUBROUTINE eos_insitu_2d( pts, pdep, prd )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE eos_insitu_2d  ***
-      !!
-      !! ** Purpose :   Compute the in situ density (ratio rho/rho0) from
-      !!      potential temperature and salinity using an equation of state
-      !!      selected in the nameos namelist. * 2D field case
-      !!
-      !! ** Action  : - prd , the in situ density (no units) (unmasked)
-      !!
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celsius]
-      !                                                           ! 2 : salinity               [psu]
-      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(in   ) ::   pdep  ! depth                      [m]
-      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(  out) ::   prd   ! in situ density
-      !
-      INTEGER  ::   ji, jj, jk                ! dummy loop indices
-      REAL(wp) ::   zt , zh , zs              ! local scalars
-      REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      -
-      !!----------------------------------------------------------------------
-      !
-      IF( ln_timing )   CALL timing_start('eos2d')
-      !
-      prd(:,:) = 0._wp
-      !
-      SELECT CASE( neos )
-      !
-      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
-         !
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               !
-               zh  = pdep(ji,jj) * r1_Z0                                  ! depth
-               zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature
-               zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
+            DO_3D_11_11( 1, jpkm1 )
+               !
+               zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth
+               zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature
+               zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
+               ztm = tmask(ji,jj,jk)                                         ! tmask
                !
                zn3 = EOS013*zt   &
@@ -569,52 +436,148 @@
                zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
                !
-               prd(ji,jj) = zn * r1_rho0 - 1._wp               ! unmasked in situ density anomaly
-               !
-            END DO
-         END DO
-         !
-         CALL lbc_lnk( 'eosbn2', prd, 'T', 1. )                    ! Lateral boundary conditions
-         !
+               prhop(ji,jj,jk) = zn0 * ztm                           ! potential density referenced at the surface
+               !
+               prd(ji,jj,jk) = (  zn * r1_rho0 - 1._wp  ) * ztm      ! density anomaly (masked)
+            END_3D
+         ENDIF
+         
       CASE( np_seos )                !==  simplified EOS  ==!
          !
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               !
-               zt    = pts  (ji,jj,jp_tem)  - 10._wp
-               zs    = pts  (ji,jj,jp_sal)  - 35._wp
-               zh    = pdep (ji,jj)                         ! depth at the partial step level
-               !
-               zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt   &
-                  &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   &
-                  &  - rn_nu * zt * zs
-                  !
-               prd(ji,jj) = zn * r1_rho0               ! unmasked in situ density anomaly
-               !
-            END DO
-         END DO
-         !
-         CALL lbc_lnk( 'eosbn2', prd, 'T', 1. )                    ! Lateral boundary conditions
+         DO_3D_11_11( 1, jpkm1 )
+            zt  = pts  (ji,jj,jk,jp_tem) - 10._wp
+            zs  = pts  (ji,jj,jk,jp_sal) - 35._wp
+            zh  = pdep (ji,jj,jk)
+            ztm = tmask(ji,jj,jk)
+            !                                                     ! potential density referenced at the surface
+            zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt ) * zt   &
+               &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs ) * zs   &
+               &  - rn_nu * zt * zs
+            prhop(ji,jj,jk) = ( rho0 + zn ) * ztm
+            !                                                     ! density anomaly (masked)
+            zn = zn - ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zh
+            prd(ji,jj,jk) = zn * r1_rho0 * ztm
+            !
+         END_3D
+         !
+      CASE( np_leos )                !==  linear ISOMIP EOS  ==!
+         !
+         DO_3D_11_11( 1, jpkm1 )
+            zt  = pts  (ji,jj,jk,jp_tem) - (-1._wp)
+            zs  = pts  (ji,jj,jk,jp_sal) - 34.2_wp
+            zh  = pdep (ji,jj,jk)
+            ztm = tmask(ji,jj,jk)
+            !                                                     ! potential density referenced at the surface
+            zn =  rho0 * ( - rn_a0 * zt + rn_b0 * zs )
+            prhop(ji,jj,jk) = ( rho0 + zn ) * ztm
+            !                                                     ! density anomaly (masked)
+            prd(ji,jj,jk) = zn * r1_rho0 * ztm
+            !
+         END_3D
+         !
+      END SELECT
+      !
+      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-pot: ', tab3d_2=prhop, clinfo2=' pot : ', kdim=jpk )
+      !
+      IF( ln_timing )   CALL timing_stop('eos-pot')
+      !
+   END SUBROUTINE eos_insitu_pot
+
+
+   SUBROUTINE eos_insitu_2d( pts, pdep, prd )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE eos_insitu_2d  ***
+      !!
+      !! ** Purpose :   Compute the in situ density (ratio rho/rho0) from
+      !!      potential temperature and salinity using an equation of state
+      !!      selected in the nameos namelist. * 2D field case
+      !!
+      !! ** Action  : - prd , the in situ density (no units) (unmasked)
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celsius]
+      !                                                           ! 2 : salinity               [psu]
+      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(in   ) ::   pdep  ! depth                      [m]
+      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(  out) ::   prd   ! in situ density
+      !
+      INTEGER  ::   ji, jj, jk                ! dummy loop indices
+      REAL(wp) ::   zt , zh , zs              ! local scalars
+      REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      -
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('eos2d')
+      !
+      prd(:,:) = 0._wp
+      !
+      SELECT CASE( neos )
+      !
+      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
+         !
+         DO_2D_11_11
+            !
+            zh  = pdep(ji,jj) * r1_Z0                                  ! depth
+            zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature
+            zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
+            !
+            zn3 = EOS013*zt   &
+               &   + EOS103*zs+EOS003
+               !
+            zn2 = (EOS022*zt   &
+               &   + EOS112*zs+EOS012)*zt   &
+               &   + (EOS202*zs+EOS102)*zs+EOS002
+               !
+            zn1 = (((EOS041*zt   &
+               &   + EOS131*zs+EOS031)*zt   &
+               &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   &
+               &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   &
+               &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
+               !
+            zn0 = (((((EOS060*zt   &
+               &   + EOS150*zs+EOS050)*zt   &
+               &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   &
+               &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   &
+               &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   &
+               &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   &
+               &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
+               !
+            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+            !
+            prd(ji,jj) = zn * r1_rho0 - 1._wp               ! unmasked in situ density anomaly
+            !
+         END_2D
+         !
+      CASE( np_seos )                !==  simplified EOS  ==!
+         !
+         DO_2D_11_11
+            !
+            zt    = pts  (ji,jj,jp_tem)  - 10._wp
+            zs    = pts  (ji,jj,jp_sal)  - 35._wp
+            zh    = pdep (ji,jj)                         ! depth at the partial step level
+            !
+            zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt   &
+               &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   &
+               &  - rn_nu * zt * zs
+               !
+            prd(ji,jj) = zn * r1_rho0               ! unmasked in situ density anomaly
+            !
+         END_2D
          !
       CASE( np_leos )                !==  ISOMIP EOS  ==!
          !
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               !
-               zt    = pts  (ji,jj,jp_tem)  - (-1._wp)
-               zs    = pts  (ji,jj,jp_sal)  - 34.2_wp
-               zh    = pdep (ji,jj)                         ! depth at the partial step level
-               !
-               zn =  rho0 * ( - rn_a0 * zt + rn_b0 * zs )
-                  !
-               prd(ji,jj) = zn * r1_rho0               ! unmasked in situ density anomaly
-               !
-            END DO
-         END DO
-         !
-         CALL lbc_lnk( 'eosbn2', prd, 'T', 1. )                    ! Lateral boundary conditions
+         DO_2D_11_11
+            !
+            zt    = pts  (ji,jj,jp_tem)  - (-1._wp)
+            zs    = pts  (ji,jj,jp_sal)  - 34.2_wp
+            zh    = pdep (ji,jj)                         ! depth at the partial step level
+            !
+            zn =  rho0 * ( - rn_a0 * zt + rn_b0 * zs )
+            !
+            prd(ji,jj) = zn * r1_rho0               ! unmasked in situ density anomaly
+            !
+         END_2D
+         !
          !
       END SELECT
       !
-      IF(ln_ctl)   CALL prt_ctl( tab2d_1=prd, clinfo1=' eos2d: ' )
+      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=prd, clinfo1=' eos2d: ' )
       !
       IF( ln_timing )   CALL timing_stop('eos2d')
@@ -648,98 +611,86 @@
       CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
          !
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  !
-                  zh  = gdept(ji,jj,jk,Kmm) * r1_Z0                                ! depth
-                  zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature
-                  zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
-                  ztm = tmask(ji,jj,jk)                                         ! tmask
-                  !
-                  ! alpha
-                  zn3 = ALP003
-                  !
-                  zn2 = ALP012*zt + ALP102*zs+ALP002
-                  !
-                  zn1 = ((ALP031*zt   &
-                     &   + ALP121*zs+ALP021)*zt   &
-                     &   + (ALP211*zs+ALP111)*zs+ALP011)*zt   &
-                     &   + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001
-                     !
-                  zn0 = ((((ALP050*zt   &
-                     &   + ALP140*zs+ALP040)*zt   &
-                     &   + (ALP230*zs+ALP130)*zs+ALP030)*zt   &
-                     &   + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt   &
-                     &   + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt   &
-                     &   + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000
-                     !
-                  zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
-                  !
-                  pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm
-                  !
-                  ! beta
-                  zn3 = BET003
-                  !
-                  zn2 = BET012*zt + BET102*zs+BET002
-                  !
-                  zn1 = ((BET031*zt   &
-                     &   + BET121*zs+BET021)*zt   &
-                     &   + (BET211*zs+BET111)*zs+BET011)*zt   &
-                     &   + ((BET301*zs+BET201)*zs+BET101)*zs+BET001
-                     !
-                  zn0 = ((((BET050*zt   &
-                     &   + BET140*zs+BET040)*zt   &
-                     &   + (BET230*zs+BET130)*zs+BET030)*zt   &
-                     &   + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt   &
-                     &   + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt   &
-                     &   + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000
-                     !
-                  zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
-                  !
-                  pab(ji,jj,jk,jp_sal) = zn / zs * r1_rho0 * ztm
-                  !
-               END DO
-            END DO
-         END DO
+         DO_3D_11_11( 1, jpkm1 )
+            !
+            zh  = gdept(ji,jj,jk,Kmm) * r1_Z0                                ! depth
+            zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature
+            zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
+            ztm = tmask(ji,jj,jk)                                         ! tmask
+            !
+            ! alpha
+            zn3 = ALP003
+            !
+            zn2 = ALP012*zt + ALP102*zs+ALP002
+            !
+            zn1 = ((ALP031*zt   &
+               &   + ALP121*zs+ALP021)*zt   &
+               &   + (ALP211*zs+ALP111)*zs+ALP011)*zt   &
+               &   + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001
+               !
+            zn0 = ((((ALP050*zt   &
+               &   + ALP140*zs+ALP040)*zt   &
+               &   + (ALP230*zs+ALP130)*zs+ALP030)*zt   &
+               &   + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt   &
+               &   + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt   &
+               &   + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000
+               !
+            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+            !
+            pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm
+            !
+            ! beta
+            zn3 = BET003
+            !
+            zn2 = BET012*zt + BET102*zs+BET002
+            !
+            zn1 = ((BET031*zt   &
+               &   + BET121*zs+BET021)*zt   &
+               &   + (BET211*zs+BET111)*zs+BET011)*zt   &
+               &   + ((BET301*zs+BET201)*zs+BET101)*zs+BET001
+               !
+            zn0 = ((((BET050*zt   &
+               &   + BET140*zs+BET040)*zt   &
+               &   + (BET230*zs+BET130)*zs+BET030)*zt   &
+               &   + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt   &
+               &   + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt   &
+               &   + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000
+               !
+            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+            !
+            pab(ji,jj,jk,jp_sal) = zn / zs * r1_rho0 * ztm
+            !
+         END_3D
          !
       CASE( np_seos )                  !==  simplified EOS  ==!
          !
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zt  = pts (ji,jj,jk,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0)
-                  zs  = pts (ji,jj,jk,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0)
-                  zh  = gdept(ji,jj,jk,Kmm)                ! depth in meters at t-point
-                  ztm = tmask(ji,jj,jk)                  ! land/sea bottom mask = surf. mask
-                  !
-                  zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs
-                  pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm   ! alpha
-                  !
-                  zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt
-                  pab(ji,jj,jk,jp_sal) = zn * r1_rho0 * ztm   ! beta
-                  !
-               END DO
-            END DO
-         END DO
+         DO_3D_11_11( 1, jpkm1 )
+            zt  = pts (ji,jj,jk,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0)
+            zs  = pts (ji,jj,jk,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0)
+            zh  = gdept(ji,jj,jk,Kmm)                ! depth in meters at t-point
+            ztm = tmask(ji,jj,jk)                  ! land/sea bottom mask = surf. mask
+            !
+            zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs
+            pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm   ! alpha
+            !
+            zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt
+            pab(ji,jj,jk,jp_sal) = zn * r1_rho0 * ztm   ! beta
+            !
+         END_3D
          !
       CASE( np_leos )                  !==  linear ISOMIP EOS  ==!
          !
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zt  = pts (ji,jj,jk,jp_tem) - (-1._wp)
-                  zs  = pts (ji,jj,jk,jp_sal) - 34.2_wp   ! abs. salinity anomaly (s-S0)
-                  zh  = gdept(ji,jj,jk,Kmm)                 ! depth in meters at t-point
-                  ztm = tmask(ji,jj,jk)                   ! land/sea bottom mask = surf. mask
-                  !
-                  zn  = rn_a0 * rho0
-                  pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm   ! alpha
-                  !
-                  zn  = rn_b0 * rho0
-                  pab(ji,jj,jk,jp_sal) = zn * r1_rho0 * ztm   ! beta
-                  !
-               END DO
-            END DO
-         END DO
+         DO_3D_11_11( 1, jpkm1 )
+            zt  = pts (ji,jj,jk,jp_tem) - (-1._wp)
+            zs  = pts (ji,jj,jk,jp_sal) - 34.2_wp   ! abs. salinity anomaly (s-S0)
+            zh  = gdept(ji,jj,jk,Kmm)                 ! depth in meters at t-point
+            ztm = tmask(ji,jj,jk)                   ! land/sea bottom mask = surf. mask
+            !
+            zn  = rn_a0 * rho0
+            pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm   ! alpha
+            !
+            zn  = rn_b0 * rho0
+            pab(ji,jj,jk,jp_sal) = zn * r1_rho0 * ztm   ! beta
+            !
+         END_3D
          !
       CASE DEFAULT
@@ -749,6 +700,6 @@
       END SELECT
       !
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pab(:,:,:,jp_tem), clinfo1=' rab_3d_t: ', &
-         &                       tab3d_2=pab(:,:,:,jp_sal), clinfo2=' rab_3d_s : ', kdim=jpk )
+      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pab(:,:,:,jp_tem), clinfo1=' rab_3d_t: ', &
+         &                                  tab3d_2=pab(:,:,:,jp_sal), clinfo2=' rab_3d_s : ', kdim=jpk )
       !
       IF( ln_timing )   CALL timing_stop('rab_3d')
@@ -783,98 +734,86 @@
       CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
          !
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               !
-               zh  = pdep(ji,jj) * r1_Z0                                  ! depth
-               zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature
-               zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
-               !
-               ! alpha
-               zn3 = ALP003
-               !
-               zn2 = ALP012*zt + ALP102*zs+ALP002
-               !
-               zn1 = ((ALP031*zt   &
-                  &   + ALP121*zs+ALP021)*zt   &
-                  &   + (ALP211*zs+ALP111)*zs+ALP011)*zt   &
-                  &   + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001
-                  !
-               zn0 = ((((ALP050*zt   &
-                  &   + ALP140*zs+ALP040)*zt   &
-                  &   + (ALP230*zs+ALP130)*zs+ALP030)*zt   &
-                  &   + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt   &
-                  &   + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt   &
-                  &   + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000
-                  !
-               zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
-               !
-               pab(ji,jj,jp_tem) = zn * r1_rho0
-               !
-               ! beta
-               zn3 = BET003
-               !
-               zn2 = BET012*zt + BET102*zs+BET002
-               !
-               zn1 = ((BET031*zt   &
-                  &   + BET121*zs+BET021)*zt   &
-                  &   + (BET211*zs+BET111)*zs+BET011)*zt   &
-                  &   + ((BET301*zs+BET201)*zs+BET101)*zs+BET001
-                  !
-               zn0 = ((((BET050*zt   &
-                  &   + BET140*zs+BET040)*zt   &
-                  &   + (BET230*zs+BET130)*zs+BET030)*zt   &
-                  &   + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt   &
-                  &   + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt   &
-                  &   + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000
-                  !
-               zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
-               !
-               pab(ji,jj,jp_sal) = zn / zs * r1_rho0
-               !
-               !
-            END DO
-         END DO
-         !                            ! Lateral boundary conditions
-         CALL lbc_lnk_multi( 'eosbn2', pab(:,:,jp_tem), 'T', 1. , pab(:,:,jp_sal), 'T', 1. )                    
+         DO_2D_11_11
+            !
+            zh  = pdep(ji,jj) * r1_Z0                                  ! depth
+            zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature
+            zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
+            !
+            ! alpha
+            zn3 = ALP003
+            !
+            zn2 = ALP012*zt + ALP102*zs+ALP002
+            !
+            zn1 = ((ALP031*zt   &
+               &   + ALP121*zs+ALP021)*zt   &
+               &   + (ALP211*zs+ALP111)*zs+ALP011)*zt   &
+               &   + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001
+               !
+            zn0 = ((((ALP050*zt   &
+               &   + ALP140*zs+ALP040)*zt   &
+               &   + (ALP230*zs+ALP130)*zs+ALP030)*zt   &
+               &   + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt   &
+               &   + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt   &
+               &   + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000
+               !
+            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+            !
+            pab(ji,jj,jp_tem) = zn * r1_rho0
+            !
+            ! beta
+            zn3 = BET003
+            !
+            zn2 = BET012*zt + BET102*zs+BET002
+            !
+            zn1 = ((BET031*zt   &
+               &   + BET121*zs+BET021)*zt   &
+               &   + (BET211*zs+BET111)*zs+BET011)*zt   &
+               &   + ((BET301*zs+BET201)*zs+BET101)*zs+BET001
+               !
+            zn0 = ((((BET050*zt   &
+               &   + BET140*zs+BET040)*zt   &
+               &   + (BET230*zs+BET130)*zs+BET030)*zt   &
+               &   + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt   &
+               &   + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt   &
+               &   + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000
+               !
+            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+            !
+            pab(ji,jj,jp_sal) = zn / zs * r1_rho0
+            !
+            !
+         END_2D
          !
       CASE( np_seos )                  !==  simplified EOS  ==!
          !
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               !
-               zt    = pts  (ji,jj,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0)
-               zs    = pts  (ji,jj,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0)
-               zh    = pdep (ji,jj)                   ! depth at the partial step level
-               !
-               zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs
-               pab(ji,jj,jp_tem) = zn * r1_rho0   ! alpha
-               !
-               zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt
-               pab(ji,jj,jp_sal) = zn * r1_rho0   ! beta
-               !
-            END DO
-         END DO
-         !                            ! Lateral boundary conditions
-         CALL lbc_lnk_multi( 'eosbn2', pab(:,:,jp_tem), 'T', 1. , pab(:,:,jp_sal), 'T', 1. )                    
+         DO_2D_11_11
+            !
+            zt    = pts  (ji,jj,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0)
+            zs    = pts  (ji,jj,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0)
+            zh    = pdep (ji,jj)                   ! depth at the partial step level
+            !
+            zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs
+            pab(ji,jj,jp_tem) = zn * r1_rho0   ! alpha
+            !
+            zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt
+            pab(ji,jj,jp_sal) = zn * r1_rho0   ! beta
+            !
+         END_2D
          !
       CASE( np_leos )                  !==  linear ISOMIP EOS  ==!
          !
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               !
-               zt    = pts  (ji,jj,jp_tem) - (-1._wp)   ! pot. temperature anomaly (t-T0)
-               zs    = pts  (ji,jj,jp_sal) - 34.2_wp   ! abs. salinity anomaly (s-S0)
-               zh    = pdep (ji,jj)                   ! depth at the partial step level
-               !
-               zn  = rn_a0 * rho0
-               pab(ji,jj,jp_tem) = zn * r1_rho0   ! alpha
-               !
-               zn  = rn_b0 * rho0
-               pab(ji,jj,jp_sal) = zn * r1_rho0   ! beta
-               !
-            END DO
-         END DO
-         !
-         CALL lbc_lnk_multi( 'eosbn2', pab(:,:,jp_tem), 'T', 1. , pab(:,:,jp_sal), 'T', 1. )                    ! Lateral boundary conditions
+         DO_2D_11_11
+            !
+            zt    = pts  (ji,jj,jp_tem) - (-1._wp)   ! pot. temperature anomaly (t-T0)
+            zs    = pts  (ji,jj,jp_sal) - 34.2_wp   ! abs. salinity anomaly (s-S0)
+            zh    = pdep (ji,jj)                   ! depth at the partial step level
+            !
+            zn  = rn_a0 * rho0
+            pab(ji,jj,jp_tem) = zn * r1_rho0   ! alpha
+            !
+            zn  = rn_b0 * rho0
+            pab(ji,jj,jp_sal) = zn * r1_rho0   ! beta
+            !
+         END_2D
          !
       CASE DEFAULT
@@ -884,6 +823,6 @@
       END SELECT
       !
-      IF(ln_ctl)   CALL prt_ctl( tab2d_1=pab(:,:,jp_tem), clinfo1=' rab_2d_t: ', &
-         &                       tab2d_2=pab(:,:,jp_sal), clinfo2=' rab_2d_s : ' )
+      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=pab(:,:,jp_tem), clinfo1=' rab_2d_t: ', &
+         &                                  tab2d_2=pab(:,:,jp_sal), clinfo2=' rab_2d_s : ' )
       !
       IF( ln_timing )   CALL timing_stop('rab_2d')
@@ -1026,21 +965,17 @@
       IF( ln_timing )   CALL timing_start('bn2')
       !
-      DO jk = 2, jpkm1           ! interior points only (2=< jk =< jpkm1 )
-         DO jj = 1, jpj          ! surface and bottom value set to zero one for all in istate.F90
-            DO ji = 1, jpi
-               zrw =   ( gdepw(ji,jj,jk  ,Kmm) - gdept(ji,jj,jk,Kmm) )   &
-                  &  / ( gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm) ) 
-                  !
-               zaw = pab(ji,jj,jk,jp_tem) * (1. - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw 
-               zbw = pab(ji,jj,jk,jp_sal) * (1. - zrw) + pab(ji,jj,jk-1,jp_sal) * zrw
-               !
-               pn2(ji,jj,jk) = grav * (  zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )     &
-                  &                    - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )  )  &
-                  &            / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-      !
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2  : ', kdim=jpk )
+      DO_3D_11_11( 2, jpkm1 )
+         zrw =   ( gdepw(ji,jj,jk  ,Kmm) - gdept(ji,jj,jk,Kmm) )   &
+            &  / ( gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm) ) 
+            !
+         zaw = pab(ji,jj,jk,jp_tem) * (1. - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw 
+         zbw = pab(ji,jj,jk,jp_sal) * (1. - zrw) + pab(ji,jj,jk-1,jp_sal) * zrw
+         !
+         pn2(ji,jj,jk) = grav * (  zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )     &
+            &                    - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )  )  &
+            &            / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)
+      END_3D
+      !
+      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2  : ', kdim=jpk )
       !
       IF( ln_timing )   CALL timing_stop('bn2')
@@ -1078,28 +1013,26 @@
       z1_T0   = 1._wp/40._wp
       !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            !
-            zt  = ctmp   (ji,jj) * z1_T0
-            zs  = SQRT( ABS( psal(ji,jj) + zdeltaS ) * r1_S0 )
-            ztm = tmask(ji,jj,1)
-            !
-            zn = ((((-2.1385727895e-01_wp*zt   &
-               &   - 2.7674419971e-01_wp*zs+1.0728094330_wp)*zt   &
-               &   + (2.6366564313_wp*zs+3.3546960647_wp)*zs-7.8012209473_wp)*zt   &
-               &   + ((1.8835586562_wp*zs+7.3949191679_wp)*zs-3.3937395875_wp)*zs-5.6414948432_wp)*zt   &
-               &   + (((3.5737370589_wp*zs-1.5512427389e+01_wp)*zs+2.4625741105e+01_wp)*zs   &
-               &      +1.9912291000e+01_wp)*zs-3.2191146312e+01_wp)*zt   &
-               &   + ((((5.7153204649e-01_wp*zs-3.0943149543_wp)*zs+9.3052495181_wp)*zs   &
-               &      -9.4528934807_wp)*zs+3.1066408996_wp)*zs-4.3504021262e-01_wp
-               !
-            zd = (2.0035003456_wp*zt   &
-               &   -3.4570358592e-01_wp*zs+5.6471810638_wp)*zt   &
-               &   + (1.5393993508_wp*zs-6.9394762624_wp)*zs+1.2750522650e+01_wp
-               !
-            ptmp(ji,jj) = ( zt / z1_T0 + zn / zd ) * ztm
-               !
-         END DO
-      END DO
+      DO_2D_11_11
+         !
+         zt  = ctmp   (ji,jj) * z1_T0
+         zs  = SQRT( ABS( psal(ji,jj) + zdeltaS ) * r1_S0 )
+         ztm = tmask(ji,jj,1)
+         !
+         zn = ((((-2.1385727895e-01_wp*zt   &
+            &   - 2.7674419971e-01_wp*zs+1.0728094330_wp)*zt   &
+            &   + (2.6366564313_wp*zs+3.3546960647_wp)*zs-7.8012209473_wp)*zt   &
+            &   + ((1.8835586562_wp*zs+7.3949191679_wp)*zs-3.3937395875_wp)*zs-5.6414948432_wp)*zt   &
+            &   + (((3.5737370589_wp*zs-1.5512427389e+01_wp)*zs+2.4625741105e+01_wp)*zs   &
+            &      +1.9912291000e+01_wp)*zs-3.2191146312e+01_wp)*zt   &
+            &   + ((((5.7153204649e-01_wp*zs-3.0943149543_wp)*zs+9.3052495181_wp)*zs   &
+            &      -9.4528934807_wp)*zs+3.1066408996_wp)*zs-4.3504021262e-01_wp
+            !
+         zd = (2.0035003456_wp*zt   &
+            &   -3.4570358592e-01_wp*zs+5.6471810638_wp)*zt   &
+            &   + (1.5393993508_wp*zs-6.9394762624_wp)*zs+1.2750522650e+01_wp
+            !
+         ptmp(ji,jj) = ( zt / z1_T0 + zn / zd ) * ztm
+            !
+      END_2D
       !
       IF( ln_timing )   CALL timing_stop('eos_pt_from_ct')
@@ -1133,16 +1066,14 @@
          !
          z1_S0 = 1._wp / 35.16504_wp
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zs= SQRT( ABS( psal(ji,jj) ) * z1_S0 )           ! square root salinity
-               ptf(ji,jj) = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs &
-                  &          - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp
-            END DO
-         END DO
+         DO_2D_11_11
+            zs= SQRT( ABS( psal(ji,jj) ) * z1_S0 )           ! square root salinity
+            ptf(ji,jj) = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs &
+               &          - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp
+         END_2D
          ptf(:,:) = ptf(:,:) * psal(:,:)
          !
          IF( PRESENT( pdep ) )   ptf(:,:) = ptf(:,:) - 7.53e-4 * pdep(:,:)
          !
-      CASE ( np_eos80, np_leos )                !==  PT,SP (UNESCO formulation)  ==!
+      CASE ( np_eos80 )                !==  PT,SP (UNESCO formulation)  ==!
          !
          ptf(:,:) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(:,:) )   &
@@ -1190,5 +1121,5 @@
          IF( PRESENT( pdep ) )   ptf = ptf - 7.53e-4 * pdep
          !
-      CASE ( np_eos80, np_leos )                !==  PT,SP (UNESCO formulation)  ==!
+      CASE ( np_eos80 )                !==  PT,SP (UNESCO formulation)  ==!
          !
          ptf = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal )   &
@@ -1242,104 +1173,92 @@
       CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
          !
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  !
-                  zh  = gdept(ji,jj,jk,Kmm) * r1_Z0                                ! depth
-                  zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature
-                  zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
-                  ztm = tmask(ji,jj,jk)                                         ! tmask
-                  !
-                  ! potential energy non-linear anomaly
-                  zn2 = (PEN012)*zt   &
-                     &   + PEN102*zs+PEN002
-                     !
-                  zn1 = ((PEN021)*zt   &
-                     &   + PEN111*zs+PEN011)*zt   &
-                     &   + (PEN201*zs+PEN101)*zs+PEN001
-                     !
-                  zn0 = ((((PEN040)*zt   &
-                     &   + PEN130*zs+PEN030)*zt   &
-                     &   + (PEN220*zs+PEN120)*zs+PEN020)*zt   &
-                     &   + ((PEN310*zs+PEN210)*zs+PEN110)*zs+PEN010)*zt   &
-                     &   + (((PEN400*zs+PEN300)*zs+PEN200)*zs+PEN100)*zs+PEN000
-                     !
-                  zn  = ( zn2 * zh + zn1 ) * zh + zn0
-                  !
-                  ppen(ji,jj,jk)  = zn * zh * r1_rho0 * ztm
-                  !
-                  ! alphaPE non-linear anomaly
-                  zn2 = APE002
-                  !
-                  zn1 = (APE011)*zt   &
-                     &   + APE101*zs+APE001
-                     !
-                  zn0 = (((APE030)*zt   &
-                     &   + APE120*zs+APE020)*zt   &
-                     &   + (APE210*zs+APE110)*zs+APE010)*zt   &
-                     &   + ((APE300*zs+APE200)*zs+APE100)*zs+APE000
-                     !
-                  zn  = ( zn2 * zh + zn1 ) * zh + zn0
-                  !                              
-                  pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_rho0 * ztm
-                  !
-                  ! betaPE non-linear anomaly
-                  zn2 = BPE002
-                  !
-                  zn1 = (BPE011)*zt   &
-                     &   + BPE101*zs+BPE001
-                     !
-                  zn0 = (((BPE030)*zt   &
-                     &   + BPE120*zs+BPE020)*zt   &
-                     &   + (BPE210*zs+BPE110)*zs+BPE010)*zt   &
-                     &   + ((BPE300*zs+BPE200)*zs+BPE100)*zs+BPE000
-                     !
-                  zn  = ( zn2 * zh + zn1 ) * zh + zn0
-                  !                              
-                  pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_rho0 * ztm
-                  !
-               END DO
-            END DO
-         END DO
+         DO_3D_11_11( 1, jpkm1 )
+            !
+            zh  = gdept(ji,jj,jk,Kmm) * r1_Z0                                ! depth
+            zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature
+            zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
+            ztm = tmask(ji,jj,jk)                                         ! tmask
+            !
+            ! potential energy non-linear anomaly
+            zn2 = (PEN012)*zt   &
+               &   + PEN102*zs+PEN002
+               !
+            zn1 = ((PEN021)*zt   &
+               &   + PEN111*zs+PEN011)*zt   &
+               &   + (PEN201*zs+PEN101)*zs+PEN001
+               !
+            zn0 = ((((PEN040)*zt   &
+               &   + PEN130*zs+PEN030)*zt   &
+               &   + (PEN220*zs+PEN120)*zs+PEN020)*zt   &
+               &   + ((PEN310*zs+PEN210)*zs+PEN110)*zs+PEN010)*zt   &
+               &   + (((PEN400*zs+PEN300)*zs+PEN200)*zs+PEN100)*zs+PEN000
+               !
+            zn  = ( zn2 * zh + zn1 ) * zh + zn0
+            !
+            ppen(ji,jj,jk)  = zn * zh * r1_rho0 * ztm
+            !
+            ! alphaPE non-linear anomaly
+            zn2 = APE002
+            !
+            zn1 = (APE011)*zt   &
+               &   + APE101*zs+APE001
+               !
+            zn0 = (((APE030)*zt   &
+               &   + APE120*zs+APE020)*zt   &
+               &   + (APE210*zs+APE110)*zs+APE010)*zt   &
+               &   + ((APE300*zs+APE200)*zs+APE100)*zs+APE000
+               !
+            zn  = ( zn2 * zh + zn1 ) * zh + zn0
+            !                              
+            pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_rho0 * ztm
+            !
+            ! betaPE non-linear anomaly
+            zn2 = BPE002
+            !
+            zn1 = (BPE011)*zt   &
+               &   + BPE101*zs+BPE001
+               !
+            zn0 = (((BPE030)*zt   &
+               &   + BPE120*zs+BPE020)*zt   &
+               &   + (BPE210*zs+BPE110)*zs+BPE010)*zt   &
+               &   + ((BPE300*zs+BPE200)*zs+BPE100)*zs+BPE000
+               !
+            zn  = ( zn2 * zh + zn1 ) * zh + zn0
+            !                              
+            pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_rho0 * ztm
+            !
+         END_3D
          !
       CASE( np_seos )                !==  Vallis (2006) simplified EOS  ==!
          !
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zt  = pts(ji,jj,jk,jp_tem) - 10._wp  ! temperature anomaly (t-T0)
-                  zs = pts (ji,jj,jk,jp_sal) - 35._wp  ! abs. salinity anomaly (s-S0)
-                  zh  = gdept(ji,jj,jk,Kmm)              ! depth in meters  at t-point
-                  ztm = tmask(ji,jj,jk)                ! tmask
-                  zn  = 0.5_wp * zh * r1_rho0 * ztm
-                  !                                    ! Potential Energy
-                  ppen(ji,jj,jk) = ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zn
-                  !                                    ! alphaPE
-                  pab_pe(ji,jj,jk,jp_tem) = - rn_a0 * rn_mu1 * zn
-                  pab_pe(ji,jj,jk,jp_sal) =   rn_b0 * rn_mu2 * zn
-                  !
-               END DO
-            END DO
-         END DO
+         DO_3D_11_11( 1, jpkm1 )
+            zt  = pts(ji,jj,jk,jp_tem) - 10._wp  ! temperature anomaly (t-T0)
+            zs = pts (ji,jj,jk,jp_sal) - 35._wp  ! abs. salinity anomaly (s-S0)
+            zh  = gdept(ji,jj,jk,Kmm)              ! depth in meters  at t-point
+            ztm = tmask(ji,jj,jk)                ! tmask
+            zn  = 0.5_wp * zh * r1_rho0 * ztm
+            !                                    ! Potential Energy
+            ppen(ji,jj,jk) = ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zn
+            !                                    ! alphaPE
+            pab_pe(ji,jj,jk,jp_tem) = - rn_a0 * rn_mu1 * zn
+            pab_pe(ji,jj,jk,jp_sal) =   rn_b0 * rn_mu2 * zn
+            !
+         END_3D
          !
       CASE( np_leos )                !==  linear ISOMIP EOS  ==!
          !
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zt  = pts(ji,jj,jk,jp_tem) - (-1._wp)  ! temperature anomaly (t-T0)
-                  zs = pts (ji,jj,jk,jp_sal) - 34.2_wp   ! abs. salinity anomaly (s-S0)
-                  zh  = gdept(ji,jj,jk,Kmm)                ! depth in meters  at t-point
-                  ztm = tmask(ji,jj,jk)                  ! tmask
-                  zn  = 0.5_wp * zh * r1_rho0 * ztm
-                  !                                    ! Potential Energy
-                  ppen(ji,jj,jk) = 0.
-                  !                                    ! alphaPE
-                  pab_pe(ji,jj,jk,jp_tem) = 0.
-                  pab_pe(ji,jj,jk,jp_sal) = 0.
-                  !
-               END DO
-            END DO
-         END DO
+         DO_3D_11_11( 1, jpkm1 )
+            zt  = pts(ji,jj,jk,jp_tem) - (-1._wp)  ! temperature anomaly (t-T0)
+            zs = pts (ji,jj,jk,jp_sal) - 34.2_wp   ! abs. salinity anomaly (s-S0)
+            zh  = gdept(ji,jj,jk,Kmm)                ! depth in meters  at t-point
+            ztm = tmask(ji,jj,jk)                  ! tmask
+            zn  = 0.5_wp * zh * r1_rho0 * ztm
+            !                                    ! Potential Energy
+            ppen(ji,jj,jk) = 0.
+            !                                    ! alphaPE
+            pab_pe(ji,jj,jk,jp_tem) = 0.
+            pab_pe(ji,jj,jk,jp_sal) = 0.
+            !
+         END_3D
          !
       CASE DEFAULT
@@ -1365,14 +1284,11 @@
       INTEGER  ::   ioptio   ! local integer
       !!
-      NAMELIST/nameos/ ln_TEOS10, ln_EOS80, ln_SEOS   , ln_LEOS, &
-         &             rn_a0    , rn_b0   , rn_lambda1, rn_mu1 , &
-         &                                  rn_lambda2, rn_mu2 , rn_nu
-      !!----------------------------------------------------------------------
-      !
-      REWIND( numnam_ref )              ! Namelist nameos in reference namelist : equation of state
+      NAMELIST/nameos/ ln_TEOS10, ln_EOS80, ln_SEOS, ln_LEOS, rn_a0, rn_b0, &
+         &             rn_lambda1, rn_mu1, rn_lambda2, rn_mu2, rn_nu
+      !!----------------------------------------------------------------------
+      !
       READ  ( numnam_ref, nameos, IOSTAT = ios, ERR = 901 )
 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nameos in reference namelist' )
       !
-      REWIND( numnam_cfg )              ! Namelist nameos in configuration namelist : equation of state
       READ  ( numnam_cfg, nameos, IOSTAT = ios, ERR = 902 )
 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nameos in configuration namelist' )
Index: NEMO/branches/2020/ticket_2459/tests/ISOMIP+/MY_SRC/isfcavgam.F90
===================================================================
--- NEMO/branches/2020/ticket_2459/tests/ISOMIP+/MY_SRC/isfcavgam.F90	(revision 12887)
+++ NEMO/branches/2020/ticket_2459/tests/ISOMIP+/MY_SRC/isfcavgam.F90	(revision 12889)
@@ -91,7 +91,7 @@
          pgs(:,:) = rn_gammas0
       CASE ( 'vel' ) ! gamma is proportional to u*
-         CALL gammats_vel      (                   zutbl, zvtbl, rCd0_top, r_ke0_top,               pgt, pgs )
+         CALL gammats_vel      (                   zutbl, zvtbl, rCd0_top, rn_vtide**2,               pgt, pgs )
       CASE ( 'vel_stab' ) ! gamma depends of stability of boundary layer and u*
-         CALL gammats_vel_stab (Kmm, pttbl, pstbl, zutbl, zvtbl, rCd0_top, r_ke0_top, pqoce, pqfwf, pgt, pgs )
+         CALL gammats_vel_stab (Kmm, pttbl, pstbl, zutbl, zvtbl, rCd0_top, rn_vtide**2, pqoce, pqfwf, pgt, pgs )
       CASE DEFAULT
          CALL ctl_stop('STOP','method to compute gamma (cn_gammablk) is unknown (should not see this)')
Index: NEMO/branches/2020/ticket_2459/tests/ISOMIP+/MY_SRC/isfstp.F90
===================================================================
--- NEMO/branches/2020/ticket_2459/tests/ISOMIP+/MY_SRC/isfstp.F90	(revision 12887)
+++ NEMO/branches/2020/ticket_2459/tests/ISOMIP+/MY_SRC/isfstp.F90	(revision 12889)
@@ -250,5 +250,5 @@
       IF ( l_isfoasis .AND. ln_isf ) THEN
          !
-         CALL ctl_stop( ' ln_ctl and ice shelf not tested' )
+         CALL ctl_stop( 'namelist combination ln_cpl and ln_isf not tested' )
          !
          ! NEMO coupled to ATMO model with isf cavity need oasis method for melt computation 
@@ -291,9 +291,7 @@
       !!----------------------------------------------------------------------
       !
-      REWIND( numnam_ref )              ! Namelist namsbc_rnf in reference namelist : Runoffs 
       READ  ( numnam_ref, namisf, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namisf in reference namelist' )
       !
-      REWIND( numnam_cfg )              ! Namelist namsbc_rnf in configuration namelist : Runoffs
       READ  ( numnam_cfg, namisf, IOSTAT = ios, ERR = 902 )
 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namisf in configuration namelist' )
Index: NEMO/branches/2020/ticket_2459/tests/ISOMIP+/MY_SRC/istate.F90
===================================================================
--- NEMO/branches/2020/ticket_2459/tests/ISOMIP+/MY_SRC/istate.F90	(revision 12887)
+++ NEMO/branches/2020/ticket_2459/tests/ISOMIP+/MY_SRC/istate.F90	(revision 12889)
@@ -41,4 +41,6 @@
    PUBLIC   istate_init   ! routine called by step.F90
 
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -75,5 +77,5 @@
       rhd  (:,:,:  ) = 0._wp   ;   rhop (:,:,:  ) = 0._wp      ! set one for all to 0 at level jpk
       rn2b (:,:,:  ) = 0._wp   ;   rn2  (:,:,:  ) = 0._wp      ! set one for all to 0 at levels 1 and jpk
-      ts   (:,:,:,:,Kaa) = 0._wp                               ! set one for all to 0 at level jpk
+      ts  (:,:,:,:,Kaa) = 0._wp                                   ! set one for all to 0 at level jpk
       rab_b(:,:,:,:) = 0._wp   ;   rab_n(:,:,:,:) = 0._wp      ! set one for all to 0 at level jpk
 #if defined key_agrif
@@ -90,5 +92,5 @@
          !                                    ! ---------------
          numror = 0                           ! define numror = 0 -> no restart file to read
-         neuler = 0                           ! Set time-step indicator at nit000 (euler forward)
+         l_1st_euler = .true.                 ! Set time-step indicator at nit000 (euler forward)
          CALL day_init                        ! model calendar (using both namelist and restart infos)
          !                                    ! Initialization of ocean to zero
@@ -103,11 +105,9 @@
                ! Apply minimum wetdepth criterion
                !
-               DO jj = 1,jpj
-                  DO ji = 1,jpi
-                     IF( ht_0(ji,jj) + ssh(ji,jj,Kbb)  < rn_wdmin1 ) THEN
-                        ssh(ji,jj,Kbb) = tmask(ji,jj,1)*( rn_wdmin1 - (ht_0(ji,jj)) )
-                     ENDIF
-                  END DO
-               END DO 
+               DO_2D_11_11
+                  IF( ht_0(ji,jj) + ssh(ji,jj,Kbb)  < rn_wdmin1 ) THEN
+                     ssh(ji,jj,Kbb) = tmask(ji,jj,1)*( rn_wdmin1 - (ht_0(ji,jj)) )
+                  ENDIF
+               END_2D
             ENDIF 
             uu  (:,:,:,Kbb) = 0._wp
@@ -159,15 +159,11 @@
       !
 !!gm  the use of umsak & vmask is not necessary below as uu(:,:,:,Kmm), vv(:,:,:,Kmm), uu(:,:,:,Kbb), vv(:,:,:,Kbb) are always masked
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               uu_b(ji,jj,Kmm) = uu_b(ji,jj,Kmm) + e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) * umask(ji,jj,jk)
-               vv_b(ji,jj,Kmm) = vv_b(ji,jj,Kmm) + e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm) * vmask(ji,jj,jk)
-               !
-               uu_b(ji,jj,Kbb) = uu_b(ji,jj,Kbb) + e3u(ji,jj,jk,Kbb) * uu(ji,jj,jk,Kbb) * umask(ji,jj,jk)
-               vv_b(ji,jj,Kbb) = vv_b(ji,jj,Kbb) + e3v(ji,jj,jk,Kbb) * vv(ji,jj,jk,Kbb) * vmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
+      DO_3D_11_11( 1, jpkm1 )
+         uu_b(ji,jj,Kmm) = uu_b(ji,jj,Kmm) + e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) * umask(ji,jj,jk)
+         vv_b(ji,jj,Kmm) = vv_b(ji,jj,Kmm) + e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm) * vmask(ji,jj,jk)
+         !
+         uu_b(ji,jj,Kbb) = uu_b(ji,jj,Kbb) + e3u(ji,jj,jk,Kbb) * uu(ji,jj,jk,Kbb) * umask(ji,jj,jk)
+         vv_b(ji,jj,Kbb) = vv_b(ji,jj,Kbb) + e3v(ji,jj,jk,Kbb) * vv(ji,jj,jk,Kbb) * vmask(ji,jj,jk)
+      END_3D
       !
       uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * r1_hu(:,:,Kmm)
Index: NEMO/branches/2020/ticket_2459/tests/ISOMIP+/MY_SRC/sbcfwb.F90
===================================================================
--- NEMO/branches/2020/ticket_2459/tests/ISOMIP+/MY_SRC/sbcfwb.F90	(revision 12887)
+++ NEMO/branches/2020/ticket_2459/tests/ISOMIP+/MY_SRC/sbcfwb.F90	(revision 12889)
@@ -151,5 +151,5 @@
          ENDIF   
          !                                         ! Update fwfold if new year start
-         ikty = 365 * 86400 / rn_Dt               !!bug  use of 365 days leap year or 360d year !!!!!!!
+         ikty = 365 * 86400 / rn_Dt                  !!bug  use of 365 days leap year or 360d year !!!!!!!
          IF( MOD( kt, ikty ) == 0 ) THEN
             a_fwb_b = a_fwb                           ! mean sea level taking into account the ice+snow
Index: NEMO/branches/2020/ticket_2459/tests/ISOMIP+/MY_SRC/tradmp.F90
===================================================================
--- NEMO/branches/2020/ticket_2459/tests/ISOMIP+/MY_SRC/tradmp.F90	(revision 12887)
+++ NEMO/branches/2020/ticket_2459/tests/ISOMIP+/MY_SRC/tradmp.F90	(revision 12889)
@@ -51,4 +51,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   resto    !: restoring coeff. on T and S (s-1)
 
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -110,41 +112,29 @@
       CASE( 0 )                        !*  newtonian damping throughout the water column  *!
          DO jn = 1, jpts
-            DO jk = 1, jpkm1
-               DO jj = 2, jpjm1
-                  DO ji = fs_2, fs_jpim1   ! vector opt.
-                     pts(ji,jj,jk,jn,Krhs) = pts(ji,jj,jk,jn,Krhs)           &
-                        &                  + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jn) - pts(ji,jj,jk,jn,Kbb) )
-                  END DO
-               END DO
-            END DO
+            DO_3D_00_00( 1, jpkm1 )
+               pts(ji,jj,jk,jn,Krhs) = pts(ji,jj,jk,jn,Krhs)           &
+                  &                  + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jn) - pts(ji,jj,jk,jn,Kbb) )
+            END_3D
          END DO
          !
       CASE ( 1 )                       !*  no damping in the turbocline (avt > 5 cm2/s)  *!
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  IF( avt(ji,jj,jk) <= avt_c ) THEN
-                     pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   &
-                        &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - pts(ji,jj,jk,jp_tem,Kbb) )
-                     pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs)   &
-                        &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - pts(ji,jj,jk,jp_sal,Kbb) )
-                  ENDIF
-               END DO
-            END DO
-         END DO
+         DO_3D_00_00( 1, jpkm1 )
+            IF( avt(ji,jj,jk) <= avt_c ) THEN
+               pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   &
+                  &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - pts(ji,jj,jk,jp_tem,Kbb) )
+               pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs)   &
+                  &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - pts(ji,jj,jk,jp_sal,Kbb) )
+            ENDIF
+         END_3D
          !
       CASE ( 2 )                       !*  no damping in the mixed layer   *!
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  IF( gdept(ji,jj,jk,Kmm) >= hmlp (ji,jj) ) THEN
-                     pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   &
-                        &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - pts(ji,jj,jk,jp_tem,Kbb) )
-                     pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs)   &
-                        &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - pts(ji,jj,jk,jp_sal,Kbb) )
-                  ENDIF
-               END DO
-            END DO
-         END DO
+         DO_3D_00_00( 1, jpkm1 )
+            IF( gdept(ji,jj,jk,Kmm) >= hmlp (ji,jj) ) THEN
+               pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   &
+                  &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - pts(ji,jj,jk,jp_tem,Kbb) )
+               pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs)   &
+                  &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - pts(ji,jj,jk,jp_sal,Kbb) )
+            ENDIF
+         END_3D
          !
       END SELECT
@@ -157,6 +147,6 @@
       ENDIF
       !                           ! Control print
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' dmp  - Ta: ', mask1=tmask,   &
-         &                       tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
+      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' dmp  - Ta: ', mask1=tmask,   &
+         &                                  tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
       !
       IF( ln_timing )   CALL timing_stop('tra_dmp')
@@ -178,9 +168,7 @@
       !!----------------------------------------------------------------------
       !
-      REWIND( numnam_ref )   ! Namelist namtra_dmp in reference namelist : T & S relaxation
       READ  ( numnam_ref, namtra_dmp, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_dmp in reference namelist' )
       !
-      REWIND( numnam_cfg )   ! Namelist namtra_dmp in configuration namelist : T & S relaxation
       READ  ( numnam_cfg, namtra_dmp, IOSTAT = ios, ERR = 902 )
 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtra_dmp in configuration namelist' )
