Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/INSTALL.rst
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/INSTALL.rst	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/INSTALL.rst	(revision 13877)
@@ -122,5 +122,5 @@
 .. code:: console
 
-   $ svn co https://forge.ipsl.jussieu.fr/nemo/svn/NEMO/trunk
+   $ svn co https://forge.ipsl.jussieu.fr/nemo/svn/NEMO/releases/r4.0/r4.0.4
 
 Description of 1\ :sup:`st` level tree structure
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/AGRIF_DEMO/EXPREF/file_def_nemo-ice.xml
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/AGRIF_DEMO/EXPREF/file_def_nemo-ice.xml	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/AGRIF_DEMO/EXPREF/file_def_nemo-ice.xml	(revision 13877)
@@ -52,7 +52,6 @@
        <field field_ref="normstr"          name="normstr" />
        <field field_ref="sheastr"          name="sheastr" />
-       <field field_ref="isig1"            name="isig1"   />
-       <field field_ref="isig2"            name="isig2"   />
-       <field field_ref="isig3"            name="isig3"   />
+       <field field_ref="sig1_pnorm"       name="sig1_pnorm"/>
+       <field field_ref="sig2_pnorm"       name="sig2_pnorm"/>
        
        <!-- heat fluxes -->
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/C1D_PAPA/EXPREF/namelist_cfg
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/C1D_PAPA/EXPREF/namelist_cfg	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/C1D_PAPA/EXPREF/namelist_cfg	(revision 13877)
@@ -428,7 +428,4 @@
 /
 !-----------------------------------------------------------------------
-!-----------------------------------------------------------------------
-/
-!-----------------------------------------------------------------------
 &namhsb        !  Heat and salt budgets                                 (default: OFF)
 !-----------------------------------------------------------------------
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/ORCA2_ICE_ABL/EXPREF/file_def_nemo-ice.xml
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/ORCA2_ICE_ABL/EXPREF/file_def_nemo-ice.xml	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/ORCA2_ICE_ABL/EXPREF/file_def_nemo-ice.xml	(revision 13877)
@@ -53,7 +53,6 @@
        <field field_ref="normstr"          name="normstr" />
        <field field_ref="sheastr"          name="sheastr" />
-       <field field_ref="isig1"            name="isig1"   />
-       <field field_ref="isig2"            name="isig2"   />
-       <field field_ref="isig3"            name="isig3"   />
+       <field field_ref="sig1_pnorm"       name="sig1_pnorm"/>
+       <field field_ref="sig2_pnorm"       name="sig2_pnorm"/>
        
        <!-- heat fluxes -->
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/ORCA2_ICE_PISCES/EXPREF/file_def_nemo-ice.xml
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/ORCA2_ICE_PISCES/EXPREF/file_def_nemo-ice.xml	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/ORCA2_ICE_PISCES/EXPREF/file_def_nemo-ice.xml	(revision 13877)
@@ -53,7 +53,6 @@
        <field field_ref="normstr"          name="normstr" />
        <field field_ref="sheastr"          name="sheastr" />
-       <field field_ref="isig1"            name="isig1"   />
-       <field field_ref="isig2"            name="isig2"   />
-       <field field_ref="isig3"            name="isig3"   />
+       <field field_ref="sig1_pnorm"       name="sig1_pnorm"/>
+       <field field_ref="sig2_pnorm"       name="sig2_pnorm"/>
        
        <!-- heat fluxes -->
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/ORCA2_OFF_PISCES/EXPREF/namelist_cfg
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/ORCA2_OFF_PISCES/EXPREF/namelist_cfg	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/ORCA2_OFF_PISCES/EXPREF/namelist_cfg	(revision 13877)
@@ -388,7 +388,4 @@
 /
 !-----------------------------------------------------------------------
-!-----------------------------------------------------------------------
-/
-!-----------------------------------------------------------------------
 &namhsb        !  Heat and salt budgets                                 (default: OFF)
 !-----------------------------------------------------------------------
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/ORCA2_OFF_TRC/EXPREF/namelist_cfg
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/ORCA2_OFF_TRC/EXPREF/namelist_cfg	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/ORCA2_OFF_TRC/EXPREF/namelist_cfg	(revision 13877)
@@ -384,7 +384,4 @@
 /
 !-----------------------------------------------------------------------
-!-----------------------------------------------------------------------
-/
-!-----------------------------------------------------------------------
 &namhsb        !  Heat and salt budgets                                 (default: OFF)
 !-----------------------------------------------------------------------
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/ORCA2_SAS_ICE/EXPREF/file_def_nemo-ice.xml
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/ORCA2_SAS_ICE/EXPREF/file_def_nemo-ice.xml	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/ORCA2_SAS_ICE/EXPREF/file_def_nemo-ice.xml	(revision 13877)
@@ -52,7 +52,6 @@
        <field field_ref="normstr"          name="normstr" />
        <field field_ref="sheastr"          name="sheastr" />
-       <field field_ref="isig1"            name="isig1"   />
-       <field field_ref="isig2"            name="isig2"   />
-       <field field_ref="isig3"            name="isig3"   />
+       <field field_ref="sig1_pnorm"       name="sig1_pnorm"/>
+       <field field_ref="sig2_pnorm"       name="sig2_pnorm"/>
        
        <!-- heat fluxes -->
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/SHARED/field_def_nemo-ice.xml
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/SHARED/field_def_nemo-ice.xml	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/SHARED/field_def_nemo-ice.xml	(revision 13877)
@@ -75,7 +75,6 @@
           <field id="utau_bi"      long_name="X-component of ocean bottom stress on sea ice -landfast" standard_name="ocean_bottom_upward_x_stress"              unit="N/m2" />
           <field id="vtau_bi"      long_name="Y-component of ocean bottom stress on sea ice -landfast" standard_name="ocean_bottom_upward_y_stress"              unit="N/m2" />
-          <field id="isig1"        long_name="1st principal stress component for EVP rhg"                                                                        unit=""     />
-          <field id="isig2"        long_name="2nd principal stress component for EVP rhg"                                                                        unit=""     />
-          <field id="isig3"        long_name="convergence measure for EVP rheology (must be around 1)"                                                           unit=""     />
+          <field id="sig1_pnorm"   long_name="P-normalized 1st principal stress component"                                                                       unit=""     />
+          <field id="sig2_pnorm"   long_name="P-normalized 2nd principal stress component"                                                                       unit=""     />
           <field id="normstr"      long_name="Average normal stress in sea ice"                        standard_name="average_normal_stress"                     unit="N/m"  />
           <field id="sheastr"      long_name="Maximum shear stress in sea ice"                         standard_name="maximum_shear_stress"                      unit="N/m"  />
@@ -165,5 +164,9 @@
 
           <!-- diags -->
-          <field id="hfxdhc"       long_name="Heat content variation in snow and ice (neg = ice cooling)"   unit="W/m2" />
+          <field id="hfxdhc"        long_name="Heat content variation in snow and ice (neg = ice cooling)"   unit="W/m2" />
+          <!-- available if ln_icediachk=T -->
+          <field id="icedrift_mass" long_name="Ice mass drift (conservation check)"   unit="kg/m2/s" />
+          <field id="icedrift_salt" long_name="Ice salt drift (conservation check)"   unit="kg/m2/s" />
+          <field id="icedrift_heat" long_name="Ice heat drift (conservation check)"   unit="W/m2"    />
 
 	  <!-- sbcssm variables -->
@@ -400,7 +403,6 @@
 	  <field field_ref="normstr"          name="normstr" />
 	  <field field_ref="sheastr"          name="sheastr" />
-	  <field field_ref="isig1"            name="isig1"   />
-	  <field field_ref="isig2"            name="isig2"   />
-	  <field field_ref="isig3"            name="isig3"   />
+	  <field field_ref="sig1_pnorm"       name="sig1_pnorm"/>
+	  <field field_ref="sig2_pnorm"       name="sig2_pnorm"/>
 	  
 	  <!-- heat fluxes -->
@@ -459,8 +461,4 @@
           <field field_ref="vfxthin"          name="vfxthin" />
 	
-	  <!-- diag error for negative ice volume after advection -->
-	  <field field_ref="iceneg_pres"      name="sineg_pres" />
-	  <field field_ref="iceneg_volu"      name="sineg_volu" />
-	  <field field_ref="iceneg_hfx"       name="sineg_hfx"  />
         </field_group>
 
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/SHARED/namelist_ref
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/SHARED/namelist_ref	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/SHARED/namelist_ref	(revision 13877)
@@ -1264,5 +1264,4 @@
 !!gm   ln_trdmld_instant = .false.         !  flag to diagnose trends of instantantaneous or mean ML T/S
 !!gm
-/
 !-----------------------------------------------------------------------
 &namhsb        !  Heat and salt budgets                                 (default: OFF)
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/SPITZ12/EXPREF/file_def_nemo-ice.xml
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/SPITZ12/EXPREF/file_def_nemo-ice.xml	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/SPITZ12/EXPREF/file_def_nemo-ice.xml	(revision 13877)
@@ -53,7 +53,6 @@
        <field field_ref="normstr"          name="normstr" />
        <field field_ref="sheastr"          name="sheastr" />
-       <field field_ref="isig1"            name="isig1"   />
-       <field field_ref="isig2"            name="isig2"   />
-       <field field_ref="isig3"            name="isig3"   />
+       <field field_ref="sig1_pnorm"       name="sig1_pnorm"/>
+       <field field_ref="sig2_pnorm"       name="sig2_pnorm"/>
        
        <!-- heat fluxes -->
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/WED025/EXPREF/file_def_nemo-ice.xml
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/WED025/EXPREF/file_def_nemo-ice.xml	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/WED025/EXPREF/file_def_nemo-ice.xml	(revision 13877)
@@ -50,7 +50,6 @@
        <field field_ref="normstr"          name="normstr" />
        <field field_ref="sheastr"          name="sheastr" />
-       <field field_ref="isig1"            name="isig1"   />
-       <field field_ref="isig2"            name="isig2"   />
-       <field field_ref="isig3"            name="isig3"   />
+       <field field_ref="sig1_pnorm"       name="sig1_pnorm"/>
+       <field field_ref="sig2_pnorm"       name="sig2_pnorm"/>
        
        <!-- heat fluxes -->
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/WED025/EXPREF/namelist_cfg
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/WED025/EXPREF/namelist_cfg	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/cfgs/WED025/EXPREF/namelist_cfg	(revision 13877)
@@ -583,7 +583,4 @@
 /
 !-----------------------------------------------------------------------
-!-----------------------------------------------------------------------
-/
-!-----------------------------------------------------------------------
 &namhsb        !  Heat and salt budgets                                 (default: OFF)
 !-----------------------------------------------------------------------
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/NEMO_manual_state.txt
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/NEMO_manual_state.txt	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/NEMO_manual_state.txt	(revision 13877)
@@ -39,5 +39,5 @@
          namdia: iiceprt jiceprt
     nam_diaharm: nit000_han nitend_han nstep_han tname(1) tname(2)
-         namdrg: ln_OFF
+         namdrg: ln_drg_OFF
      namdrg_bot: rn_Cd0 rn_Uc0 rn_Cdmax
      namdrg_top: rn_Cd0 rn_Uc0 rn_Cdmax
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/latex/NEMO/subfiles/chap_TRA.tex
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/latex/NEMO/subfiles/chap_TRA.tex	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/latex/NEMO/subfiles/chap_TRA.tex	(revision 13877)
@@ -1229,5 +1229,5 @@
 In the computer code, a density anomaly, $d_a = \rho / \rho_o - 1$, is computed,
 with $\rho_o$ a reference density.
-Called \textit{rau0} in the code,
+Called \textit{rho0} in the code,
 $\rho_o$ is set in \mdl{phycst} to a value of \texttt{1,026} $Kg/m^3$.
 This is a sensible choice for the reference density used in a Boussinesq ocean climate model,
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/latex/NEMO/subfiles/chap_ZDF.tex
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/latex/NEMO/subfiles/chap_ZDF.tex	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/latex/NEMO/subfiles/chap_ZDF.tex	(revision 13877)
@@ -1160,5 +1160,5 @@
 \]
 When \np[=.true.]{ln_lin}{ln\_lin}, the value of $r$ used is \np{rn_Uc0}{rn\_Uc0}*\np{rn_Cd0}{rn\_Cd0}.
-Setting \np[=.true.]{ln_OFF}{ln\_OFF} (and \forcode{ln_lin=.true.}) is equivalent to setting $r=0$ and leads to a free-slip boundary condition.
+Setting \np[=.true.]{ln_drg_OFF}{ln\_OFF} (and \forcode{ln_lin=.true.}) is equivalent to setting $r=0$ and leads to a free-slip boundary condition.
 
 These values are assigned in \mdl{zdfdrg}.
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/nambdy_dta
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/nambdy_dta	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/nambdy_dta	(revision 13877)
@@ -29,4 +29,5 @@
    bn_aip      = 'NOT USED'              ,         24.       , 'siapnd'  ,    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     ''
    bn_hip      = 'NOT USED'              ,         24.       , 'sihpnd'  ,    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     ''
+   bn_hil      = 'NOT USED'              ,         24.       , 'sihlid'  ,    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     ''
    ! if bn_t_i etc are "not used", then define arbitrary temperatures and salinity and ponds
    rn_ice_tem  = 270.         !  arbitrary temperature               of incoming sea ice
@@ -35,3 +36,4 @@
    rn_ice_apnd = 0.2          !       --   pond fraction = a_ip/a_i            --
    rn_ice_hpnd = 0.05         !       --   pond depth                          --
+   rn_ice_hlid = 0.0          !       --   pond lid depth                      --
 /
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namdia
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namdia	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namdia	(revision 13877)
@@ -8,5 +8,5 @@
    ln_icediahsb     = .false.         !  output the heat, mass & salt budgets (T) or not (F)
    ln_icectl        = .false.         !  ice points output for debug (T or F)
-   iiceprt          =  10             !  i-index for debug
-   jiceprt          =  10             !  j-index for debug
+      iiceprt       =  10             !     i-index for debug
+      jiceprt       =  10             !     j-index for debug
 /
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namdrg
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namdrg	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namdrg	(revision 13877)
@@ -2,5 +2,5 @@
 &namdrg        !   top/bottom drag coefficient                          (default: NO selection)
 !-----------------------------------------------------------------------
-   ln_OFF      = .false.   !  free-slip       : Cd = 0                  (F => fill namdrg_bot
+   ln_drg_OFF  = .false.   !  free-slip       : Cd = 0                  (F => fill namdrg_bot
    ln_lin      = .false.   !      linear  drag: Cd = Cd0 Uc0                   &   namdrg_top)
    ln_non_lin  = .false.   !  non-linear  drag: Cd = Cd0 |U|
@@ -8,3 +8,4 @@
    !
    ln_drgimp   = .true.    !  implicit top/bottom friction flag
+      ln_drgice_imp = .false. ! implicit ice-ocean drag
 /
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namdrg_bot
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namdrg_bot	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namdrg_bot	(revision 13877)
@@ -1,4 +1,4 @@
 !-----------------------------------------------------------------------
-&namdrg_bot    !   BOTTOM friction                                      (ln_OFF =F)
+&namdrg_bot    !   BOTTOM friction                                      (ln_drg_OFF =F)
 !-----------------------------------------------------------------------
    rn_Cd0      =  1.e-3    !  drag coefficient [-]
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namdrg_top
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namdrg_top	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namdrg_top	(revision 13877)
@@ -1,4 +1,4 @@
 !-----------------------------------------------------------------------
-&namdrg_top    !   TOP friction                                         (ln_OFF =F & ln_isfcav=T)
+&namdrg_top    !   TOP friction                                         (ln_drg_OFF =F & ln_isfcav=T)
 !-----------------------------------------------------------------------
    rn_Cd0      =  1.e-3    !  drag coefficient [-]
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namdyn
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namdyn	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namdyn	(revision 13877)
@@ -10,8 +10,8 @@
    rn_ishlat        =   2.            !  lbc : free slip (0) ; partial slip (0-2) ; no slip (2) ; strong slip (>2)
    ln_landfast_L16  = .false.         !  landfast: parameterization from Lemieux 2016
-      rn_depfra     =   0.125         !        fraction of ocean depth that ice must reach to initiate landfast
+      rn_lf_depfra  =   0.125         !        fraction of ocean depth that ice must reach to initiate landfast
                                       !          recommended range: [0.1 ; 0.25]
-      rn_icebfr     =  15.            !        maximum bottom stress per unit volume [N/m3]
-      rn_lfrelax    =   1.e-5         !        relaxation time scale to reach static friction [s-1]
-      rn_tensile    =   0.05          !        isotropic tensile strength [0-0.5??]
+      rn_lf_bfr     =  15.            !        maximum bottom stress per unit volume [N/m3]
+      rn_lf_relax   =   1.e-5         !        relaxation time scale to reach static friction [s-1]
+      rn_lf_tensile =   0.05          !        isotropic tensile strength [0-0.5??]
 /
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namdyn_rhg
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namdyn_rhg	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namdyn_rhg	(revision 13877)
@@ -9,3 +9,4 @@
       rn_relast     =   0.333         !     ratio of elastic timescale to ice time step: Telast = dt_ice * rn_relast 
                                       !        advised value: 1/3 (rn_nevp=120) or 1/9 (rn_nevp=300)
+   ln_rhg_chkcvg    = .false.         !  check convergence of rheology (outputs: file ice_cvg.nc & variable uice_cvg)
 /
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/nameos
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/nameos	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/nameos	(revision 13877)
@@ -7,5 +7,5 @@
                                  !
    !                     ! S-EOS coefficients (ln_seos=T):
-   !                             !  rd(T,S,Z)*rau0 = -a0*(1+.5*lambda*dT+mu*Z+nu*dS)*dT+b0*dS
+   !                             !  rd(T,S,Z)*rho0 = -a0*(1+.5*lambda*dT+mu*Z+nu*dS)*dT+b0*dS
    rn_a0       =  1.6550e-1      !  thermal expension coefficient
    rn_b0       =  7.6554e-1      !  saline  expension coefficient
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namini
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namini	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namini	(revision 13877)
@@ -3,5 +3,7 @@
 !------------------------------------------------------------------------------
    ln_iceini        = .true.          !  activate ice initialization (T) or not (F)
-   ln_iceini_file   = .false.         !  netcdf file provided for initialization (T) or not (F)
+   nn_iceini_file   =   0             !     0 = Initialise sea ice based on SSTs
+                                      !     1 = Initialise sea ice from single category netcdf file
+                                      !     2 = Initialise sea ice from multi category restart file
    rn_thres_sst     =   2.0           !  max temp. above Tfreeze with initial ice = (sst - tfreeze)
    rn_hti_ini_n     =   3.0           !  initial ice thickness       (m), North
@@ -23,5 +25,7 @@
    rn_hpd_ini_n     =   0.05          !  initial pond depth          (m), North
    rn_hpd_ini_s     =   0.05          !        "            "             South
-   ! -- for ln_iceini_file = T
+   rn_hld_ini_n     =   0.0           !  initial pond lid depth      (m), North
+   rn_hld_ini_s     =   0.0           !        "            "             South
+   ! -- for nn_iceini_file = 1
    sn_hti = 'Ice_initialization'    , -12 ,'hti'   ,  .false.  , .true., 'yearly'  , '' , '', ''
    sn_hts = 'Ice_initialization'    , -12 ,'hts'   ,  .false.  , .true., 'yearly'  , '' , '', ''
@@ -34,4 +38,5 @@
    sn_apd = 'NOT USED'              , -12 ,'apd'   ,  .false.  , .true., 'yearly'  , '' , '', ''
    sn_hpd = 'NOT USED'              , -12 ,'hpd'   ,  .false.  , .true., 'yearly'  , '' , '', ''
+   sn_hld = 'NOT USED'              , -12 ,'hld'   ,  .false.  , .true., 'yearly'  , '' , '', ''
    cn_dir='./'
 /
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namsbc_blk
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namsbc_blk	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namsbc_blk	(revision 13877)
@@ -35,8 +35,11 @@
    sn_tair     = 't_10.15JUNE2009_fill'       ,    6.        , 'T_10_MOD',   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , ''
    sn_humi     = 'q_10.15JUNE2009_fill'       ,    6.        , 'Q_10_MOD',   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , ''
-   sn_hpgi     = 'NONE'                       ,   24.        , 'uhpg'    ,   .false.   , .false., 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bicubic', 'UG'     , ''
-   sn_hpgj     = 'NONE'                       ,   24.        , 'vhpg'    ,   .false.   , .false., 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bicubic', 'VG'     , ''
    sn_prec     = 'ncar_precip.15JUNE2009_fill',   -1.        , 'PRC_MOD1',   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , ''
    sn_snow     = 'ncar_precip.15JUNE2009_fill',   -1.        , 'SNOW'    ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , ''
    sn_slp      = 'slp.15JUNE2009_fill'        ,    6.        , 'SLP'     ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , ''
+   sn_uoatm    = 'NOT USED'                   ,    6.        , 'UOATM'   ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , 'Uoceatm', ''
+   sn_voatm    = 'NOT USED'                   ,    6.        , 'VOATM'   ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , 'Voceatm', ''
+   sn_cc       = 'NOT USED'                   ,   24.        , 'CC'      ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , ''
+   sn_hpgi     = 'NOT USED'                   ,   24.        , 'uhpg'    ,   .false.   , .false., 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bicubic', 'UG'     , ''
+   sn_hpgj     = 'NOT USED'                   ,   24.        , 'vhpg'    ,   .false.   , .false., 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bicubic', 'VG'     , ''
 /
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namsbc_cpl
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namsbc_cpl	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namsbc_cpl	(revision 13877)
@@ -2,9 +2,9 @@
 &namsbc_cpl    !   coupled ocean/atmosphere model                       ("key_oasis3")
 !-----------------------------------------------------------------------
-   nn_cplmodel   =     1   !  Maximum number of models to/from which NEMO is potentially sending/receiving data
-   ln_usecplmask = .false. !  use a coupling mask file to merge data received from several models
-   !                       !   -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel)
-   nn_cats_cpl   =     5   !  Number of sea ice categories over which coupling is to be carried out (if not 1)
-
+   nn_cplmodel       =     1   !  Maximum number of models to/from which NEMO is potentially sending/receiving data
+   ln_usecplmask     = .false. !  use a coupling mask file to merge data received from several models
+   !                           !   -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel)
+   ln_scale_ice_flux = .false. !  use ice fluxes that are already "ice weighted" ( i.e. multiplied ice concentration)
+   nn_cats_cpl       =     5   !  Number of sea ice categories over which coupling is to be carried out (if not 1)
    !_____________!__________________________!____________!_____________!______________________!________!
    !             !        description       !  multiple  !    vector   !       vector         ! vector !
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namthd
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namthd	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namthd	(revision 13877)
@@ -6,3 +6,5 @@
    ln_icedO         = .true.          !  activate ice growth in open-water (T) or not (F)
    ln_icedS         = .true.          !  activate brine drainage (T) or not (F)
+   !                                   
+   ln_leadhfx       = .true.          !  heat in the leads is used to melt sea-ice before warming the ocean
 /
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namthd_pnd
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namthd_pnd	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namthd_pnd	(revision 13877)
@@ -2,9 +2,12 @@
 &namthd_pnd     !   Melt ponds
 !------------------------------------------------------------------------------
-   ln_pnd           = .false.         !  activate melt ponds or not
-     ln_pnd_H12     = .false.         !  activate evolutive melt ponds (from Holland et al 2012)
-     ln_pnd_CST     = .false.         !  activate constant  melt ponds
-       rn_apnd      =   0.2           !     prescribed pond fraction, at Tsu=0 degC
-       rn_hpnd      =   0.05          !     prescribed pond depth,    at Tsu=0 degC
-     ln_pnd_alb     = .false.         !  melt ponds affect albedo or not
+   ln_pnd            = .false.        !  activate melt ponds or not
+      ln_pnd_LEV     = .false.        !  level ice melt ponds (from Flocco et al 2007,2010 & Holland et al 2012)
+         rn_apnd_min =   0.15         !     minimum ice fraction that contributes to melt pond. range: 0.0 -- 0.15 ??
+         rn_apnd_max =   0.85         !     maximum ice fraction that contributes to melt pond. range: 0.7 -- 0.85 ??
+      ln_pnd_CST     = .false.        !  constant  melt ponds
+         rn_apnd     =   0.2          !     prescribed pond fraction, at Tsu=0 degC
+         rn_hpnd     =   0.05         !     prescribed pond depth,    at Tsu=0 degC
+      ln_pnd_lids    = .true.         !  frozen lids on top of the ponds (only for ln_pnd_LEV)
+      ln_pnd_alb     = .true.         !  effect of melt ponds on ice albedo
 /
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namthd_zdf
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namthd_zdf	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namthd_zdf	(revision 13877)
@@ -7,4 +7,8 @@
    rn_cnd_s         =   0.31          !  thermal conductivity of the snow (0.31 W/m/K, Maykut and Untersteiner, 1971)
                                       !     Obs: 0.1-0.5 (Lecomte et al, JAMES 2013)
-   rn_kappa_i       =   1.0           !  radiation attenuation coefficient in sea ice [1/m]
+   rn_kappa_i       =   1.0           !  radiation attenuation coefficient in sea ice                     [1/m]
+   rn_kappa_s       =  10.0           !  nn_qtrice = 0: radiation attenuation coefficient in snow         [1/m]
+   rn_kappa_smlt    =   7.0           !  nn_qtrice = 1: radiation attenuation coefficient in melting snow [1/m]
+   rn_kappa_sdry    =  10.0           !                 radiation attenuation coefficient in dry snow     [1/m]
+   ln_zdf_chkcvg    = .false.         !  check convergence of heat diffusion scheme (output variable: tice_cvg)
 /
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namzdf_gls
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namzdf_gls	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namzdf_gls	(revision 13877)
@@ -13,4 +13,9 @@
    nn_z0_met     =     2   !  Method for surface roughness computation (0/1/2/3)
    !                             ! =3 requires ln_wave=T
+   nn_z0_ice     =   1     !  attenutaion of surface wave breaking under ice
+   !                       !           = 0 no impact of ice cover
+   !                       !           = 1 roughness uses rn_hsri and is weigthed by 1-TANH(10*fr_i)
+   !                       !           = 2 roughness uses rn_hsri and is weighted by 1-fr_i
+   !                       !           = 3 roughness uses rn_hsri and is weighted by 1-MIN(1,4*fr_i)
    nn_bc_surf    =     1   !  surface condition (0/1=Dir/Neum)
    nn_bc_bot     =     1   !  bottom condition (0/1=Dir/Neum)
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namzdf_tke
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namzdf_tke	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/doc/namelists/namzdf_tke	(revision 13877)
@@ -15,5 +15,4 @@
    ln_mxl0     = .true.    !  surface mixing length scale = F(wind stress) (T) or not (F)
    rn_mxl0     =   0.04    !  surface  buoyancy lenght scale minimum value
-   ln_drg      = .false.   !  top/bottom friction added as boundary condition of TKE
    ln_lc       = .true.    !  Langmuir cell parameterisation (Axell 2002)
       rn_lc       =   0.15    !  coef. associated to Langmuir cells
@@ -26,4 +25,8 @@
                               !        = 0  constant 10 m length scale
                               !        = 1  0.5m at the equator to 30m poleward of 40 degrees
-      rn_eice     =   4       !  below sea ice: =0 ON ; =4 OFF when ice fraction > 1/4   
+   nn_eice     =   1       !  attenutaion of langmuir & surface wave breaking under ice
+   !                       !           = 0 no impact of ice cover on langmuir & surface wave breaking
+   !                       !           = 1 weigthed by 1-TANH(10*fr_i)
+   !                       !           = 2 weighted by 1-fr_i
+   !                       !           = 3 weighted by 1-MIN(1,4*fr_i)
 /
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/ice.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/ice.F90	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/ice.F90	(revision 13877)
@@ -391,4 +391,9 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_vice         !: ice volume variation   [m/s] 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_vsnw         !: snw volume variation   [m/s] 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_aice         !: ice conc.  variation   [s-1] 
+   !
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_adv_mass     !: advection of mass (kg/m2/s)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_adv_salt     !: advection of salt (g/m2/s)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_adv_heat     !: advection of heat (W/m2)
    !
    !!----------------------------------------------------------------------
@@ -493,7 +498,8 @@
       ! * Ice diagnostics
       ii = ii + 1
-      ALLOCATE( diag_trp_vi(jpi,jpj) , diag_trp_vs (jpi,jpj) , diag_trp_ei(jpi,jpj),   & 
-         &      diag_trp_es(jpi,jpj) , diag_trp_sv (jpi,jpj) , diag_heat  (jpi,jpj),   &
-         &      diag_sice  (jpi,jpj) , diag_vice   (jpi,jpj) , diag_vsnw  (jpi,jpj), STAT=ierr(ii) )
+      ALLOCATE( diag_trp_vi(jpi,jpj) , diag_trp_vs (jpi,jpj) , diag_trp_ei(jpi,jpj),                      & 
+         &      diag_trp_es(jpi,jpj) , diag_trp_sv (jpi,jpj) , diag_heat  (jpi,jpj),                      &
+         &      diag_sice  (jpi,jpj) , diag_vice   (jpi,jpj) , diag_vsnw  (jpi,jpj), diag_aice(jpi,jpj),  &
+         &      diag_adv_mass(jpi,jpj), diag_adv_salt(jpi,jpj), diag_adv_heat(jpi,jpj), STAT=ierr(ii) )
 
       ! * Ice conservation
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/ice1d.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/ice1d.F90	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/ice1d.F90	(revision 13877)
@@ -145,4 +145,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sst_1d
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sss_1d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   frq_m_1d
 
    ! convergence check
@@ -225,5 +226,5 @@
       !
       ii = ii + 1
-      ALLOCATE( sst_1d(jpij) , sss_1d(jpij) , STAT=ierr(ii) )
+      ALLOCATE( sst_1d(jpij) , sss_1d(jpij) , frq_m_1d(jpij) , STAT=ierr(ii) )
       !
       ii = ii + 1
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icecor.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icecor.F90	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icecor.F90	(revision 13877)
@@ -55,5 +55,4 @@
       INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
       REAL(wp) ::   zsal, zzc
-      REAL(wp), DIMENSION(jpi,jpj) ::   zafx   ! concentration trends diag
       !!----------------------------------------------------------------------
       ! controls
@@ -95,12 +94,16 @@
                zsal = sv_i(ji,jj,jl)
                sv_i(ji,jj,jl) = MIN(  MAX( rn_simin*v_i(ji,jj,jl) , sv_i(ji,jj,jl) ) , rn_simax*v_i(ji,jj,jl)  )
-               sfx_res(ji,jj) = sfx_res(ji,jj) - ( sv_i(ji,jj,jl) - zsal ) * zzc   ! associated salt flux
+               IF( kn /= 0 ) & ! no ice-ocean exchanges if kn=0 (for bdy for instance) otherwise conservation diags will fail
+                  &   sfx_res(ji,jj) = sfx_res(ji,jj) - ( sv_i(ji,jj,jl) - zsal ) * zzc   ! associated salt flux
             END_2D
          END DO
       ENDIF
-      !                             !-----------------------------------------------------
-      CALL ice_var_zapsmall         !  Zap small values                                  !
-      !                             !-----------------------------------------------------
 
+      IF( kn /= 0 ) THEN   ! no zapsmall if kn=0 (for bdy for instance) because we do not want ice-ocean exchanges (wfx,sfx,hfx)
+         !                                                              otherwise conservation diags will fail
+         !                          !-----------------------------------------------------
+         CALL ice_var_zapsmall      !  Zap small values                                  !
+         !                          !-----------------------------------------------------
+      ENDIF
       !                             !-----------------------------------------------------
       IF( kn == 2 ) THEN            !  Ice drift case: Corrections to avoid wrong values !
@@ -119,47 +122,4 @@
 #endif
       ENDIF
-
-      !                             !-----------------------------------------------------
-      SELECT CASE( kn )             !  Diagnostics                                       !
-      !                             !-----------------------------------------------------
-      CASE( 1 )                        !--- dyn trend diagnostics
-         !
-         IF( ln_icediachk .OR. iom_use('hfxdhc') ) THEN
-            diag_heat(:,:) = - SUM(SUM( e_i (:,:,1:nlay_i,:) - e_i_b (:,:,1:nlay_i,:), dim=4 ), dim=3 ) * r1_Dt_ice &      ! W.m-2
-               &             - SUM(SUM( e_s (:,:,1:nlay_s,:) - e_s_b (:,:,1:nlay_s,:), dim=4 ), dim=3 ) * r1_Dt_ice
-            diag_sice(:,:) =   SUM(     sv_i(:,:,:)          - sv_i_b(:,:,:)                  , dim=3 ) * r1_Dt_ice * rhoi
-            diag_vice(:,:) =   SUM(     v_i (:,:,:)          - v_i_b (:,:,:)                  , dim=3 ) * r1_Dt_ice * rhoi
-            diag_vsnw(:,:) =   SUM(     v_s (:,:,:)          - v_s_b (:,:,:)                  , dim=3 ) * r1_Dt_ice * rhos
-         ENDIF
-         !                       ! concentration tendency (dynamics)
-         IF( iom_use('afxdyn') .OR. iom_use('afxthd') .OR. iom_use('afxtot') ) THEN 
-            zafx(:,:) = SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice 
-            CALL iom_put( 'afxdyn' , zafx )
-         ENDIF
-         !
-      CASE( 2 )                        !--- thermo trend diagnostics & ice aging
-         !
-         oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rDt_ice   ! ice natural aging incrementation
-         !
-         IF( ln_icediachk .OR. iom_use('hfxdhc') ) THEN
-            diag_heat(:,:) = diag_heat(:,:) &
-               &             - SUM(SUM( e_i (:,:,1:nlay_i,:) - e_i_b (:,:,1:nlay_i,:), dim=4 ), dim=3 ) * r1_Dt_ice &
-               &             - SUM(SUM( e_s (:,:,1:nlay_s,:) - e_s_b (:,:,1:nlay_s,:), dim=4 ), dim=3 ) * r1_Dt_ice
-            diag_sice(:,:) = diag_sice(:,:) &
-               &             + SUM(     sv_i(:,:,:)          - sv_i_b(:,:,:)                  , dim=3 ) * r1_Dt_ice * rhoi
-            diag_vice(:,:) = diag_vice(:,:) &
-               &             + SUM(     v_i (:,:,:)          - v_i_b (:,:,:)                  , dim=3 ) * r1_Dt_ice * rhoi
-            diag_vsnw(:,:) = diag_vsnw(:,:) &
-               &             + SUM(     v_s (:,:,:)          - v_s_b (:,:,:)                  , dim=3 ) * r1_Dt_ice * rhos
-            CALL iom_put ( 'hfxdhc' , diag_heat ) 
-         ENDIF
-         !                       ! concentration tendency (total + thermo)
-         IF( iom_use('afxdyn') .OR. iom_use('afxthd') .OR. iom_use('afxtot') ) THEN 
-            zafx(:,:) = zafx(:,:) + SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice
-            CALL iom_put( 'afxthd' , SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice )
-            CALL iom_put( 'afxtot' , zafx )
-         ENDIF
-         !
-      END SELECT
       !
       ! controls
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icectl.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icectl.F90	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icectl.F90	(revision 13877)
@@ -43,4 +43,6 @@
    PUBLIC   ice_prt
    PUBLIC   ice_prt3D
+   PUBLIC   ice_drift_wri
+   PUBLIC   ice_drift_init
 
    ! thresold rates for conservation
@@ -49,4 +51,10 @@
    REAL(wp), PARAMETER ::   zchk_s   = 2.5e-6   ! g/m2/s  <=> 1e-6 m of ice per hour spuriously gained/lost (considering s=10g/kg)
    REAL(wp), PARAMETER ::   zchk_t   = 7.5e-2   ! W/m2    <=> 1e-6 m of ice per hour spuriously gained/lost (considering Lf=3e5J/kg)
+
+   ! for drift outputs
+   CHARACTER(LEN=50)   ::   clname="icedrift_diagnostics.ascii"   ! ascii filename
+   INTEGER             ::   numicedrift                           ! outfile unit
+   REAL(wp)            ::   rdiag_icemass, rdiag_icesalt, rdiag_iceheat 
+   REAL(wp)            ::   rdiag_adv_icemass, rdiag_adv_icesalt, rdiag_adv_iceheat 
    
    !! * Substitutions
@@ -132,6 +140,6 @@
 
          ! -- advection scheme is conservative? -- !
-         zvtrp = glob_sum( 'icectl', ( diag_trp_vi * rhoi + diag_trp_vs * rhos ) * e1e2t ) ! must be close to 0 (only for Prather)
-         zetrp = glob_sum( 'icectl', ( diag_trp_ei        + diag_trp_es        ) * e1e2t ) ! must be close to 0 (only for Prather)
+         zvtrp = glob_sum( 'icectl', diag_adv_mass * e1e2t )
+         zetrp = glob_sum( 'icectl', diag_adv_heat * e1e2t )
 
          ! ice area (+epsi10 to set a threshold > 0 when there is no ice) 
@@ -156,9 +164,8 @@
                &                   WRITE(numout,*)   cd_routine,' : violation a_i > amax      = ',zdiag_amax
             ! check if advection scheme is conservative
-            !    only check for Prather because Ultimate-Macho uses corrective fluxes (wfx etc)
-            !    so the formulation for conservation is different (and not coded) 
-            !    it does not mean UM is not conservative (it is checked with above prints) => update (09/2019): same for Prather now
-            !IF( ln_adv_Pra .AND. ABS(zvtrp) > zchk_m * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) &
-            !   &                   WRITE(numout,*)   cd_routine,' : violation adv scheme [kg] = ',zvtrp * rDt_ice
+            IF( ABS(zvtrp) > zchk_m * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) &
+               &                   WRITE(numout,*)   cd_routine,' : violation adv scheme [kg] = ',zvtrp * rdt_ice
+            IF( ABS(zetrp) > zchk_t * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) &
+               &                   WRITE(numout,*)   cd_routine,' : violation adv scheme [J]  = ',zetrp * rdt_ice
          ENDIF
          !
@@ -186,13 +193,16 @@
       ! water flux
       ! -- mass diag -- !
-      zdiag_mass = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t )
+      zdiag_mass = glob_sum( 'icectl', (  wfx_ice   + wfx_snw   + wfx_spr + wfx_sub &
+         &                              + diag_vice + diag_vsnw - diag_adv_mass ) * e1e2t )
 
       ! -- salt diag -- !
-      zdiag_salt = glob_sum( 'icectl', ( sfx + diag_sice ) * e1e2t )
+      zdiag_salt = glob_sum( 'icectl', ( sfx + diag_sice - diag_adv_salt ) * e1e2t )
 
       ! -- heat diag -- !
-      ! clem: not the good formulation
-      !!zdiag_heat  = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat + hfx_thd + hfx_dyn + hfx_res + hfx_sub + hfx_spr  &
-      !!   &                              ) * e1e2t )
+      zdiag_heat  = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat - diag_adv_heat ) * e1e2t )
+      ! equivalent to this:
+      !!zdiag_heat = glob_sum( 'icectl', ( -diag_heat + hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw &
+      !!   &                                          - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr &
+      !!   &                                          ) * e1e2t )
 
       ! ice area (+epsi10 to set a threshold > 0 when there is no ice) 
@@ -204,5 +214,6 @@
          IF( ABS(zdiag_salt) > zchk_s * rn_icechk_glo * zarea ) &
             &                   WRITE(numout,*) cd_routine,' : violation salt cons. [g]  = ',zdiag_salt * rDt_ice
-         !!IF( ABS(zdiag_heat) > zchk_t * rn_icechk_glo * zarea ) WRITE(numout,*) cd_routine,' : violation heat cons. [J]  = ',zdiag_heat * rDt_ice
+         IF( ABS(zdiag_heat) > zchk_t * rn_icechk_glo * zarea ) &
+            &                   WRITE(numout,*) cd_routine,' : violation heat cons. [J]  = ',zdiag_heat * rDt_ice
       ENDIF
       !
@@ -725,4 +736,133 @@
       
    END SUBROUTINE ice_prt3D
+
+
+   SUBROUTINE ice_drift_wri( kt )
+      !!-------------------------------------------------------------------
+      !!                     ***  ROUTINE ice_drift_wri ***
+      !!
+      !! ** Purpose : conservation of mass, salt and heat
+      !!              write the drift in a ascii file at each time step
+      !!              and the total run drifts
+      !!-------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ice time-step index
+      !
+      INTEGER  ::   ji, jj
+      REAL(wp) ::   zdiag_mass, zdiag_salt, zdiag_heat, zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat
+      REAL(wp), DIMENSION(jpi,jpj) ::   zdiag_mass2D, zdiag_salt2D, zdiag_heat2D
+      !!-------------------------------------------------------------------
+      !
+      IF( kt == nit000 .AND. lwp ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'ice_drift_wri: sea-ice drifts'
+         WRITE(numout,*) '~~~~~~~~~~~~~'
+      ENDIF
+      !
+      ! 2D budgets (must be close to 0)
+      IF( iom_use('icedrift_mass') .OR. iom_use('icedrift_salt') .OR. iom_use('icedrift_heat') ) THEN
+         DO_2D( 1, 1, 1, 1 )
+            zdiag_mass2D(ji,jj) =   wfx_ice(ji,jj)   + wfx_snw(ji,jj)   + wfx_spr(ji,jj) + wfx_sub(ji,jj) &
+               &                  + diag_vice(ji,jj) + diag_vsnw(ji,jj) - diag_adv_mass(ji,jj)
+            zdiag_salt2D(ji,jj) = sfx(ji,jj) + diag_sice(ji,jj) - diag_adv_salt(ji,jj)
+            zdiag_heat2D(ji,jj) = qt_oce_ai(ji,jj) - qt_atm_oi(ji,jj) + diag_heat(ji,jj) - diag_adv_heat(ji,jj)
+         END_2D
+         !
+         ! write outputs
+         CALL iom_put( 'icedrift_mass', zdiag_mass2D )
+         CALL iom_put( 'icedrift_salt', zdiag_salt2D )
+         CALL iom_put( 'icedrift_heat', zdiag_heat2D )
+      ENDIF
+
+      ! -- mass diag -- !
+      zdiag_mass     = glob_sum( 'icectl', (  wfx_ice   + wfx_snw   + wfx_spr + wfx_sub &
+         &                                  + diag_vice + diag_vsnw - diag_adv_mass ) * e1e2t ) * rdt_ice
+      zdiag_adv_mass = glob_sum( 'icectl', diag_adv_mass * e1e2t ) * rDt_ice
+
+      ! -- salt diag -- !
+      zdiag_salt     = glob_sum( 'icectl', ( sfx + diag_sice - diag_adv_salt ) * e1e2t ) * rdt_ice * 1.e-3
+      zdiag_adv_salt = glob_sum( 'icectl', diag_adv_salt * e1e2t ) * rDt_ice * 1.e-3
+
+      ! -- heat diag -- !
+      zdiag_heat     = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat - diag_adv_heat ) * e1e2t )
+      zdiag_adv_heat = glob_sum( 'icectl', diag_adv_heat * e1e2t )
+
+      !                    ! write out to file
+      IF( lwp ) THEN
+         ! check global drift (must be close to 0)
+         WRITE(numicedrift,FMT='(2x,i6,3x,a19,4x,f25.5)') kt, 'mass drift     [kg]', zdiag_mass
+         WRITE(numicedrift,FMT='(11x,     a19,4x,f25.5)')     'salt drift     [kg]', zdiag_salt
+         WRITE(numicedrift,FMT='(11x,     a19,4x,f25.5)')     'heat drift     [W] ', zdiag_heat
+         ! check drift from advection scheme (can be /=0 with bdy but not sure why)
+         WRITE(numicedrift,FMT='(11x,     a19,4x,f25.5)')     'mass drift adv [kg]', zdiag_adv_mass
+         WRITE(numicedrift,FMT='(11x,     a19,4x,f25.5)')     'salt drift adv [kg]', zdiag_adv_salt
+         WRITE(numicedrift,FMT='(11x,     a19,4x,f25.5)')     'heat drift adv [W] ', zdiag_adv_heat
+      ENDIF
+      !                    ! drifts
+      rdiag_icemass = rdiag_icemass + zdiag_mass
+      rdiag_icesalt = rdiag_icesalt + zdiag_salt
+      rdiag_iceheat = rdiag_iceheat + zdiag_heat
+      rdiag_adv_icemass = rdiag_adv_icemass + zdiag_adv_mass
+      rdiag_adv_icesalt = rdiag_adv_icesalt + zdiag_adv_salt
+      rdiag_adv_iceheat = rdiag_adv_iceheat + zdiag_adv_heat
+      !
+      !                    ! output drifts and close ascii file
+      IF( kt == nitend - nn_fsbc + 1 .AND. lwp ) THEN
+         ! to ascii file
+         WRITE(numicedrift,*) '******************************************'
+         WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run mass drift     [kg]', rdiag_icemass
+         WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run mass drift adv [kg]', rdiag_adv_icemass
+         WRITE(numicedrift,*) '******************************************'
+         WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run salt drift     [kg]', rdiag_icesalt
+         WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run salt drift adv [kg]', rdiag_adv_icesalt
+         WRITE(numicedrift,*) '******************************************'
+         WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run heat drift     [W] ', rdiag_iceheat
+         WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run heat drift adv [W] ', rdiag_adv_iceheat
+         CLOSE( numicedrift )
+         !
+         ! to ocean output
+         WRITE(numout,*)
+         WRITE(numout,*) 'ice_drift_wri: ice drifts information for the run '
+         WRITE(numout,*) '~~~~~~~~~~~~~'
+         ! check global drift (must be close to 0)
+         WRITE(numout,*) '   sea-ice mass drift     [kg] = ', rdiag_icemass
+         WRITE(numout,*) '   sea-ice salt drift     [kg] = ', rdiag_icesalt
+         WRITE(numout,*) '   sea-ice heat drift     [W]  = ', rdiag_iceheat
+         ! check drift from advection scheme (can be /=0 with bdy but not sure why)
+         WRITE(numout,*) '   sea-ice mass drift adv [kg] = ', rdiag_adv_icemass
+         WRITE(numout,*) '   sea-ice salt drift adv [kg] = ', rdiag_adv_icesalt
+         WRITE(numout,*) '   sea-ice heat drift adv [W]  = ', rdiag_adv_iceheat
+      ENDIF
+      !
+   END SUBROUTINE ice_drift_wri
+
+   SUBROUTINE ice_drift_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE ice_drift_init  ***
+      !!                   
+      !! ** Purpose :   create output file, initialise arrays
+      !!----------------------------------------------------------------------
+      !
+      IF( .NOT.ln_icediachk ) RETURN ! exit
+      !
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'ice_drift_init: Output ice drifts to ',TRIM(clname), ' file'
+         WRITE(numout,*) '~~~~~~~~~~~~~'
+         WRITE(numout,*)
+         !
+         ! create output ascii file
+         CALL ctl_opn( numicedrift, clname, 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', 1, numout, lwp, narea )
+         WRITE(numicedrift,*) 'Timestep  Drifts'
+         WRITE(numicedrift,*) '******************************************'
+      ENDIF
+      !
+      rdiag_icemass = 0._wp
+      rdiag_icesalt = 0._wp
+      rdiag_iceheat = 0._wp
+      rdiag_adv_icemass = 0._wp
+      rdiag_adv_icesalt = 0._wp
+      rdiag_adv_iceheat = 0._wp
+      !
+   END SUBROUTINE ice_drift_init
       
 #else
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icedyn_adv_pra.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icedyn_adv_pra.F90	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icedyn_adv_pra.F90	(revision 13877)
@@ -88,5 +88,5 @@
       INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices
       INTEGER  ::   icycle                  ! number of sub-timestep for the advection
-      REAL(wp) ::   zdt                     !   -      -
+      REAL(wp) ::   zdt, z1_dt              !   -      -
       REAL(wp), DIMENSION(1)                  ::   zcflprv, zcflnow   ! for global communication
       REAL(wp), DIMENSION(jpi,jpj)            ::   zati1, zati2
@@ -100,4 +100,6 @@
       REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) ::   z0es
       REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) ::   z0ei
+      !! diagnostics
+      REAL(wp), DIMENSION(jpi,jpj)            ::   zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat      
       !!----------------------------------------------------------------------
       !
@@ -109,24 +111,8 @@
       ELSEWHERE                      ; zs_i(:,:,:) = 0._wp
       END WHERE
-      DO jl = 1, jpl
-         DO_2D( 0, 0, 0, 0 )
-            zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj  ,jl), ph_ip(ji  ,jj+1,jl), &
-               &                                               ph_ip(ji-1,jj  ,jl), ph_ip(ji  ,jj-1,jl), &
-               &                                               ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), &
-               &                                               ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) )
-            zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj  ,jl), ph_i (ji  ,jj+1,jl), &
-               &                                               ph_i (ji-1,jj  ,jl), ph_i (ji  ,jj-1,jl), &
-               &                                               ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), &
-               &                                               ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) )
-            zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj  ,jl), ph_s (ji  ,jj+1,jl), &
-               &                                               ph_s (ji-1,jj  ,jl), ph_s (ji  ,jj-1,jl), &
-               &                                               ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), &
-               &                                               ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) )
-            zsi_max (ji,jj,jl) = MAX( epsi20, zs_i (ji,jj,jl), zs_i (ji+1,jj  ,jl), zs_i (ji  ,jj+1,jl), &
-               &                                               zs_i (ji-1,jj  ,jl), zs_i (ji  ,jj-1,jl), &
-               &                                               zs_i (ji+1,jj+1,jl), zs_i (ji-1,jj-1,jl), &
-               &                                               zs_i (ji+1,jj-1,jl), zs_i (ji-1,jj+1,jl) )
-         END_2D
-      END DO
+      CALL icemax3D( ph_i , zhi_max )
+      CALL icemax3D( ph_s , zhs_max )
+      CALL icemax3D( ph_ip, zhip_max)
+      CALL icemax3D( zs_i , zsi_max )
 #if defined key_mpi3
       CALL lbc_lnk_nc_multi( 'icedyn_adv_pra', zhi_max, 'T', 1._wp, zhs_max, 'T', 1._wp, zhip_max, 'T', 1._wp, zsi_max, 'T', 1._wp )
@@ -145,27 +131,13 @@
          ELSEWHERE                      ; ze_s(:,:,jk,:) = 0._wp
          END WHERE
-      END DO
-      DO jl = 1, jpl
-         DO_3D( 0, 0, 0, 0, 1, nlay_i )
-            zei_max(ji,jj,jk,jl) = MAX( epsi20, ze_i(ji,jj,jk,jl), ze_i(ji+1,jj  ,jk,jl), ze_i(ji  ,jj+1,jk,jl), &
-               &                                                   ze_i(ji-1,jj  ,jk,jl), ze_i(ji  ,jj-1,jk,jl), &
-               &                                                   ze_i(ji+1,jj+1,jk,jl), ze_i(ji-1,jj-1,jk,jl), &
-               &                                                   ze_i(ji+1,jj-1,jk,jl), ze_i(ji-1,jj+1,jk,jl) )
-         END_3D
-      END DO
-      DO jl = 1, jpl
-         DO_3D( 0, 0, 0, 0, 1, nlay_s )
-            zes_max(ji,jj,jk,jl) = MAX( epsi20, ze_s(ji,jj,jk,jl), ze_s(ji+1,jj  ,jk,jl), ze_s(ji  ,jj+1,jk,jl), &
-               &                                                   ze_s(ji-1,jj  ,jk,jl), ze_s(ji  ,jj-1,jk,jl), &
-               &                                                   ze_s(ji+1,jj+1,jk,jl), ze_s(ji-1,jj-1,jk,jl), &
-               &                                                   ze_s(ji+1,jj-1,jk,jl), ze_s(ji-1,jj+1,jk,jl) )
-         END_3D
-      END DO
+      END DO   
+      CALL icemax4D( ze_i , zei_max )
+      CALL icemax4D( ze_s , zes_max )
 #if defined key_mpi3
-      CALL lbc_lnk_nc_multi( 'icedyn_adv_pra', zei_max, 'T', 1. )
-      CALL lbc_lnk_nc_multi( 'icedyn_adv_pra', zes_max, 'T', 1. )
+      CALL lbc_lnk_nc_multi( 'icedyn_adv_pra', zei_max, 'T', 1._wp )
+      CALL lbc_lnk_nc_multi( 'icedyn_adv_pra', zes_max, 'T', 1._wp )
 #else
-      CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1. )
-      CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1. )
+      CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1._wp )
+      CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1._wp )
 #endif
       !
@@ -184,4 +156,5 @@
       ENDIF
       zdt = rDt_ice / REAL(icycle)
+      z1_dt = 1._wp / zdt
       
       ! --- transport --- !
@@ -190,4 +163,10 @@
 
       DO jt = 1, icycle
+
+         ! diagnostics
+         zdiag_adv_mass(:,:) =   SUM(  pv_i(:,:,:) , dim=3 ) * rhoi + SUM(  pv_s(:,:,:) , dim=3 ) * rhos
+         zdiag_adv_salt(:,:) =   SUM( psv_i(:,:,:) , dim=3 ) * rhoi
+         zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) &
+            &                  - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 )
 
          ! record at_i before advection (for open water)
@@ -288,6 +267,33 @@
                   CALL adv_x( zdt , zudy , 0._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl ) 
                ENDIF
-           ENDIF
-            !
+            ENDIF
+            !
+         ENDIF
+         
+         ! --- Lateral boundary conditions --- !
+         !     caution: for gradients (sx and sy) the sign changes
+         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ice , 'T', 1._wp, sxice , 'T', -1._wp, syice , 'T', -1._wp  & ! ice volume
+            &                                , sxxice, 'T', 1._wp, syyice, 'T',  1._wp, sxyice, 'T',  1._wp  &
+            &                                , z0snw , 'T', 1._wp, sxsn  , 'T', -1._wp, sysn  , 'T', -1._wp  & ! snw volume
+            &                                , sxxsn , 'T', 1._wp, syysn , 'T',  1._wp, sxysn , 'T',  1._wp  )
+         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0smi , 'T', 1._wp, sxsal , 'T', -1._wp, sysal , 'T', -1._wp  & ! ice salinity
+            &                                , sxxsal, 'T', 1._wp, syysal, 'T',  1._wp, sxysal, 'T',  1._wp  &
+            &                                , z0ai  , 'T', 1._wp, sxa   , 'T', -1._wp, sya   , 'T', -1._wp  & ! ice concentration
+            &                                , sxxa  , 'T', 1._wp, syya  , 'T',  1._wp, sxya  , 'T',  1._wp  )
+         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0oi  , 'T', 1._wp, sxage , 'T', -1._wp, syage , 'T', -1._wp  & ! ice age
+            &                                , sxxage, 'T', 1._wp, syyage, 'T',  1._wp, sxyage, 'T',  1._wp  )
+         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0es  , 'T', 1._wp, sxc0  , 'T', -1._wp, syc0  , 'T', -1._wp  & ! snw enthalpy
+            &                                , sxxc0 , 'T', 1._wp, syyc0 , 'T',  1._wp, sxyc0 , 'T',  1._wp  ) 
+         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ei  , 'T', 1._wp, sxe   , 'T', -1._wp, sye   , 'T', -1._wp  & ! ice enthalpy
+            &                                , sxxe  , 'T', 1._wp, syye  , 'T',  1._wp, sxye  , 'T',  1._wp  )
+         IF ( ln_pnd_LEV ) THEN
+            CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ap , 'T', 1._wp, sxap , 'T', -1._wp, syap , 'T', -1._wp  & ! melt pond fraction
+               &                                , sxxap, 'T', 1._wp, syyap, 'T',  1._wp, sxyap, 'T',  1._wp  &
+               &                                , z0vp , 'T', 1._wp, sxvp , 'T', -1._wp, syvp , 'T', -1._wp  & ! melt pond volume
+               &                                , sxxvp, 'T', 1._wp, syyvp, 'T',  1._wp, sxyvp, 'T',  1._wp  ) 
+            IF ( ln_pnd_lids ) THEN
+               CALL lbc_lnk_multi( 'icedyn_adv_pra', z0vl ,'T', 1._wp, sxvl ,'T', -1._wp, syvl ,'T', -1._wp  & ! melt pond lid volume
+                  &                                , sxxvl,'T', 1._wp, syyvl,'T',  1._wp, sxyvl,'T',  1._wp  ) 
+            ENDIF
          ENDIF
 
@@ -325,4 +331,13 @@
          CALL lbc_lnk( 'icedyn_adv_pra', pato_i, 'T',  1.0_wp )
 #endif
+         !
+         ! --- diagnostics --- !
+         diag_adv_mass(:,:) = diag_adv_mass(:,:) + (   SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos &
+            &                                        - zdiag_adv_mass(:,:) ) * z1_dt
+         diag_adv_salt(:,:) = diag_adv_salt(:,:) + (   SUM( psv_i(:,:,:) , dim=3 ) * rhoi &
+            &                                        - zdiag_adv_salt(:,:) ) * z1_dt
+         diag_adv_heat(:,:) = diag_adv_heat(:,:) + ( - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) &
+            &                                        - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) &
+            &                                        - zdiag_adv_heat(:,:) ) * z1_dt
          !
          ! --- Ensure non-negative fields --- !
@@ -363,11 +378,18 @@
       !! 
       INTEGER  ::   ji, jj, jl, jcat                     ! dummy loop indices
+      INTEGER  ::   jj0                                  ! dummy loop indices
       REAL(wp) ::   zs1max, zslpmax, ztemp               ! local scalars
       REAL(wp) ::   zs1new, zalf , zalfq , zbt           !   -      -
       REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !   -      -
+      REAL(wp) ::   zpsm, zps0
+      REAL(wp) ::   zpsx, zpsy, zpsxx, zpsyy, zpsxy
       REAL(wp), DIMENSION(jpi,jpj) ::   zf0 , zfx  , zfy   , zbet   ! 2D workspace
       REAL(wp), DIMENSION(jpi,jpj) ::   zfm , zfxx , zfyy  , zfxy   !  -      -
       REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q         !  -      -
       !-----------------------------------------------------------------------
+      ! in order to avoid lbc_lnk (communications):
+      !    jj loop must be 1:jpj   if adv_x is called first
+      !                and 2:jpj-1 if adv_x is called second
+      jj0 = NINT(pcrh)
       !
       jcat = SIZE( ps0 , 3 )   ! size of input arrays
@@ -376,134 +398,159 @@
          !
          ! Limitation of moments.                                           
-         DO_2D( 0, 0, 1, 1 )
-            !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                     
-            psm (ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 )
-            !
-            zslpmax = MAX( 0._wp, ps0(ji,jj,jl) )
-            zs1max  = 1.5 * zslpmax
-            zs1new  = MIN( zs1max, MAX( -zs1max, psx(ji,jj,jl) ) )
-            zs2new  = MIN(  2.0 * zslpmax - 0.3334 * ABS( zs1new ),      &
-               &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj,jl) )  )
-            rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
-
-            ps0 (ji,jj,jl) = zslpmax  
-            psx (ji,jj,jl) = zs1new         * rswitch
-            psxx(ji,jj,jl) = zs2new         * rswitch
-            psy (ji,jj,jl) = psy (ji,jj,jl) * rswitch
-            psyy(ji,jj,jl) = psyy(ji,jj,jl) * rswitch
-            psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch
-         END_2D
-
-         !  Calculate fluxes and moments between boxes i<-->i+1              
-         DO_2D( 0, 0, 1, 1 )                   !  Flux from i to i+1 WHEN u GT 0
-            zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) )
-            zalf         =  MAX( 0._wp, put(ji,jj) ) * pdt / psm(ji,jj,jl)
-            zalfq        =  zalf * zalf
-            zalf1        =  1.0 - zalf
-            zalf1q       =  zalf1 * zalf1
-            !
-            zfm (ji,jj)  =  zalf  *   psm (ji,jj,jl)
-            zf0 (ji,jj)  =  zalf  * ( ps0 (ji,jj,jl) + zalf1 * ( psx(ji,jj,jl) + (zalf1 - zalf) * psxx(ji,jj,jl) ) )
-            zfx (ji,jj)  =  zalfq * ( psx (ji,jj,jl) + 3.0 * zalf1 * psxx(ji,jj,jl) )
-            zfxx(ji,jj)  =  zalf  *   psxx(ji,jj,jl) * zalfq
-            zfy (ji,jj)  =  zalf  * ( psy (ji,jj,jl) + zalf1 * psxy(ji,jj,jl) )
-            zfxy(ji,jj)  =  zalfq *   psxy(ji,jj,jl)
-            zfyy(ji,jj)  =  zalf  *   psyy(ji,jj,jl)
-
-            !  Readjust moments remaining in the box.
-            psm (ji,jj,jl)  =  psm (ji,jj,jl) - zfm(ji,jj)
-            ps0 (ji,jj,jl)  =  ps0 (ji,jj,jl) - zf0(ji,jj)
-            psx (ji,jj,jl)  =  zalf1q * ( psx(ji,jj,jl) - 3.0 * zalf * psxx(ji,jj,jl) )
-            psxx(ji,jj,jl)  =  zalf1  * zalf1q * psxx(ji,jj,jl)
-            psy (ji,jj,jl)  =  psy (ji,jj,jl) - zfy(ji,jj)
-            psyy(ji,jj,jl)  =  psyy(ji,jj,jl) - zfyy(ji,jj)
-            psxy(ji,jj,jl)  =  zalf1q * psxy(ji,jj,jl)
-         END_2D
-
-         DO_2D( 0, 0, 1, 0 )                   !  Flux from i+1 to i when u LT 0.
-            zalf          = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl) 
-            zalg  (ji,jj) = zalf
-            zalfq         = zalf * zalf
-            zalf1         = 1.0 - zalf
-            zalg1 (ji,jj) = zalf1
-            zalf1q        = zalf1 * zalf1
-            zalg1q(ji,jj) = zalf1q
-            !
-            zfm   (ji,jj) = zfm (ji,jj) + zalf  *    psm (ji+1,jj,jl)
-            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * (  ps0 (ji+1,jj,jl) &
-               &                                   - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) )
-            zfx   (ji,jj) = zfx (ji,jj) + zalfq * (  psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) )
-            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *    psxx(ji+1,jj,jl) * zalfq
-            zfy   (ji,jj) = zfy (ji,jj) + zalf  * (  psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) )
-            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *    psxy(ji+1,jj,jl)
-            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *    psyy(ji+1,jj,jl)
-         END_2D
-
-         DO_2D( 0, 0, 0, 0 )                   !  Readjust moments remaining in the box.
-            zbt  =       zbet(ji-1,jj)
-            zbt1 = 1.0 - zbet(ji-1,jj)
-            !
-            psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji-1,jj) )
-            ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji-1,jj) )
-            psx (ji,jj,jl) = zalg1q(ji-1,jj) * ( psx(ji,jj,jl) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj,jl) )
-            psxx(ji,jj,jl) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj,jl)
-            psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) - zfy (ji-1,jj) )
-            psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) - zfyy(ji-1,jj) )
-            psxy(ji,jj,jl) = zalg1q(ji-1,jj) * psxy(ji,jj,jl)
-         END_2D
-
-         !   Put the temporary moments into appropriate neighboring boxes.    
-         DO_2D( 0, 0, 0, 0 )                   !   Flux from i to i+1 IF u GT 0.
-            zbt  =       zbet(ji-1,jj)
-            zbt1 = 1.0 - zbet(ji-1,jj)
-            psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj,jl)
-            zalf          = zbt * zfm(ji-1,jj) / psm(ji,jj,jl)
-            zalf1         = 1.0 - zalf
-            ztemp         = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji-1,jj)
-            !
-            ps0 (ji,jj,jl) =  zbt  * ( ps0(ji,jj,jl) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj,jl)
-            psx (ji,jj,jl) =  zbt  * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) + zbt1 * psx(ji,jj,jl)
-            psxx(ji,jj,jl) =  zbt  * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj,jl)                             &
-               &                     + 5.0 * ( zalf * zalf1 * ( psx (ji,jj,jl) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp )  ) &
-               &            + zbt1 * psxx(ji,jj,jl)
-            psxy(ji,jj,jl) =  zbt  * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj,jl)             &
-               &                     + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * psy(ji,jj,jl) ) )   &
-               &            + zbt1 * psxy(ji,jj,jl)
-            psy (ji,jj,jl) =  zbt  * ( psy (ji,jj,jl) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj,jl)
-            psyy(ji,jj,jl) =  zbt  * ( psyy(ji,jj,jl) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj,jl)
-         END_2D
-
-         DO_2D( 0, 0, 0, 0 )                   !  Flux from i+1 to i IF u LT 0.
-            zbt  =       zbet(ji,jj)
-            zbt1 = 1.0 - zbet(ji,jj)
-            psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) )
-            zalf          = zbt1 * zfm(ji,jj) / psm(ji,jj,jl)
-            zalf1         = 1.0 - zalf
-            ztemp         = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj)
-            !
-            ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) )
-            psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp )
-            psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) &
-               &                                           + 5.0 * ( zalf * zalf1 * ( - psx(ji,jj,jl) + zfx(ji,jj) )    &
-               &                                           + ( zalf1 - zalf ) * ztemp ) )
-            psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl)  &
-               &                                           + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj,jl) ) )
-            psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) + zfy (ji,jj) )
-            psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) + zfyy(ji,jj) )
-         END_2D
-
+         DO jj = Njs0 - jj0, Nje0 + jj0
+            
+            DO ji = Nis0 - 1, Nie0 + 1
+
+               zpsm  = psm (ji,jj,jl) ! optimization
+               zps0  = ps0 (ji,jj,jl)
+               zpsx  = psx (ji,jj,jl)
+               zpsxx = psxx(ji,jj,jl)
+               zpsy  = psy (ji,jj,jl)
+               zpsyy = psyy(ji,jj,jl)
+               zpsxy = psxy(ji,jj,jl)
+
+               !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                     
+               zpsm = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20 )
+               !
+               zslpmax = MAX( 0._wp, zps0 )
+               zs1max  = 1.5 * zslpmax
+               zs1new  = MIN( zs1max, MAX( -zs1max, zpsx ) )
+               zs2new  = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), MAX( ABS( zs1new ) - zslpmax, zpsxx ) )
+               rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
+
+               zps0  = zslpmax  
+               zpsx  = zs1new  * rswitch
+               zpsxx = zs2new  * rswitch
+               zpsy  = zpsy    * rswitch
+               zpsyy = zpsyy   * rswitch
+               zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch
+
+               !  Calculate fluxes and moments between boxes i<-->i+1              
+               !                                !  Flux from i to i+1 WHEN u GT 0 
+               zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) )
+               zalf         =  MAX( 0._wp, put(ji,jj) ) * pdt / zpsm
+               zalfq        =  zalf * zalf
+               zalf1        =  1.0 - zalf
+               zalf1q       =  zalf1 * zalf1
+               !
+               zfm (ji,jj)  =  zalf  *   zpsm 
+               zf0 (ji,jj)  =  zalf  * ( zps0  + zalf1 * ( zpsx + (zalf1 - zalf) * zpsxx ) )
+               zfx (ji,jj)  =  zalfq * ( zpsx  + 3.0 * zalf1 * zpsxx )
+               zfxx(ji,jj)  =  zalf  *   zpsxx * zalfq
+               zfy (ji,jj)  =  zalf  * ( zpsy  + zalf1 * zpsxy )
+               zfxy(ji,jj)  =  zalfq *   zpsxy
+               zfyy(ji,jj)  =  zalf  *   zpsyy
+
+               !                                !  Readjust moments remaining in the box.
+               zpsm  =  zpsm  - zfm(ji,jj)
+               zps0  =  zps0  - zf0(ji,jj)
+               zpsx  =  zalf1q * ( zpsx - 3.0 * zalf * zpsxx )
+               zpsxx =  zalf1  * zalf1q * zpsxx
+               zpsy  =  zpsy  - zfy (ji,jj)
+               zpsyy =  zpsyy - zfyy(ji,jj)
+               zpsxy =  zalf1q * zpsxy
+               !
+               psm (ji,jj,jl) = zpsm ! optimization
+               ps0 (ji,jj,jl) = zps0 
+               psx (ji,jj,jl) = zpsx 
+               psxx(ji,jj,jl) = zpsxx
+               psy (ji,jj,jl) = zpsy 
+               psyy(ji,jj,jl) = zpsyy
+               psxy(ji,jj,jl) = zpsxy
+               !
+            END DO
+            
+            DO ji = Nis0 - 1, Nie0
+               !                                !  Flux from i+1 to i when u LT 0.
+               zalf          = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl) 
+               zalg  (ji,jj) = zalf
+               zalfq         = zalf * zalf
+               zalf1         = 1.0 - zalf
+               zalg1 (ji,jj) = zalf1
+               zalf1q        = zalf1 * zalf1
+               zalg1q(ji,jj) = zalf1q
+               !
+               zfm   (ji,jj) = zfm (ji,jj) + zalf  *    psm (ji+1,jj,jl)
+               zf0   (ji,jj) = zf0 (ji,jj) + zalf  * (  ps0 (ji+1,jj,jl) &
+                  &                                   - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) )
+               zfx   (ji,jj) = zfx (ji,jj) + zalfq * (  psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) )
+               zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *    psxx(ji+1,jj,jl) * zalfq
+               zfy   (ji,jj) = zfy (ji,jj) + zalf  * (  psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) )
+               zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *    psxy(ji+1,jj,jl)
+               zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *    psyy(ji+1,jj,jl)
+            END DO
+
+            DO ji = Nis0, Nie0
+               !
+               zpsm  = psm (ji,jj,jl) ! optimization
+               zps0  = ps0 (ji,jj,jl)
+               zpsx  = psx (ji,jj,jl)
+               zpsxx = psxx(ji,jj,jl)
+               zpsy  = psy (ji,jj,jl)
+               zpsyy = psyy(ji,jj,jl)
+               zpsxy = psxy(ji,jj,jl)
+               !                                !  Readjust moments remaining in the box.
+               zbt  =       zbet(ji-1,jj)
+               zbt1 = 1.0 - zbet(ji-1,jj)
+               !
+               zpsm  = zbt * zpsm + zbt1 * ( zpsm - zfm(ji-1,jj) )
+               zps0  = zbt * zps0 + zbt1 * ( zps0 - zf0(ji-1,jj) )
+               zpsx  = zalg1q(ji-1,jj) * ( zpsx + 3.0 * zalg(ji-1,jj) * zpsxx )
+               zpsxx = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * zpsxx
+               zpsy  = zbt * zpsy  + zbt1 * ( zpsy  - zfy (ji-1,jj) )
+               zpsyy = zbt * zpsyy + zbt1 * ( zpsyy - zfyy(ji-1,jj) )
+               zpsxy = zalg1q(ji-1,jj) * zpsxy
+
+               !   Put the temporary moments into appropriate neighboring boxes.    
+               !                                !   Flux from i to i+1 IF u GT 0.
+               zbt   =       zbet(ji-1,jj)
+               zbt1  = 1.0 - zbet(ji-1,jj)
+               zpsm  = zbt * ( zpsm + zfm(ji-1,jj) ) + zbt1 * zpsm
+               zalf  = zbt * zfm(ji-1,jj) / zpsm
+               zalf1 = 1.0 - zalf
+               ztemp = zalf * zps0 - zalf1 * zf0(ji-1,jj)
+               !
+               zps0  =  zbt  * ( zps0 + zf0(ji-1,jj) ) + zbt1 * zps0
+               zpsx  =  zbt  * ( zalf * zfx(ji-1,jj) + zalf1 * zpsx + 3.0 * ztemp ) + zbt1 * zpsx
+               zpsxx =  zbt  * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * zpsxx                            &
+                  &            + 5.0 * ( zalf * zalf1 * ( zpsx  - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) &
+                  &            + zbt1 * zpsxx
+               zpsxy =  zbt  * ( zalf * zfxy(ji-1,jj) + zalf1 * zpsxy            &
+                  &            + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * zpsy ) )  &
+                  &            + zbt1 * zpsxy
+               zpsy  =  zbt  * ( zpsy  + zfy (ji-1,jj) ) + zbt1 * zpsy 
+               zpsyy =  zbt  * ( zpsyy + zfyy(ji-1,jj) ) + zbt1 * zpsyy
+
+               !                                !  Flux from i+1 to i IF u LT 0.
+               zbt   =       zbet(ji,jj)
+               zbt1  = 1.0 - zbet(ji,jj)
+               zpsm  = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) )
+               zalf  = zbt1 * zfm(ji,jj) / zpsm
+               zalf1 = 1.0 - zalf
+               ztemp = - zalf * zps0 + zalf1 * zf0(ji,jj)
+               !
+               zps0  = zbt * zps0  + zbt1 * ( zps0 + zf0(ji,jj) )
+               zpsx  = zbt * zpsx  + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * zpsx + 3.0 * ztemp )
+               zpsxx = zbt * zpsxx + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * zpsxx &
+                  &                         + 5.0 * ( zalf * zalf1 * ( - zpsx + zfx(ji,jj) )    &
+                  &                         + ( zalf1 - zalf ) * ztemp ) )
+               zpsxy = zbt * zpsxy + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * zpsxy  &
+                  &                         + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * zpsy ) )
+               zpsy  = zbt * zpsy  + zbt1 * ( zpsy  + zfy (ji,jj) )
+               zpsyy = zbt * zpsyy + zbt1 * ( zpsyy + zfyy(ji,jj) )
+               !
+               psm (ji,jj,jl) = zpsm  ! optimization
+               ps0 (ji,jj,jl) = zps0 
+               psx (ji,jj,jl) = zpsx 
+               psxx(ji,jj,jl) = zpsxx
+               psy (ji,jj,jl) = zpsy 
+               psyy(ji,jj,jl) = zpsyy
+               psxy(ji,jj,jl) = zpsxy
+            END DO
+            !
+         END DO
+         !
       END DO
-
-      !-- Lateral boundary conditions
-#if defined key_mpi3
-      CALL lbc_lnk_nc_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T',  1.0_wp, ps0 , 'T',  1.0_wp   &
-         &                                , psx             , 'T', -1.0_wp, psy , 'T', -1.0_wp   &   ! caution gradient ==> the sign changes
-         &                                , psxx            , 'T',  1.0_wp, psyy, 'T',  1.0_wp , psxy, 'T',  1.0_wp )
-#else
-      CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T',  1.0_wp, ps0 , 'T',  1.0_wp   &
-         &                                , psx             , 'T', -1.0_wp, psy , 'T', -1.0_wp   &   ! caution gradient ==> the sign changes
-         &                                , psxx            , 'T',  1.0_wp, psyy, 'T',  1.0_wp , psxy, 'T',  1.0_wp )
-#endif
-      !
+      !      
    END SUBROUTINE adv_x
 
@@ -526,11 +573,18 @@
       !!
       INTEGER  ::   ji, jj, jl, jcat                     ! dummy loop indices
+      INTEGER  ::   ji0                                  ! dummy loop indices
       REAL(wp) ::   zs1max, zslpmax, ztemp               ! temporary scalars
       REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         -
       REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         -
+      REAL(wp) ::   zpsm, zps0
+      REAL(wp) ::   zpsx, zpsy, zpsxx, zpsyy, zpsxy
       REAL(wp), DIMENSION(jpi,jpj) ::   zf0, zfx , zfy , zbet   ! 2D workspace
       REAL(wp), DIMENSION(jpi,jpj) ::   zfm, zfxx, zfyy, zfxy   !  -      -
       REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q     !  -      -
       !---------------------------------------------------------------------
+      ! in order to avoid lbc_lnk (communications):
+      !    ji loop must be 1:jpi   if adv_y is called first
+      !                and 2:jpi-1 if adv_y is called second
+      ji0 = NINT(pcrh)
       !
       jcat = SIZE( ps0 , 3 )   ! size of input arrays
@@ -539,50 +593,66 @@
          !
          ! Limitation of moments.
-         DO_2D( 1, 1, 0, 0 )
-            !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
-            psm(ji,jj,jl) = MAX(  pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20  )
-            !
-            zslpmax = MAX( 0._wp, ps0(ji,jj,jl) )
+         DO_2D( 1, 1, ji0, ji0 )
+            !
+            zpsm  = psm (ji,jj,jl) ! optimization
+            zps0  = ps0 (ji,jj,jl)
+            zpsx  = psx (ji,jj,jl)
+            zpsxx = psxx(ji,jj,jl)
+            zpsy  = psy (ji,jj,jl)
+            zpsyy = psyy(ji,jj,jl)
+            zpsxy = psxy(ji,jj,jl)
+            !
+            !  Initialize volumes of boxes (=area if adv_y first called, =psm otherwise)
+            zpsm = MAX(  pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20  )
+            !
+            zslpmax = MAX( 0._wp, zps0 )
             zs1max  = 1.5 * zslpmax
-            zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj,jl) ) )
-            zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   &
-               &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj,jl) )  )
+            zs1new  = MIN( zs1max, MAX( -zs1max, zpsy ) )
+            zs2new  = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), MAX( ABS( zs1new )-zslpmax, zpsyy ) )
             rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
             !
-            ps0 (ji,jj,jl) = zslpmax  
-            psx (ji,jj,jl) = psx (ji,jj,jl) * rswitch
-            psxx(ji,jj,jl) = psxx(ji,jj,jl) * rswitch
-            psy (ji,jj,jl) = zs1new         * rswitch
-            psyy(ji,jj,jl) = zs2new         * rswitch
-            psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch
-         END_2D
- 
-         !  Calculate fluxes and moments between boxes j<-->j+1              
-         DO_2D( 1, 1, 0, 0 )                   !  Flux from j to j+1 WHEN v GT 0
+            zps0  = zslpmax  
+            zpsx  = zpsx  * rswitch
+            zpsxx = zpsxx * rswitch
+            zpsy  = zs1new         * rswitch
+            zpsyy = zs2new         * rswitch
+            zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch
+
+            !  Calculate fluxes and moments between boxes j<-->j+1              
+            !                                !  Flux from j to j+1 WHEN v GT 0   
             zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) )
-            zalf         =  MAX( 0._wp, pvt(ji,jj) ) * pdt / psm(ji,jj,jl)
+            zalf         =  MAX( 0._wp, pvt(ji,jj) ) * pdt / zpsm
             zalfq        =  zalf * zalf
             zalf1        =  1.0 - zalf
             zalf1q       =  zalf1 * zalf1
             !
-            zfm (ji,jj)  =  zalf  * psm(ji,jj,jl)
-            zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj,jl) + zalf1 * ( psy(ji,jj,jl)  + (zalf1-zalf) * psyy(ji,jj,jl) ) ) 
-            zfy (ji,jj)  =  zalfq *( psy(ji,jj,jl) + 3.0*zalf1*psyy(ji,jj,jl) )
-            zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj,jl)
-            zfx (ji,jj)  =  zalf  * ( psx(ji,jj,jl) + zalf1 * psxy(ji,jj,jl) )
-            zfxy(ji,jj)  =  zalfq * psxy(ji,jj,jl)
-            zfxx(ji,jj)  =  zalf  * psxx(ji,jj,jl)
-            !
-            !  Readjust moments remaining in the box.
-            psm (ji,jj,jl)  =  psm (ji,jj,jl) - zfm(ji,jj)
-            ps0 (ji,jj,jl)  =  ps0 (ji,jj,jl) - zf0(ji,jj)
-            psy (ji,jj,jl)  =  zalf1q * ( psy(ji,jj,jl) -3.0 * zalf * psyy(ji,jj,jl) )
-            psyy(ji,jj,jl)  =  zalf1 * zalf1q * psyy(ji,jj,jl)
-            psx (ji,jj,jl)  =  psx (ji,jj,jl) - zfx(ji,jj)
-            psxx(ji,jj,jl)  =  psxx(ji,jj,jl) - zfxx(ji,jj)
-            psxy(ji,jj,jl)  =  zalf1q * psxy(ji,jj,jl)
+            zfm (ji,jj)  =  zalf  * zpsm
+            zf0 (ji,jj)  =  zalf  * ( zps0 + zalf1 * ( zpsy  + (zalf1-zalf) * zpsyy ) ) 
+            zfy (ji,jj)  =  zalfq *( zpsy + 3.0*zalf1*zpsyy )
+            zfyy(ji,jj)  =  zalf  * zalfq * zpsyy
+            zfx (ji,jj)  =  zalf  * ( zpsx + zalf1 * zpsxy )
+            zfxy(ji,jj)  =  zalfq * zpsxy
+            zfxx(ji,jj)  =  zalf  * zpsxx
+            !
+            !                                !  Readjust moments remaining in the box.
+            zpsm   =  zpsm  - zfm(ji,jj)
+            zps0   =  zps0  - zf0(ji,jj)
+            zpsy   =  zalf1q * ( zpsy -3.0 * zalf * zpsyy )
+            zpsyy  =  zalf1 * zalf1q * zpsyy
+            zpsx   =  zpsx  - zfx(ji,jj)
+            zpsxx  =  zpsxx - zfxx(ji,jj)
+            zpsxy  =  zalf1q * zpsxy
+            !
+            psm (ji,jj,jl) = zpsm ! optimization
+            ps0 (ji,jj,jl) = zps0 
+            psx (ji,jj,jl) = zpsx 
+            psxx(ji,jj,jl) = zpsxx
+            psy (ji,jj,jl) = zpsy 
+            psyy(ji,jj,jl) = zpsyy
+            psxy(ji,jj,jl) = zpsxy
          END_2D
          !
-         DO_2D( 1, 0, 0, 0 )                   !  Flux from j+1 to j when v LT 0.
+         DO_2D( 1, 0, ji0, ji0 )
+            !                                !  Flux from j+1 to j when v LT 0.
             zalf          = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1,jl) 
             zalg  (ji,jj) = zalf
@@ -603,71 +673,74 @@
          END_2D
 
-         !  Readjust moments remaining in the box. 
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D( 0, 0, ji0, ji0 )
+            !                                !  Readjust moments remaining in the box.
             zbt  =         zbet(ji,jj-1)
             zbt1 = ( 1.0 - zbet(ji,jj-1) )
             !
-            psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji,jj-1) )
-            ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji,jj-1) )
-            psy (ji,jj,jl) = zalg1q(ji,jj-1) * ( psy(ji,jj,jl) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj,jl) )
-            psyy(ji,jj,jl) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj,jl)
-            psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) - zfx (ji,jj-1) )
-            psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) - zfxx(ji,jj-1) )
-            psxy(ji,jj,jl) = zalg1q(ji,jj-1) * psxy(ji,jj,jl)
+            zpsm  = psm (ji,jj,jl) ! optimization
+            zps0  = ps0 (ji,jj,jl)
+            zpsx  = psx (ji,jj,jl)
+            zpsxx = psxx(ji,jj,jl)
+            zpsy  = psy (ji,jj,jl)
+            zpsyy = psyy(ji,jj,jl)
+            zpsxy = psxy(ji,jj,jl)
+            !
+            zpsm  = zbt * zpsm + zbt1 * ( zpsm - zfm(ji,jj-1) )
+            zps0  = zbt * zps0 + zbt1 * ( zps0 - zf0(ji,jj-1) )
+            zpsy  = zalg1q(ji,jj-1) * ( zpsy + 3.0 * zalg(ji,jj-1) * zpsyy )
+            zpsyy = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * zpsyy
+            zpsx  = zbt * zpsx  + zbt1 * ( zpsx  - zfx (ji,jj-1) )
+            zpsxx = zbt * zpsxx + zbt1 * ( zpsxx - zfxx(ji,jj-1) )
+            zpsxy = zalg1q(ji,jj-1) * zpsxy
+
+            !   Put the temporary moments into appropriate neighboring boxes.    
+            !                                !   Flux from j to j+1 IF v GT 0.
+            zbt   =       zbet(ji,jj-1)
+            zbt1  = 1.0 - zbet(ji,jj-1)
+            zpsm  = zbt * ( zpsm + zfm(ji,jj-1) ) + zbt1 * zpsm 
+            zalf  = zbt * zfm(ji,jj-1) / zpsm 
+            zalf1 = 1.0 - zalf
+            ztemp = zalf * zps0 - zalf1 * zf0(ji,jj-1)
+            !
+            zps0  =   zbt  * ( zps0 + zf0(ji,jj-1) ) + zbt1 * zps0
+            zpsy  =   zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * zpsy + 3.0 * ztemp )  &
+               &             + zbt1 * zpsy  
+            zpsyy =   zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * zpsyy                           &
+               &             + 5.0 * ( zalf * zalf1 * ( zpsy - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) & 
+               &             + zbt1 * zpsyy
+            zpsxy =   zbt  * ( zalf * zfxy(ji,jj-1) + zalf1 * zpsxy             &
+               &             + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * zpsx ) )  &
+               &             + zbt1 * zpsxy
+            zpsx  =   zbt * ( zpsx  + zfx (ji,jj-1) ) + zbt1 * zpsx 
+            zpsxx =   zbt * ( zpsxx + zfxx(ji,jj-1) ) + zbt1 * zpsxx
+
+            !                                !  Flux from j+1 to j IF v LT 0.
+            zbt   =       zbet(ji,jj)
+            zbt1  = 1.0 - zbet(ji,jj)
+            zpsm  = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) )
+            zalf  = zbt1 * zfm(ji,jj) / zpsm
+            zalf1 = 1.0 - zalf
+            ztemp = - zalf * zps0 + zalf1 * zf0(ji,jj)
+            !
+            zps0  = zbt * zps0  + zbt1 * (  zps0 + zf0(ji,jj) )
+            zpsy  = zbt * zpsy  + zbt1 * (  zalf * zfy(ji,jj) + zalf1 * zpsy + 3.0 * ztemp )
+            zpsyy = zbt * zpsyy + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * zpsyy &
+               &                         + 5.0 * ( zalf * zalf1 * ( - zpsy + zfy(ji,jj) )     &
+               &                         + ( zalf1 - zalf ) * ztemp ) )
+            zpsxy = zbt * zpsxy + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * zpsxy  &
+               &                         + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * zpsx ) )
+            zpsx  = zbt * zpsx  + zbt1 * ( zpsx  + zfx (ji,jj) )
+            zpsxx = zbt * zpsxx + zbt1 * ( zpsxx + zfxx(ji,jj) )
+            !
+            psm (ji,jj,jl) = zpsm ! optimization
+            ps0 (ji,jj,jl) = zps0 
+            psx (ji,jj,jl) = zpsx 
+            psxx(ji,jj,jl) = zpsxx
+            psy (ji,jj,jl) = zpsy 
+            psyy(ji,jj,jl) = zpsyy
+            psxy(ji,jj,jl) = zpsxy
          END_2D
-
-         !   Put the temporary moments into appropriate neighboring boxes.    
-         DO_2D( 0, 0, 0, 0 )                   !  Flux from j to j+1 IF v GT 0.
-            zbt  =       zbet(ji,jj-1)
-            zbt1 = 1.0 - zbet(ji,jj-1)
-            psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj,jl) 
-            zalf          = zbt * zfm(ji,jj-1) / psm(ji,jj,jl) 
-            zalf1         = 1.0 - zalf
-            ztemp         = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji,jj-1)
-            !
-            ps0(ji,jj,jl)  =   zbt  * ( ps0(ji,jj,jl) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj,jl)
-            psy(ji,jj,jl)  =   zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp )  &
-               &             + zbt1 * psy(ji,jj,jl)  
-            psyy(ji,jj,jl) =   zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj,jl)                           &
-               &                      + 5.0 * ( zalf * zalf1 * ( psy(ji,jj,jl) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) & 
-               &             + zbt1 * psyy(ji,jj,jl)
-            psxy(ji,jj,jl) =   zbt  * (  zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj,jl)            &
-               &                      + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj,jl) ) )  &
-               &             + zbt1 * psxy(ji,jj,jl)
-            psx (ji,jj,jl) =   zbt * ( psx (ji,jj,jl) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj,jl)
-            psxx(ji,jj,jl) =   zbt * ( psxx(ji,jj,jl) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj,jl)
-         END_2D
-
-         DO_2D( 0, 0, 0, 0 )                   !  Flux from j+1 to j IF v LT 0.
-            zbt  =       zbet(ji,jj)
-            zbt1 = 1.0 - zbet(ji,jj)
-            psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) )
-            zalf          = zbt1 * zfm(ji,jj) / psm(ji,jj,jl)
-            zalf1         = 1.0 - zalf
-            ztemp         = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj)
-            !
-            ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * (  ps0(ji,jj,jl) + zf0(ji,jj) )
-            psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * (  zalf * zfy(ji,jj) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp )
-            psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj,jl) &
-               &                                            + 5.0 * ( zalf * zalf1 * ( - psy(ji,jj,jl) + zfy(ji,jj) )    &
-               &                                            + ( zalf1 - zalf ) * ztemp ) )
-            psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl)  &
-               &                                            + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj,jl) ) )
-            psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) + zfx (ji,jj) )
-            psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) + zfxx(ji,jj) )
-         END_2D
-
+         !
       END DO
-
-      !-- Lateral boundary conditions
-#if defined key_mpi3
-      CALL lbc_lnk_nc_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T',  1.0_wp, ps0 , 'T',  1.0_wp   &
-         &                                , psx             , 'T', -1.0_wp, psy , 'T', -1.0_wp   &   ! caution gradient ==> the sign changes
-         &                                , psxx            , 'T',  1.0_wp, psyy, 'T',  1.0_wp , psxy, 'T',  1.0_wp )
-#else
-      CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T',  1.0_wp, ps0 , 'T',  1.0_wp   &
-         &                                , psx             , 'T', -1.0_wp, psy , 'T', -1.0_wp   &   ! caution gradient ==> the sign changes
-         &                                , psxx            , 'T',  1.0_wp, psyy, 'T',  1.0_wp , psxy, 'T',  1.0_wp )
-#endif
       !
    END SUBROUTINE adv_y
@@ -897,30 +970,30 @@
             !
             !                                                        ! ice thickness
-            CALL iom_get( numrir, jpdom_auto, 'sxice' , sxice  )
-            CALL iom_get( numrir, jpdom_auto, 'syice' , syice  )
+            CALL iom_get( numrir, jpdom_auto, 'sxice' , sxice , psgn = -1._wp )
+            CALL iom_get( numrir, jpdom_auto, 'syice' , syice , psgn = -1._wp )
             CALL iom_get( numrir, jpdom_auto, 'sxxice', sxxice )
             CALL iom_get( numrir, jpdom_auto, 'syyice', syyice )
             CALL iom_get( numrir, jpdom_auto, 'sxyice', sxyice )
             !                                                        ! snow thickness
-            CALL iom_get( numrir, jpdom_auto, 'sxsn'  , sxsn   )
-            CALL iom_get( numrir, jpdom_auto, 'sysn'  , sysn   )
+            CALL iom_get( numrir, jpdom_auto, 'sxsn'  , sxsn  , psgn = -1._wp )
+            CALL iom_get( numrir, jpdom_auto, 'sysn'  , sysn  , psgn = -1._wp )
             CALL iom_get( numrir, jpdom_auto, 'sxxsn' , sxxsn  )
             CALL iom_get( numrir, jpdom_auto, 'syysn' , syysn  )
             CALL iom_get( numrir, jpdom_auto, 'sxysn' , sxysn  )
             !                                                        ! ice concentration
-            CALL iom_get( numrir, jpdom_auto, 'sxa'   , sxa    )
-            CALL iom_get( numrir, jpdom_auto, 'sya'   , sya    )
+            CALL iom_get( numrir, jpdom_auto, 'sxa'   , sxa   , psgn = -1._wp )
+            CALL iom_get( numrir, jpdom_auto, 'sya'   , sya   , psgn = -1._wp )
             CALL iom_get( numrir, jpdom_auto, 'sxxa'  , sxxa   )
             CALL iom_get( numrir, jpdom_auto, 'syya'  , syya   )
             CALL iom_get( numrir, jpdom_auto, 'sxya'  , sxya   )
             !                                                        ! ice salinity
-            CALL iom_get( numrir, jpdom_auto, 'sxsal' , sxsal  )
-            CALL iom_get( numrir, jpdom_auto, 'sysal' , sysal  )
+            CALL iom_get( numrir, jpdom_auto, 'sxsal' , sxsal , psgn = -1._wp )
+            CALL iom_get( numrir, jpdom_auto, 'sysal' , sysal , psgn = -1._wp )
             CALL iom_get( numrir, jpdom_auto, 'sxxsal', sxxsal )
             CALL iom_get( numrir, jpdom_auto, 'syysal', syysal )
             CALL iom_get( numrir, jpdom_auto, 'sxysal', sxysal )
             !                                                        ! ice age
-            CALL iom_get( numrir, jpdom_auto, 'sxage' , sxage  )
-            CALL iom_get( numrir, jpdom_auto, 'syage' , syage  )
+            CALL iom_get( numrir, jpdom_auto, 'sxage' , sxage , psgn = -1._wp )
+            CALL iom_get( numrir, jpdom_auto, 'syage' , syage , psgn = -1._wp )
             CALL iom_get( numrir, jpdom_auto, 'sxxage', sxxage )
             CALL iom_get( numrir, jpdom_auto, 'syyage', syyage )
@@ -929,6 +1002,6 @@
             DO jk = 1, nlay_s
                WRITE(zchar1,'(I2.2)') jk
-               znam = 'sxc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxc0 (:,:,jk,:) = z3d(:,:,:)
-               znam = 'syc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   syc0 (:,:,jk,:) = z3d(:,:,:)
+               znam = 'sxc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp )   ;   sxc0 (:,:,jk,:) = z3d(:,:,:)
+               znam = 'syc0'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp )   ;   syc0 (:,:,jk,:) = z3d(:,:,:)
                znam = 'sxxc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxxc0(:,:,jk,:) = z3d(:,:,:)
                znam = 'syyc0'//'_l'//zchar1 ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   syyc0(:,:,jk,:) = z3d(:,:,:)
@@ -938,6 +1011,6 @@
             DO jk = 1, nlay_i
                WRITE(zchar1,'(I2.2)') jk
-               znam = 'sxe'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxe (:,:,jk,:) = z3d(:,:,:)
-               znam = 'sye'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sye (:,:,jk,:) = z3d(:,:,:)
+               znam = 'sxe'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp )   ;   sxe (:,:,jk,:) = z3d(:,:,:)
+               znam = 'sye'//'_l'//zchar1   ;   CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp )   ;   sye (:,:,jk,:) = z3d(:,:,:)
                znam = 'sxxe'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxxe(:,:,jk,:) = z3d(:,:,:)
                znam = 'syye'//'_l'//zchar1  ;   CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   syye(:,:,jk,:) = z3d(:,:,:)
@@ -946,13 +1019,13 @@
             !
             IF( ln_pnd_LEV ) THEN                                    ! melt pond fraction
-               IF( iom_varid( numror, 'sxap', ldstop = .FALSE. ) > 0 ) THEN
-                  CALL iom_get( numrir, jpdom_auto, 'sxap' , sxap  )
-                  CALL iom_get( numrir, jpdom_auto, 'syap' , syap  )
+               IF( iom_varid( numrir, 'sxap', ldstop = .FALSE. ) > 0 ) THEN
+                  CALL iom_get( numrir, jpdom_auto, 'sxap' , sxap , psgn = -1._wp )
+                  CALL iom_get( numrir, jpdom_auto, 'syap' , syap , psgn = -1._wp )
                   CALL iom_get( numrir, jpdom_auto, 'sxxap', sxxap )
                   CALL iom_get( numrir, jpdom_auto, 'syyap', syyap )
                   CALL iom_get( numrir, jpdom_auto, 'sxyap', sxyap )
                   !                                                     ! melt pond volume
-                  CALL iom_get( numrir, jpdom_auto, 'sxvp' , sxvp  )
-                  CALL iom_get( numrir, jpdom_auto, 'syvp' , syvp  )
+                  CALL iom_get( numrir, jpdom_auto, 'sxvp' , sxvp , psgn = -1._wp )
+                  CALL iom_get( numrir, jpdom_auto, 'syvp' , syvp , psgn = -1._wp )
                   CALL iom_get( numrir, jpdom_auto, 'sxxvp', sxxvp )
                   CALL iom_get( numrir, jpdom_auto, 'syyvp', syyvp )
@@ -964,7 +1037,7 @@
                   !
                IF ( ln_pnd_lids ) THEN                               ! melt pond lid volume
-                  IF( iom_varid( numror, 'sxvl', ldstop = .FALSE. ) > 0 ) THEN
-                     CALL iom_get( numrir, jpdom_auto, 'sxvl' , sxvl  )
-                     CALL iom_get( numrir, jpdom_auto, 'syvl' , syvl  )
+                  IF( iom_varid( numrir, 'sxvl', ldstop = .FALSE. ) > 0 ) THEN
+                     CALL iom_get( numrir, jpdom_auto, 'sxvl' , sxvl , psgn = -1._wp )
+                     CALL iom_get( numrir, jpdom_auto, 'syvl' , syvl , psgn = -1._wp )
                      CALL iom_get( numrir, jpdom_auto, 'sxxvl', sxxvl )
                      CALL iom_get( numrir, jpdom_auto, 'syyvl', syyvl )
@@ -1081,4 +1154,55 @@
    END SUBROUTINE adv_pra_rst
 
+   SUBROUTINE icemax3D( pice , pmax )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE icemax3D ***                     
+      !! ** Purpose :  compute the max of the 9 points around
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:,:)      , INTENT(in ) ::   pice   ! input
+      REAL(wp), DIMENSION(:,:,:)      , INTENT(out) ::   pmax   ! output
+      REAL(wp), DIMENSION(2:jpim1,jpj)              ::   zmax   ! temporary array
+      INTEGER  ::   ji, jj, jl   ! dummy loop indices
+      !!----------------------------------------------------------------------
+      DO jl = 1, jpl
+         DO jj = Njs0-1, Nje0+1    
+            DO ji = Nis0, Nie0
+               zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jl), pice(ji-1,jj,jl), pice(ji+1,jj,jl) )
+            END DO
+         END DO
+         DO jj = Njs0, Nje0    
+            DO ji = Nis0, Nie0
+               pmax(ji,jj,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) )
+            END DO
+         END DO
+      END DO
+   END SUBROUTINE icemax3D
+
+   SUBROUTINE icemax4D( pice , pmax )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE icemax4D ***                     
+      !! ** Purpose :  compute the max of the 9 points around
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:,:,:)    , INTENT(in ) ::   pice   ! input
+      REAL(wp), DIMENSION(:,:,:,:)    , INTENT(out) ::   pmax   ! output
+      REAL(wp), DIMENSION(2:jpim1,jpj)              ::   zmax   ! temporary array
+      INTEGER  ::   jlay, ji, jj, jk, jl   ! dummy loop indices
+      !!----------------------------------------------------------------------
+      jlay = SIZE( pice , 3 )   ! size of input arrays
+      DO jl = 1, jpl
+         DO jk = 1, jlay
+            DO jj = Njs0-1, Nje0+1    
+               DO ji = Nis0, Nie0
+                  zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jk,jl), pice(ji-1,jj,jk,jl), pice(ji+1,jj,jk,jl) )
+               END DO
+            END DO
+            DO jj = Njs0, Nje0    
+               DO ji = Nis0, Nie0
+                  pmax(ji,jj,jk,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) )
+               END DO
+            END DO
+         END DO
+      END DO
+   END SUBROUTINE icemax4D
+
 #else
    !!----------------------------------------------------------------------
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icedyn_adv_umx.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icedyn_adv_umx.F90	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icedyn_adv_umx.F90	(revision 13877)
@@ -92,5 +92,5 @@
       INTEGER  ::   icycle                  ! number of sub-timestep for the advection
       REAL(wp) ::   zamsk                   ! 1 if advection of concentration, 0 if advection of other tracers
-      REAL(wp) ::   zdt, zvi_cen
+      REAL(wp) ::   zdt, z1_dt, zvi_cen
       REAL(wp), DIMENSION(1)                  ::   zcflprv, zcflnow   ! for global communication
       REAL(wp), DIMENSION(jpi,jpj)            ::   zudy, zvdx, zcu_box, zcv_box
@@ -104,4 +104,6 @@
       !
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zuv_ho, zvv_ho, zuv_ups, zvv_ups, z1_vi, z1_vs 
+      !! diagnostics
+      REAL(wp), DIMENSION(jpi,jpj)            ::   zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat      
       !!----------------------------------------------------------------------
       !
@@ -113,24 +115,8 @@
       ELSEWHERE                      ; zs_i(:,:,:) = 0._wp
       END WHERE
-      DO jl = 1, jpl
-         DO_2D( 0, 0, 0, 0 )
-            zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj  ,jl), ph_ip(ji  ,jj+1,jl), &
-               &                                               ph_ip(ji-1,jj  ,jl), ph_ip(ji  ,jj-1,jl), &
-               &                                               ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), &
-               &                                               ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) )
-            zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj  ,jl), ph_i (ji  ,jj+1,jl), &
-               &                                               ph_i (ji-1,jj  ,jl), ph_i (ji  ,jj-1,jl), &
-               &                                               ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), &
-               &                                               ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) )
-            zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj  ,jl), ph_s (ji  ,jj+1,jl), &
-               &                                               ph_s (ji-1,jj  ,jl), ph_s (ji  ,jj-1,jl), &
-               &                                               ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), &
-               &                                               ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) )
-            zsi_max (ji,jj,jl) = MAX( epsi20, zs_i (ji,jj,jl), zs_i (ji+1,jj  ,jl), zs_i (ji  ,jj+1,jl), &
-               &                                               zs_i (ji-1,jj  ,jl), zs_i (ji  ,jj-1,jl), &
-               &                                               zs_i (ji+1,jj+1,jl), zs_i (ji-1,jj-1,jl), &
-               &                                               zs_i (ji+1,jj-1,jl), zs_i (ji-1,jj+1,jl) )
-         END_2D
-      END DO
+      CALL icemax3D( ph_i , zhi_max )
+      CALL icemax3D( ph_s , zhs_max )
+      CALL icemax3D( ph_ip, zhip_max)
+      CALL icemax3D( zs_i , zsi_max )
 #if defined key_mpi3
       CALL lbc_lnk_nc_multi( 'icedyn_adv_umx', zhi_max, 'T', 1._wp, zhs_max, 'T', 1._wp, zhip_max, 'T', 1._wp, zsi_max, 'T', 1._wp )
@@ -149,28 +135,9 @@
          ELSEWHERE                      ; ze_s(:,:,jk,:) = 0._wp
          END WHERE
-      END DO
-      DO jl = 1, jpl
-         DO_3D( 0, 0, 0, 0, 1, nlay_i )
-            zei_max(ji,jj,jk,jl) = MAX( epsi20, ze_i(ji,jj,jk,jl), ze_i(ji+1,jj  ,jk,jl), ze_i(ji  ,jj+1,jk,jl), &
-               &                                                   ze_i(ji-1,jj  ,jk,jl), ze_i(ji  ,jj-1,jk,jl), &
-               &                                                   ze_i(ji+1,jj+1,jk,jl), ze_i(ji-1,jj-1,jk,jl), &
-               &                                                   ze_i(ji+1,jj-1,jk,jl), ze_i(ji-1,jj+1,jk,jl) )
-         END_3D
-      END DO
-      DO jl = 1, jpl
-         DO_3D( 0, 0, 0, 0, 1, nlay_s )
-            zes_max(ji,jj,jk,jl) = MAX( epsi20, ze_s(ji,jj,jk,jl), ze_s(ji+1,jj  ,jk,jl), ze_s(ji  ,jj+1,jk,jl), &
-               &                                                   ze_s(ji-1,jj  ,jk,jl), ze_s(ji  ,jj-1,jk,jl), &
-               &                                                   ze_s(ji+1,jj+1,jk,jl), ze_s(ji-1,jj-1,jk,jl), &
-               &                                                   ze_s(ji+1,jj-1,jk,jl), ze_s(ji-1,jj+1,jk,jl) )
-         END_3D
-      END DO
-#if defined key_mpi3
-      CALL lbc_lnk_nc_multi( 'icedyn_adv_pra', zei_max, 'T', 1. )
-      CALL lbc_lnk_nc_multi( 'icedyn_adv_pra', zes_max, 'T', 1. )
-#else
-      CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1. )
-      CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1. )
-#endif
+      END DO   
+      CALL icemax4D( ze_i , zei_max )
+      CALL icemax4D( ze_s , zes_max )
+      CALL lbc_lnk( 'icedyn_adv_umx', zei_max, 'T', 1._wp )
+      CALL lbc_lnk( 'icedyn_adv_umx', zes_max, 'T', 1._wp )
       !
       !
@@ -188,4 +155,5 @@
       ENDIF
       zdt = rDt_ice / REAL(icycle)
+      z1_dt = 1._wp / zdt
 
       ! --- transport --- !
@@ -216,4 +184,10 @@
       !---------------!
       DO jt = 1, icycle
+
+         ! diagnostics
+         zdiag_adv_mass(:,:) =   SUM(  pv_i(:,:,:) , dim=3 ) * rhoi + SUM(  pv_s(:,:,:) , dim=3 ) * rhos
+         zdiag_adv_salt(:,:) =   SUM( psv_i(:,:,:) , dim=3 ) * rhoi
+         zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) &
+            &                  - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 )
 
          ! record at_i before advection (for open water)
@@ -386,4 +360,17 @@
             ENDIF
          ENDIF
+
+         ! --- Lateral boundary conditions --- !
+         IF    ( ln_pnd_LEV .AND. ln_pnd_lids ) THEN
+            CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp &
+               &                                , pa_ip,'T',1._wp, pv_ip,'T',1._wp, pv_il,'T',1._wp )
+         ELSEIF( ln_pnd_LEV .AND. .NOT.ln_pnd_lids ) THEN
+            CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp &
+               &                                , pa_ip,'T',1._wp, pv_ip,'T',1._wp )
+         ELSE
+            CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp )
+         ENDIF
+         CALL lbc_lnk( 'icedyn_adv_umx', pe_i, 'T', 1._wp )
+         CALL lbc_lnk( 'icedyn_adv_umx', pe_s, 'T', 1._wp )
          !
          !== Open water area ==!
@@ -399,4 +386,12 @@
 #endif
          !
+         ! --- diagnostics --- !
+         diag_adv_mass(:,:) = diag_adv_mass(:,:) + (   SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos &
+            &                                        - zdiag_adv_mass(:,:) ) * z1_dt
+         diag_adv_salt(:,:) = diag_adv_salt(:,:) + (   SUM( psv_i(:,:,:) , dim=3 ) * rhoi &
+            &                                        - zdiag_adv_salt(:,:) ) * z1_dt
+         diag_adv_heat(:,:) = diag_adv_heat(:,:) + ( - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) &
+            &                                        - SUM(SUM( pe_s(:,:,1:nlay_s,:) , dim=4 ), dim=3 ) &
+            &                                        - zdiag_adv_heat(:,:) ) * z1_dt
          !
          ! --- Ensure non-negative fields and in-bound thicknesses --- !
@@ -458,6 +453,5 @@
       !!             work on H (and not V). It is partly related to the multi-category approach
       !!             Therefore, after advection we limit the thickness to the largest value of the 9-points around (only if ice
-      !!             concentration is small). Since we do not limit S and T, large values can occur at the edge but it does not really matter
-      !!             since sv_i and e_i are still good.
+      !!             concentration is small). We also limit S and T.
       !!----------------------------------------------------------------------
       REAL(wp)                        , INTENT(in   )           ::   pamsk            ! advection of concentration (1) or other tracers (0)
@@ -503,5 +497,5 @@
       IF( pamsk == 0._wp ) THEN
          DO jl = 1, jpl
-            DO_2D( 1, 0, 1, 0 )
+            DO_2D( 0, 0, 1, 0 )
                IF( ABS( pu(ji,jj) ) > epsi10 ) THEN
                   zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc    (ji,jj,jl) / pu(ji,jj)
@@ -512,4 +506,6 @@
                ENDIF
                !
+            END_2D
+            DO_2D( 1, 0, 0, 0 )
                IF( ABS( pv(ji,jj) ) > epsi10 ) THEN
                   zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc    (ji,jj,jl) / pv(ji,jj)
@@ -550,7 +546,11 @@
       IF( PRESENT( pua_ho ) ) THEN
          DO jl = 1, jpl
-            DO_2D( 1, 0, 1, 0 )
-               pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) ; pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl)
-               pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) ; pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl)
+            DO_2D( 0, 0, 1, 0 )
+               pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl)
+               pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl)
+            END_2D
+            DO_2D( 1, 0, 0, 0 )
+               pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl)
+               pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl)
             END_2D
          END DO
@@ -566,9 +566,4 @@
          END_2D
       END DO
-#if defined key_mpi3
-      CALL lbc_lnk_nc_multi( 'icedyn_adv_umx', ptc, 'T',  1.0_wp )
-#else
-      CALL lbc_lnk( 'icedyn_adv_umx', ptc, 'T',  1.0_wp )
-#endif
       !
    END SUBROUTINE adv_umx
@@ -609,5 +604,5 @@
             !
             DO jl = 1, jpl              !-- flux in x-direction
-               DO_2D( 1, 0, 1, 0 )
+               DO_2D( 1, 1, 1, 0 )
                   pfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj,jl)
                END_2D
@@ -615,5 +610,5 @@
             !
             DO jl = 1, jpl              !-- first guess of tracer from u-flux
-               DO_2D( 0, 0, 0, 0 )
+               DO_2D( 1, 1, 0, 0 )
                   ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj,jl) )              &
                      &   + ( pu     (ji,jj   ) - pu     (ji-1,jj   ) ) * pt(ji,jj,jl) * (1.-pamsk)
@@ -622,12 +617,7 @@
                END_2D
             END DO
-#if defined key_mpi3
-            CALL lbc_lnk_nc_multi( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )
-#else
-            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )
-#endif
             !
             DO jl = 1, jpl              !-- flux in y-direction
-               DO_2D( 1, 0, 1, 0 )
+               DO_2D( 1, 0, 0, 0 )
                   pfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * zpt(ji,jj+1,jl)
                END_2D
@@ -637,5 +627,5 @@
             !
             DO jl = 1, jpl              !-- flux in y-direction
-               DO_2D( 1, 0, 1, 0 )
+               DO_2D( 1, 0, 1, 1 )
                   pfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1,jl)
                END_2D
@@ -643,5 +633,5 @@
             !
             DO jl = 1, jpl              !-- first guess of tracer from v-flux
-               DO_2D( 0, 0, 0, 0 )
+               DO_2D( 0, 0, 1, 1 )
                   ztra = - ( pfv_ups(ji,jj,jl) - pfv_ups(ji,jj-1,jl) )  &
                      &   + ( pv     (ji,jj   ) - pv     (ji,jj-1   ) ) * pt(ji,jj,jl) * (1.-pamsk)
@@ -650,12 +640,7 @@
                END_2D
             END DO
-#if defined key_mpi3
-            CALL lbc_lnk_nc_multi( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )
-#else
-            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )
-#endif
             !
             DO jl = 1, jpl              !-- flux in x-direction
-               DO_2D( 1, 0, 1, 0 )
+               DO_2D( 0, 0, 1, 0 )
                   pfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * zpt(ji+1,jj,jl)
                END_2D
@@ -710,6 +695,8 @@
          !
          DO jl = 1, jpl
-            DO_2D( 1, 0, 1, 0 )
+            DO_2D( 1, 1, 1, 0 )
                pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj  ,jl) )
+            END_2D
+            DO_2D( 1, 0, 1, 1 )
                pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji  ,jj+1,jl) )
             END_2D
@@ -728,5 +715,5 @@
             !
             DO jl = 1, jpl              !-- flux in x-direction
-               DO_2D( 1, 0, 1, 0 )
+               DO_2D( 1, 1, 1, 0 )
                   pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) )
                END_2D
@@ -735,5 +722,5 @@
 
             DO jl = 1, jpl              !-- first guess of tracer from u-flux
-               DO_2D( 0, 0, 0, 0 )
+               DO_2D( 1, 1, 0, 0 )
                   ztra = - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) )              &
                      &   + ( pu    (ji,jj   ) - pu    (ji-1,jj   ) ) * pt(ji,jj,jl) * (1.-pamsk)
@@ -742,12 +729,7 @@
                END_2D
             END DO
-#if defined key_mpi3
-            CALL lbc_lnk_nc_multi( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )
-#else
-            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )
-#endif
 
             DO jl = 1, jpl              !-- flux in y-direction
-               DO_2D( 1, 0, 1, 0 )
+               DO_2D( 1, 0, 0, 0 )
                   pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji,jj+1,jl) )
                END_2D
@@ -758,5 +740,5 @@
             !
             DO jl = 1, jpl              !-- flux in y-direction
-               DO_2D( 1, 0, 1, 0 )
+               DO_2D( 1, 0, 1, 1 )
                   pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) )
                END_2D
@@ -765,5 +747,5 @@
             !
             DO jl = 1, jpl              !-- first guess of tracer from v-flux
-               DO_2D( 0, 0, 0, 0 )
+               DO_2D( 0, 0, 1, 1 )
                   ztra = - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) )  &
                      &   + ( pv    (ji,jj   ) - pv    (ji,jj-1   ) ) * pt(ji,jj,jl) * (1.-pamsk)
@@ -772,12 +754,7 @@
                END_2D
             END DO
-#if defined key_mpi3
-            CALL lbc_lnk_nc_multi( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )
-#else
-            CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )
-#endif
             !
             DO jl = 1, jpl              !-- flux in x-direction
-               DO_2D( 1, 0, 1, 0 )
+               DO_2D( 0, 0, 1, 0 )
                   pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji+1,jj,jl) )
                END_2D
@@ -952,5 +929,5 @@
          !        
          DO jl = 1, jpl
-            DO_2D( 1, 0, 1, 0 )
+            DO_2D( 0, 0, 1, 0 )
                pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   &
                   &                                         - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) )
@@ -961,5 +938,5 @@
          !
          DO jl = 1, jpl
-            DO_2D( 1, 0, 1, 0 )
+            DO_2D( 0, 0, 1, 0 )
                zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
                pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   &
@@ -971,5 +948,5 @@
          !
          DO jl = 1, jpl
-            DO_2D( 1, 0, 1, 0 )
+            DO_2D( 0, 0, 1, 0 )
                zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
                zdx2 = e1u(ji,jj) * e1u(ji,jj)
@@ -985,5 +962,5 @@
          !
          DO jl = 1, jpl
-            DO_2D( 1, 0, 1, 0 )
+            DO_2D( 0, 0, 1, 0 )
                zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
                zdx2 = e1u(ji,jj) * e1u(ji,jj)
@@ -999,5 +976,5 @@
          !
          DO jl = 1, jpl
-            DO_2D( 1, 0, 1, 0 )
+            DO_2D( 0, 0, 1, 0 )
                zcu  = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj)
                zdx2 = e1u(ji,jj) * e1u(ji,jj)
@@ -1020,5 +997,5 @@
       IF( ll_neg ) THEN
          DO jl = 1, jpl
-            DO_2D( 1, 0, 1, 0 )
+            DO_2D( 0, 0, 1, 0 )
                IF( pt_u(ji,jj,jl) < 0._wp .OR. ( imsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN
                   pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (                                pt(ji+1,jj,jl) + pt(ji,jj,jl)   &
@@ -1030,5 +1007,5 @@
       !                                                     !-- High order flux in i-direction  --!
       DO jl = 1, jpl
-         DO_2D( 1, 0, 1, 0 )
+         DO_2D( 0, 0, 1, 0 )
             pfu_ho(ji,jj,jl) = pu(ji,jj) * pt_u(ji,jj,jl)
          END_2D
@@ -1096,5 +1073,5 @@
       CASE( 1 )                                                !==  1st order central TIM  ==! (Eq. 21)
          DO jl = 1, jpl
-            DO_2D( 1, 0, 1, 0 )
+            DO_2D( 1, 0, 0, 0 )
                pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                                pt(ji,jj+1,jl) + pt(ji,jj,jl)   &
                   &                                         - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) )
@@ -1104,5 +1081,5 @@
       CASE( 2 )                                                !==  2nd order central TIM  ==! (Eq. 23)
          DO jl = 1, jpl
-            DO_2D( 1, 0, 1, 0 )
+            DO_2D( 1, 0, 0, 0 )
                zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
                pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                                pt(ji,jj+1,jl) + pt(ji,jj,jl)   &
@@ -1113,5 +1090,5 @@
       CASE( 3 )                                                !==  3rd order central TIM  ==! (Eq. 24)
          DO jl = 1, jpl
-            DO_2D( 1, 0, 1, 0 )
+            DO_2D( 1, 0, 0, 0 )
                zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
                zdy2 = e2v(ji,jj) * e2v(ji,jj)
@@ -1126,5 +1103,5 @@
       CASE( 4 )                                                !==  4th order central TIM  ==! (Eq. 27)
          DO jl = 1, jpl
-            DO_2D( 1, 0, 1, 0 )
+            DO_2D( 1, 0, 0, 0 )
                zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
                zdy2 = e2v(ji,jj) * e2v(ji,jj)
@@ -1139,5 +1116,5 @@
       CASE( 5 )                                                !==  5th order central TIM  ==! (Eq. 29)
          DO jl = 1, jpl
-            DO_2D( 1, 0, 1, 0 )
+            DO_2D( 1, 0, 0, 0 )
                zcv  = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj)
                zdy2 = e2v(ji,jj) * e2v(ji,jj)
@@ -1160,5 +1137,5 @@
       IF( ll_neg ) THEN
          DO jl = 1, jpl
-            DO_2D( 1, 0, 1, 0 )
+            DO_2D( 1, 0, 0, 0 )
                IF( pt_v(ji,jj,jl) < 0._wp .OR. ( jmsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN
                   pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt(ji,jj+1,jl) + pt(ji,jj,jl) )  &
@@ -1170,5 +1147,5 @@
       !                                                     !-- High order flux in j-direction  --!
       DO jl = 1, jpl
-         DO_2D( 1, 0, 1, 0 )
+         DO_2D( 1, 0, 0, 0 )
             pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl)
          END_2D
@@ -1206,6 +1183,8 @@
       ! --------------------------------------------------
       DO jl = 1, jpl
-         DO_2D( 1, 0, 1, 0 )
+         DO_2D( 0, 0, 1, 0 )
             pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl)
+         END_2D
+         DO_2D( 1, 0, 0, 0 )
             pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl)
          END_2D
@@ -1325,5 +1304,5 @@
       ! ---------------------------------
       DO jl = 1, jpl
-         DO_2D( 1, 0, 1, 0 )
+         DO_2D( 0, 0, 1, 0 )
             zau = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji+1,jj,jl) )
             zbu = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji+1,jj,jl) )
@@ -1336,5 +1315,5 @@
          END_2D
 
-         DO_2D( 1, 0, 1, 0 )
+         DO_2D( 1, 0, 0, 0 )
             zav = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji,jj+1,jl) )
             zbv = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji,jj+1,jl) )
@@ -1709,4 +1688,54 @@
    END SUBROUTINE Hsnow
 
+   SUBROUTINE icemax3D( pice , pmax )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE icemax3D ***                     
+      !! ** Purpose :  compute the max of the 9 points around
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:,:)      , INTENT(in ) ::   pice   ! input
+      REAL(wp), DIMENSION(:,:,:)      , INTENT(out) ::   pmax   ! output
+      REAL(wp), DIMENSION(2:jpim1,jpj)              ::   zmax   ! temporary array
+      INTEGER  ::   ji, jj, jl   ! dummy loop indices
+      !!----------------------------------------------------------------------
+      DO jl = 1, jpl
+         DO jj = Njs0-1, Nje0+1    
+            DO ji = Nis0, Nie0
+               zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jl), pice(ji-1,jj,jl), pice(ji+1,jj,jl) )
+            END DO
+         END DO
+         DO jj = Njs0, Nje0    
+            DO ji = Nis0, Nie0
+               pmax(ji,jj,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) )
+            END DO
+         END DO
+      END DO
+   END SUBROUTINE icemax3D
+
+   SUBROUTINE icemax4D( pice , pmax )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE icemax4D ***                     
+      !! ** Purpose :  compute the max of the 9 points around
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:,:,:)    , INTENT(in ) ::   pice   ! input
+      REAL(wp), DIMENSION(:,:,:,:)    , INTENT(out) ::   pmax   ! output
+      REAL(wp), DIMENSION(2:jpim1,jpj)              ::   zmax   ! temporary array
+      INTEGER  ::   jlay, ji, jj, jk, jl   ! dummy loop indices
+      !!----------------------------------------------------------------------
+      jlay = SIZE( pice , 3 )   ! size of input arrays
+      DO jl = 1, jpl
+         DO jk = 1, jlay
+            DO jj = Njs0-1, Nje0+1    
+               DO ji = Nis0, Nie0
+                  zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jk,jl), pice(ji-1,jj,jk,jl), pice(ji+1,jj,jk,jl) )
+               END DO
+            END DO
+            DO jj = Njs0, Nje0    
+               DO ji = Nis0, Nie0
+                  pmax(ji,jj,jk,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) )
+               END DO
+            END DO
+         END DO
+      END DO
+   END SUBROUTINE icemax4D
 
 #else
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icedyn_rdgrft.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icedyn_rdgrft.F90	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icedyn_rdgrft.F90	(revision 13877)
@@ -349,5 +349,5 @@
                ELSEIF( zGsum(ji,jl-1) < rn_gstar ) THEN
                   apartf(ji,jl) = z1_gstar * ( rn_gstar     - zGsum(ji,jl-1) ) *  &
-                     &                       ( 2._wp - ( zGsum(ji,jl-1) + rn_gstar        ) * z1_gstar )
+                     &                       ( 2._wp - ( zGsum(ji,jl-1) + rn_gstar     ) * z1_gstar )
                ELSE
                   apartf(ji,jl) = 0._wp
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icedyn_rhg_evp.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icedyn_rhg_evp.F90	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icedyn_rhg_evp.F90	(revision 13877)
@@ -148,4 +148,5 @@
       !
       REAL(wp), DIMENSION(jpi,jpj) ::   zds                             ! shear
+      REAL(wp), DIMENSION(jpi,jpj) ::   zten_i                          ! tension
       REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12                  ! stress tensor components
       REAL(wp), DIMENSION(jpi,jpj) ::   zsshdyn                         ! array used for the calculation of ice surface slope:
@@ -170,5 +171,6 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice
       !! --- diags
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig1, zsig2, zsig3
+      REAL(wp) ::   zsig1, zsig2, zsig12, zfac, z1_strength
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig_I, zsig_II, zsig1_p, zsig2_p         
       !! --- SIMIP diags
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s)
@@ -760,4 +762,6 @@
             &   ) * r1_e1e2t(ji,jj)
          zdt2 = zdt * zdt
+
+         zten_i(ji,jj) = zdt
          
          ! shear**2 at T points (doc eq. A16)
@@ -775,21 +779,18 @@
          
          ! delta at T points
-         zdelta(ji,jj)   = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )  
-         rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta(ji,jj) ) ) ! 0 if delta=0
-         pdelta_i(ji,jj) = zdelta(ji,jj) + rn_creepl * rswitch
+         zfac            = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) ! delta  
+         rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zfac ) ) ! 0 if delta=0
+         pdelta_i(ji,jj) = zfac + rn_creepl * rswitch ! delta+creepl
 
       END_2D
 #if defined key_mpi3
-      CALL lbc_lnk_nc_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1.0_wp, pdivu_i, 'T', 1.0_wp, pdelta_i, 'T', 1.0_wp )
-#else
-      CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1.0_wp, pdivu_i, 'T', 1.0_wp, pdelta_i, 'T', 1.0_wp )
+      CALL lbc_lnk_nc_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1._wp, pdivu_i, 'T', 1._wp, pdelta_i, 'T', 1._wp, zten_i, 'T', 1._wp, &
+         &                                  zs1     , 'T', 1._wp, zs2    , 'T', 1._wp, zs12    , 'F', 1._wp )
+#else
+      CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1._wp, pdivu_i, 'T', 1._wp, pdelta_i, 'T', 1._wp, zten_i, 'T', 1._wp, &
+         &                                  zs1     , 'T', 1._wp, zs2    , 'T', 1._wp, zs12    , 'F', 1._wp )
 #endif
       
       ! --- Store the stress tensor for the next time step --- !
-#if defined key_mpi3
-      CALL lbc_lnk_nc_multi( 'icedyn_rhg_evp', zs1, 'T', 1.0_wp, zs2, 'T', 1.0_wp, zs12, 'F', 1.0_wp )
-#else
-      CALL lbc_lnk_multi( 'icedyn_rhg_evp', zs1, 'T', 1.0_wp, zs2, 'T', 1.0_wp, zs12, 'F', 1.0_wp )
-#endif
       pstress1_i (:,:) = zs1 (:,:)
       pstress2_i (:,:) = zs2 (:,:)
@@ -825,41 +826,67 @@
       IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength
 
-      ! --- stress tensor --- !
-      IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') .OR. iom_use('normstr') .OR. iom_use('sheastr') ) THEN
-         !
-         ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) )
+      ! --- Stress tensor invariants (SIMIP diags) --- !
+      IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN
+         !
+         ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
          !         
-         DO_2D( 0, 0, 0, 0 )
-            zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji  ,jj-1) * pstress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point
-               &      zmsk00(ji  ,jj) * pstress12_i(ji  ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) )  &
-               &    / MAX( 1._wp, zmsk00(ji-1,jj) + zmsk00(ji,jj-1) + zmsk00(ji,jj) + zmsk00(ji-1,jj-1) )
-
-            zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress  
-
-            zdum2 = zmsk00(ji,jj) / MAX( 1._wp, strength(ji,jj) )
-
-!!               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002)
-!!               zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002)
-!!               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress
-!!                                                                                                               ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11))
-            zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015
-            zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear )                     ! shear stress
-            zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 )
-         END_2D
-#if defined key_mpi3
-         CALL lbc_lnk_nc_multi( 'icedyn_rhg_evp', zsig1, 'T', 1.0_wp, zsig2, 'T', 1.0_wp, zsig3, 'T', 1.0_wp )
-#else
-         CALL lbc_lnk_multi( 'icedyn_rhg_evp', zsig1, 'T', 1.0_wp, zsig2, 'T', 1.0_wp, zsig3, 'T', 1.0_wp )
-#endif
-         !
-         CALL iom_put( 'isig1' , zsig1 )
-         CALL iom_put( 'isig2' , zsig2 )
-         CALL iom_put( 'isig3' , zsig3 )
-         !
-         ! Stress tensor invariants (normal and shear stress N/m)
-         IF( iom_use('normstr') )   CALL iom_put( 'normstr' ,       ( zs1(:,:) + zs2(:,:) )                       * zmsk00(:,:) ) ! Normal stress
-         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress
-
-         DEALLOCATE( zsig1 , zsig2 , zsig3 )
+         DO_2D( 1, 1, 1, 1 )
+            
+            ! Ice stresses
+            ! sigma1, sigma2, sigma12 are some useful recombination of the stresses (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013)
+            ! These are NOT stress tensor components, neither stress invariants, neither stress principal components
+            ! I know, this can be confusing...
+            zfac             =   strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl ) 
+            zsig1            =   zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) )
+            zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj)
+            zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj)
+            
+            ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008)
+            zsig_I (ji,jj)   =   zsig1 * 0.5_wp                                           ! 1st stress invariant, aka average normal stress, aka negative pressure
+            zsig_II(ji,jj)   =   SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) )  ! 2nd  ''       '', aka maximum shear stress
+               
+         END_2D         
+         !
+         ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags - definitions following Coon (1974) and Feltham (2008)
+         IF( iom_use('normstr') )   CALL iom_put( 'normstr', zsig_I (:,:) * zmsk00(:,:) ) ! Normal stress
+         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr', zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress
+         
+         DEALLOCATE ( zsig_I, zsig_II )
+         
+      ENDIF
+
+      ! --- Normalized stress tensor principal components --- !
+      ! This are used to plot the normalized yield curve, see Lemieux & Dupont, 2020
+      ! Recommendation 1 : we use ice strength, not replacement pressure
+      ! Recommendation 2 : need to use deformations at PREVIOUS iterate for viscosities
+      IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN
+         !
+         ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )         
+         !         
+         DO_2D( 1, 1, 1, 1 )
+            
+            ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates 
+            !                        and **deformations** at current iterates
+            !                        following Lemieux & Dupont (2020)
+            zfac             =   zp_delt(ji,jj)
+            zsig1            =   zfac * ( pdivu_i(ji,jj) - ( zdelta(ji,jj) + rn_creepl ) )
+            zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj)
+            zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj)
+            
+            ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point
+            zsig_I(ji,jj)    =   zsig1 * 0.5_wp                                            ! 1st stress invariant, aka average normal stress, aka negative pressure
+            zsig_II(ji,jj)   =   SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) )   ! 2nd  ''       '', aka maximum shear stress
+            
+            ! Normalized  principal stresses (used to display the ellipse)
+            z1_strength      =   1._wp / MAX( 1._wp, strength(ji,jj) )
+            zsig1_p(ji,jj)   =   ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength
+            zsig2_p(ji,jj)   =   ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength
+         END_2D              
+         !
+         CALL iom_put( 'sig1_pnorm' , zsig1_p ) 
+         CALL iom_put( 'sig2_pnorm' , zsig2_p ) 
+
+         DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II )
+         
       ENDIF
 
@@ -1008,5 +1035,5 @@
          istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) )
          ! close file
-         IF( kt == nitend )   istatus = NF90_CLOSE(ncvgid)
+         IF( kt == nitend - nn_fsbc + 1 )   istatus = NF90_CLOSE(ncvgid)
       ENDIF
       
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/iceitd.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/iceitd.F90	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/iceitd.F90	(revision 13877)
@@ -627,23 +627,23 @@
          END_2D
          !
-!!clem   CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,jl) )
-         CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl) )
-         CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl) )
-         !
-         DO ji = 1, npti
-            jdonor(ji,jl)  = jl 
-            ! how much of a_i you send in cat sup is somewhat arbitrary
-!!clem: these do not work properly after a restart (I do not know why) => not sure it is still true
-!!          zdaice(ji,jl)  = a_i_1d(ji) * ( h_i_1d(ji) - hi_max(jl) + epsi10 ) / h_i_1d(ji)  
-!!          zdvice(ji,jl)  = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 )
-!!clem: these do not work properly after a restart (I do not know why) => not sure it is still true
-!!          zdaice(ji,jl)  = a_i_1d(ji)
-!!          zdvice(ji,jl)  = v_i_1d(ji)
-!!clem: these are from UCL and work ok
-            zdaice(ji,jl)  = a_i_1d(ji) * 0.5_wp
-            zdvice(ji,jl)  = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp
-         END DO
-         !
-         IF( npti > 0 ) THEN
+         IF( npti > 0 ) THEN            
+            !!clem   CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,jl) )
+            CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl) )
+            CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl) )
+            !
+            DO ji = 1, npti
+               jdonor(ji,jl)  = jl 
+               ! how much of a_i you send in cat sup is somewhat arbitrary
+               !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true
+               !!          zdaice(ji,jl)  = a_i_1d(ji) * ( h_i_1d(ji) - hi_max(jl) + epsi10 ) / h_i_1d(ji)  
+               !!          zdvice(ji,jl)  = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 )
+               !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true
+               !!          zdaice(ji,jl)  = a_i_1d(ji)
+               !!          zdvice(ji,jl)  = v_i_1d(ji)
+               !!clem: these are from UCL and work ok
+               zdaice(ji,jl)  = a_i_1d(ji) * 0.5_wp
+               zdvice(ji,jl)  = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp
+            END DO
+            !
             CALL itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )  ! Shift jl=>jl+1
             ! Reset shift parameters
@@ -666,13 +666,13 @@
          END_2D
          !
-         CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl+1) ) ! jl+1 is ok
-         CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl+1) ) ! jl+1 is ok
-         DO ji = 1, npti
-            jdonor(ji,jl) = jl + 1
-            zdaice(ji,jl) = a_i_1d(ji) 
-            zdvice(ji,jl) = v_i_1d(ji)
-         END DO
-         !
          IF( npti > 0 ) THEN
+            CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl+1) ) ! jl+1 is ok
+            CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl+1) ) ! jl+1 is ok
+            DO ji = 1, npti
+               jdonor(ji,jl) = jl + 1
+               zdaice(ji,jl) = a_i_1d(ji) 
+               zdvice(ji,jl) = v_i_1d(ji)
+            END DO
+            !
             CALL itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )  ! Shift jl+1=>jl
             ! Reset shift parameters
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icestp.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icestp.F90	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icestp.F90	(revision 13877)
@@ -55,5 +55,4 @@
    USE icedyn         ! sea-ice: dynamics
    USE icethd         ! sea-ice: thermodynamics
-   USE icecor         ! sea-ice: corrections
    USE iceupdate      ! sea-ice: sea surface boundary condition update
    USE icedia         ! sea-ice: budget diagnostics
@@ -86,4 +85,6 @@
    PUBLIC   ice_init   ! called by sbcmod.F90
 
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/ICE 4.0 , NEMO Consortium (2018)
@@ -160,4 +161,6 @@
          IF( ln_icedyn .AND. .NOT.lk_c1d )   &
             &                           CALL ice_dyn( kt, Kmm )       ! -- Ice dynamics
+         !
+                                        CALL diag_trends( 1 )         ! record dyn trends
          !
          !                          !==  lateral boundary conditions  ==!
@@ -188,6 +191,5 @@
          IF( ln_icethd )                CALL ice_thd( kt )            ! -- Ice thermodynamics      
          !
-                                        CALL ice_cor( kt , 2 )        ! -- Corrections
-         !
+                                        CALL diag_trends( 2 )         ! record thermo trends
                                         CALL ice_var_glo2eqv          ! necessary calls (at least for coupling)
                                         CALL ice_var_agg( 2 )         ! necessary calls (at least for coupling)
@@ -196,4 +198,6 @@
          !
          IF( ln_icediahsb )             CALL ice_dia( kt )            ! -- Diagnostics outputs 
+         !
+         IF( ln_icediachk )             CALL ice_drift_wri( kt )      ! -- Diagnostics outputs for conservation 
          !
                                         CALL ice_wri( kt )            ! -- Ice outputs 
@@ -208,5 +212,5 @@
       ! --- Ocean time step --- !
       !-------------------------!
-      IF( ln_icedyn )                   CALL ice_update_tau( kt, uu(:,:,1,Kbb), vv(:,:,1,Kbb) )   ! -- update surface ocean stresses
+      CALL ice_update_tau( kt, uu(:,:,1,Kbb), vv(:,:,1,Kbb) )         ! -- update surface ocean stresses
 !!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!!
       !
@@ -281,4 +285,6 @@
       !
       CALL ice_dia_init                ! initialization for diags
+      !
+      CALL ice_drift_init              ! initialization for diags of conservation
       !
       fr_i  (:,:)   = at_i(:,:)        ! initialisation of sea-ice fraction
@@ -341,6 +347,4 @@
       ENDIF
       !
-      IF( ln_bdy .AND. ln_icediachk )   CALL ctl_warn('par_init: online conservation check does not work with BDY')
-      !
       rDt_ice   = REAL(nn_fsbc) * rn_Dt          !--- sea-ice timestep and its inverse
       r1_Dt_ice = 1._wp / rDt_ice
@@ -392,52 +396,107 @@
       !!               of the time step
       !!----------------------------------------------------------------------
-      INTEGER  ::   ji, jj      ! dummy loop index
-      !!----------------------------------------------------------------------
-      sfx    (:,:) = 0._wp   ;
-      sfx_bri(:,:) = 0._wp   ;   sfx_lam(:,:) = 0._wp
-      sfx_sni(:,:) = 0._wp   ;   sfx_opw(:,:) = 0._wp
-      sfx_bog(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp
-      sfx_bom(:,:) = 0._wp   ;   sfx_sum(:,:) = 0._wp
-      sfx_res(:,:) = 0._wp   ;   sfx_sub(:,:) = 0._wp
-      !
-      wfx_snw(:,:) = 0._wp   ;   wfx_ice(:,:) = 0._wp
-      wfx_sni(:,:) = 0._wp   ;   wfx_opw(:,:) = 0._wp
-      wfx_bog(:,:) = 0._wp   ;   wfx_dyn(:,:) = 0._wp
-      wfx_bom(:,:) = 0._wp   ;   wfx_sum(:,:) = 0._wp
-      wfx_res(:,:) = 0._wp   ;   wfx_sub(:,:) = 0._wp
-      wfx_spr(:,:) = 0._wp   ;   wfx_lam(:,:) = 0._wp  
-      wfx_snw_dyn(:,:) = 0._wp ; wfx_snw_sum(:,:) = 0._wp
-      wfx_snw_sub(:,:) = 0._wp ; wfx_ice_sub(:,:) = 0._wp
-      wfx_snw_sni(:,:) = 0._wp 
-      wfx_pnd(:,:) = 0._wp
-
-      hfx_thd(:,:) = 0._wp   ;
-      hfx_snw(:,:) = 0._wp   ;   hfx_opw(:,:) = 0._wp
-      hfx_bog(:,:) = 0._wp   ;   hfx_dyn(:,:) = 0._wp
-      hfx_bom(:,:) = 0._wp   ;   hfx_sum(:,:) = 0._wp
-      hfx_res(:,:) = 0._wp   ;   hfx_sub(:,:) = 0._wp
-      hfx_spr(:,:) = 0._wp   ;   hfx_dif(:,:) = 0._wp
-      hfx_err_dif(:,:) = 0._wp
-      wfx_err_sub(:,:) = 0._wp
-      !
-      diag_heat(:,:) = 0._wp ;   diag_sice(:,:) = 0._wp
-      diag_vice(:,:) = 0._wp ;   diag_vsnw(:,:) = 0._wp
-
-      ! SIMIP diagnostics
-      qcn_ice_bot(:,:,:) = 0._wp ; qcn_ice_top(:,:,:) = 0._wp ! conductive fluxes
-      t_si       (:,:,:) = rt0   ! temp at the ice-snow interface
-
-      tau_icebfr (:,:)   = 0._wp   ! landfast ice param only (clem: important to keep the init here)
-      cnd_ice    (:,:,:) = 0._wp   ! initialisation: effective conductivity at the top of ice/snow (ln_cndflx=T)
-      qcn_ice    (:,:,:) = 0._wp   ! initialisation: conductive flux (ln_cndflx=T & ln_cndemule=T)
-      qtr_ice_bot(:,:,:) = 0._wp   ! initialization: part of solar radiation transmitted through the ice needed at least for outputs
-      qsb_ice_bot(:,:)   = 0._wp   ! (needed if ln_icethd=F)
-      !
-      ! for control checks (ln_icediachk)
-      diag_trp_vi(:,:) = 0._wp   ;   diag_trp_vs(:,:) = 0._wp
-      diag_trp_ei(:,:) = 0._wp   ;   diag_trp_es(:,:) = 0._wp
-      diag_trp_sv(:,:) = 0._wp
-      
+      INTEGER  ::   ji, jj, jl      ! dummy loop index
+      !!----------------------------------------------------------------------
+
+      DO_2D( 1, 1, 1, 1 )
+         sfx    (ji,jj) = 0._wp   ;
+         sfx_bri(ji,jj) = 0._wp   ;   sfx_lam(ji,jj) = 0._wp
+         sfx_sni(ji,jj) = 0._wp   ;   sfx_opw(ji,jj) = 0._wp
+         sfx_bog(ji,jj) = 0._wp   ;   sfx_dyn(ji,jj) = 0._wp
+         sfx_bom(ji,jj) = 0._wp   ;   sfx_sum(ji,jj) = 0._wp
+         sfx_res(ji,jj) = 0._wp   ;   sfx_sub(ji,jj) = 0._wp
+         !
+         wfx_snw(ji,jj) = 0._wp   ;   wfx_ice(ji,jj) = 0._wp
+         wfx_sni(ji,jj) = 0._wp   ;   wfx_opw(ji,jj) = 0._wp
+         wfx_bog(ji,jj) = 0._wp   ;   wfx_dyn(ji,jj) = 0._wp
+         wfx_bom(ji,jj) = 0._wp   ;   wfx_sum(ji,jj) = 0._wp
+         wfx_res(ji,jj) = 0._wp   ;   wfx_sub(ji,jj) = 0._wp
+         wfx_spr(ji,jj) = 0._wp   ;   wfx_lam(ji,jj) = 0._wp  
+         wfx_snw_dyn(ji,jj) = 0._wp ; wfx_snw_sum(ji,jj) = 0._wp
+         wfx_snw_sub(ji,jj) = 0._wp ; wfx_ice_sub(ji,jj) = 0._wp
+         wfx_snw_sni(ji,jj) = 0._wp 
+         wfx_pnd(ji,jj) = 0._wp
+
+         hfx_thd(ji,jj) = 0._wp   ;
+         hfx_snw(ji,jj) = 0._wp   ;   hfx_opw(ji,jj) = 0._wp
+         hfx_bog(ji,jj) = 0._wp   ;   hfx_dyn(ji,jj) = 0._wp
+         hfx_bom(ji,jj) = 0._wp   ;   hfx_sum(ji,jj) = 0._wp
+         hfx_res(ji,jj) = 0._wp   ;   hfx_sub(ji,jj) = 0._wp
+         hfx_spr(ji,jj) = 0._wp   ;   hfx_dif(ji,jj) = 0._wp
+         hfx_err_dif(ji,jj) = 0._wp
+         wfx_err_sub(ji,jj) = 0._wp
+         !
+         diag_heat(ji,jj) = 0._wp ;   diag_sice(ji,jj) = 0._wp
+         diag_vice(ji,jj) = 0._wp ;   diag_vsnw(ji,jj) = 0._wp
+
+         tau_icebfr (ji,jj) = 0._wp   ! landfast ice param only (clem: important to keep the init here)
+         qsb_ice_bot(ji,jj) = 0._wp   ! (needed if ln_icethd=F)
+
+         fhld(ji,jj) = 0._wp   ! needed if ln_icethd=F
+
+         ! for control checks (ln_icediachk)
+         diag_trp_vi(ji,jj) = 0._wp   ;   diag_trp_vs(ji,jj) = 0._wp
+         diag_trp_ei(ji,jj) = 0._wp   ;   diag_trp_es(ji,jj) = 0._wp
+         diag_trp_sv(ji,jj) = 0._wp
+         !
+         diag_adv_mass(ji,jj) = 0._wp
+         diag_adv_salt(ji,jj) = 0._wp
+         diag_adv_heat(ji,jj) = 0._wp
+      END_2D
+
+      DO jl = 1, jpl
+         DO_2D( 1, 1, 1, 1 )
+            ! SIMIP diagnostics
+            t_si       (ji,jj,jl) = rt0     ! temp at the ice-snow interface
+            qcn_ice_bot(ji,jj,jl) = 0._wp
+            qcn_ice_top(ji,jj,jl) = 0._wp   ! conductive fluxes
+            cnd_ice    (ji,jj,jl) = 0._wp   ! effective conductivity at the top of ice/snow (ln_cndflx=T)
+            qcn_ice    (ji,jj,jl) = 0._wp   ! conductive flux (ln_cndflx=T & ln_cndemule=T)
+            qtr_ice_bot(ji,jj,jl) = 0._wp   ! part of solar radiation transmitted through the ice needed at least for outputs
+         END_2D
+      ENDDO
+
    END SUBROUTINE diag_set0
+
+
+   SUBROUTINE diag_trends( kn )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE diag_trends  ***
+      !!
+      !! ** purpose : diagnostics of the trends. Used for conservation purposes
+      !!              and outputs
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kn    ! 1 = after dyn ; 2 = after thermo
+      !!----------------------------------------------------------------------
+      !
+      ! --- trends of heat, salt, mass (used for conservation controls)
+      IF( ln_icediachk .OR. iom_use('hfxdhc') ) THEN
+         !
+         diag_heat(:,:) = diag_heat(:,:) &
+            &             - SUM(SUM( e_i (:,:,1:nlay_i,:) - e_i_b (:,:,1:nlay_i,:), dim=4 ), dim=3 ) * r1_Dt_ice &
+            &             - SUM(SUM( e_s (:,:,1:nlay_s,:) - e_s_b (:,:,1:nlay_s,:), dim=4 ), dim=3 ) * r1_Dt_ice
+         diag_sice(:,:) = diag_sice(:,:) &
+            &             + SUM(     sv_i(:,:,:)          - sv_i_b(:,:,:)                  , dim=3 ) * r1_Dt_ice * rhoi
+         diag_vice(:,:) = diag_vice(:,:) &
+            &             + SUM(     v_i (:,:,:)          - v_i_b (:,:,:)                  , dim=3 ) * r1_Dt_ice * rhoi
+         diag_vsnw(:,:) = diag_vsnw(:,:) &
+            &             + SUM(     v_s (:,:,:)          - v_s_b (:,:,:)                  , dim=3 ) * r1_Dt_ice * rhos
+         !
+         IF( kn == 2 )    CALL iom_put ( 'hfxdhc' , diag_heat )   ! output of heat trend
+         !
+      ENDIF
+      !
+      ! --- trends of concentration (used for simip outputs)
+      IF( iom_use('afxdyn') .OR. iom_use('afxthd') .OR. iom_use('afxtot') ) THEN
+         !
+         diag_aice(:,:) = diag_aice(:,:) + SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice
+         !
+         IF( kn == 1 )   CALL iom_put( 'afxdyn' , diag_aice )                                           ! dyn trend
+         IF( kn == 2 )   CALL iom_put( 'afxthd' , SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_Dt_ice ) ! thermo trend
+         IF( kn == 2 )   CALL iom_put( 'afxtot' , diag_aice )                                           ! total trend
+         !
+      ENDIF
+      !
+   END SUBROUTINE diag_trends
 
 #else
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icetab.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icetab.F90	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icetab.F90	(revision 13877)
@@ -40,5 +40,5 @@
       INTEGER , DIMENSION(ndim1d)     , INTENT(in   ) ::   tab_ind  ! input index
       REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(in   ) ::   tab2d    ! input 2D field
-      REAL(wp), DIMENSION(ndim1d,jpl) , INTENT(  out) ::   tab1d    ! output 1D field
+      REAL(wp), DIMENSION(ndim1d,jpl) , INTENT(inout) ::   tab1d    ! output 1D field
       !
       INTEGER ::   jl, jn, jid, jjd
@@ -61,5 +61,5 @@
       INTEGER , DIMENSION(ndim1d) , INTENT(in   ) ::   tab_ind  ! input index
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   tab2d    ! input 2D field
-      REAL(wp), DIMENSION(ndim1d) , INTENT(  out) ::   tab1d    ! output 1D field
+      REAL(wp), DIMENSION(ndim1d) , INTENT(inout) ::   tab1d    ! output 1D field
       !
       INTEGER ::   jn , jid, jjd
@@ -80,5 +80,5 @@
       INTEGER , DIMENSION(ndim1d)     , INTENT(in   ) ::   tab_ind   ! input index
       REAL(wp), DIMENSION(ndim1d,jpl) , INTENT(in   ) ::   tab1d     ! input 1D field
-      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   tab2d     ! output 2D field
+      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(inout) ::   tab2d     ! output 2D field
       !
       INTEGER ::   jl, jn, jid, jjd
@@ -101,5 +101,5 @@
       INTEGER , DIMENSION(ndim1d) , INTENT(in   ) ::   tab_ind   ! input index
       REAL(wp), DIMENSION(ndim1d) , INTENT(in   ) ::   tab1d     ! input 1D field
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   tab2d     ! output 2D field
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   tab2d     ! output 2D field
       !
       INTEGER ::   jn , jid, jjd
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icethd.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icethd.F90	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icethd.F90	(revision 13877)
@@ -18,5 +18,5 @@
    USE ice            ! sea-ice: variables
 !!gm list trop longue ==>>> why not passage en argument d'appel ?
-   USE sbc_oce , ONLY : sss_m, sst_m, e3t_m, utau, vtau, ssu_m, ssv_m, frq_m, qns_tot, qsr_tot, sprecip, ln_cpl
+   USE sbc_oce , ONLY : sss_m, sst_m, e3t_m, utau, vtau, ssu_m, ssv_m, frq_m, sprecip, ln_cpl
    USE sbc_ice , ONLY : qsr_oce, qns_oce, qemp_oce, qsr_ice, qns_ice, dqns_ice, evap_ice, qprec_ice, qevap_ice, &
       &                 qml_ice, qcn_ice, qtr_ice_top
@@ -30,4 +30,5 @@
    USE icethd_pnd     ! sea-ice: melt ponds
    USE iceitd         ! sea-ice: remapping thickness distribution
+   USE icecor         ! sea-ice: corrections
    USE icetab         ! sea-ice: 1D <==> 2D transformation
    USE icevar         ! sea-ice: operations
@@ -52,5 +53,5 @@
    LOGICAL ::   ln_icedO         ! activate ice growth in open-water (T) or not (F)
    LOGICAL ::   ln_icedS         ! activate gravity drainage and flushing (T) or not (F)
-   LOGICAL ::   ln_leadhfx       !  heat in the leads is used to melt sea-ice before warming the ocean
+   LOGICAL ::   ln_leadhfx       ! heat in the leads is used to melt sea-ice before warming the ocean
 
    !! for convergence tests
@@ -91,8 +92,8 @@
       !
       INTEGER  :: ji, jj, jk, jl   ! dummy loop indices
-      REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg
-      REAL(wp), PARAMETER :: zfric_umin = 0._wp           ! lower bound for the friction velocity (cice value=5.e-04)
-      REAL(wp), PARAMETER :: zch        = 0.0057_wp       ! heat transfer coefficient
-      REAL(wp), DIMENSION(jpi,jpj) ::   zu_io, zv_io, zfric   ! ice-ocean velocity (m/s) and frictional velocity (m2/s2)
+      REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg, zqfr_pos
+      REAL(wp), PARAMETER :: zfric_umin = 0._wp       ! lower bound for the friction velocity (cice value=5.e-04)
+      REAL(wp), PARAMETER :: zch        = 0.0057_wp   ! heat transfer coefficient
+      REAL(wp), DIMENSION(jpi,jpj) ::   zu_io, zv_io, zfric, zvel   ! ice-ocean velocity (m/s) and frictional velocity (m2/s2)
       !
       !!-------------------------------------------------------------------
@@ -124,4 +125,6 @@
                &                    (  zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj)   &
                &                     + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1) ) ) * tmask(ji,jj,1)
+            zvel(ji,jj) = 0.5_wp * SQRT( ( u_ice(ji-1,jj) + u_ice(ji,jj) ) * ( u_ice(ji-1,jj) + u_ice(ji,jj) ) + &
+               &                         ( v_ice(ji,jj-1) + v_ice(ji,jj) ) * ( v_ice(ji,jj-1) + v_ice(ji,jj) ) )
          END_2D
       ELSE      !  if no ice dynamics => transmit directly the atmospheric stress to the ocean
@@ -130,10 +133,11 @@
                &                         (  utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj)   &
                &                          + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) ) * tmask(ji,jj,1)
+            zvel(ji,jj) = 0._wp
          END_2D
       ENDIF
 #if defined key_mpi3
-      CALL lbc_lnk_nc_multi( 'icethd', zfric, 'T',  1.0_wp )
+      CALL lbc_lnk_nc_multi( 'icethd', zfric, 'T',  1.0_wp, zvel, 'T', 1.0_wp )
 #else
-      CALL lbc_lnk( 'icethd', zfric, 'T',  1.0_wp )
+      CALL lbc_lnk_multi( 'icethd', zfric, 'T',  1.0_wp, zvel, 'T', 1.0_wp )
 #endif
       !
@@ -144,10 +148,4 @@
          rswitch  = tmask(ji,jj,1) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice
          !
-         !           !  solar irradiance transmission at the mixed layer bottom and used in the lead heat budget
-         !           !  practically no "direct lateral ablation"
-         !           
-         !           !  net downward heat flux from the ice to the ocean, expressed as a function of ocean 
-         !           !  temperature and turbulent mixing (McPhee, 1992)
-         !
          ! --- Energy received in the lead from atm-oce exchanges, zqld is defined everywhere (J.m-2) --- !
          zqld =  tmask(ji,jj,1) * rDt_ice *  &
@@ -155,34 +153,56 @@
             &      ( 1._wp - at_i_b(ji,jj) ) * qns_oce(ji,jj) + qemp_oce(ji,jj) )
 
-         ! --- Energy needed to bring ocean surface layer until its freezing (mostly<0 but >0 if supercooling, J.m-2) --- !
+         ! --- Energy needed to bring ocean surface layer until its freezing, zqfr is defined everywhere (J.m-2) --- !
+         !     (mostly<0 but >0 if supercooling)
          zqfr     = rho0 * rcp * e3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) * tmask(ji,jj,1)  ! both < 0 (t_bo < sst) and > 0 (t_bo > sst)
          zqfr_neg = MIN( zqfr , 0._wp )                                                                    ! only < 0
-
-         ! --- Sensible ocean-to-ice heat flux (mostly>0 but <0 if supercooling, W/m2)
+         zqfr_pos = MAX( zqfr , 0._wp )                                                                    ! only > 0
+
+         ! --- Sensible ocean-to-ice heat flux (W/m2) --- !
+         !     (mostly>0 but <0 if supercooling)
          zfric_u            = MAX( SQRT( zfric(ji,jj) ), zfric_umin ) 
-         qsb_ice_bot(ji,jj) = rswitch * rho0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ! W.m-2
-
-         qsb_ice_bot(ji,jj) = rswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) )
+         qsb_ice_bot(ji,jj) = rswitch * rho0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) )
+         
          ! upper bound for qsb_ice_bot: the heat retrieved from the ocean must be smaller than the heat necessary to reach 
          !                              the freezing point, so that we do not have SST < T_freeze
-         !                              This implies: - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rtdice ) - zqfr >= 0
-
-         !-- Energy Budget of the leads (J.m-2), source of ice growth in open water. Must be < 0 to form ice
-         qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rDt_ice ) - zqfr )
-
-         ! If there is ice and leads are warming => transfer energy from the lead budget and use it for bottom melting 
-         ! If the grid cell is fully covered by ice (no leads) => transfer energy from the lead budget to the ice bottom budget
-         IF( ( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) .OR. at_i(ji,jj) >= (1._wp - epsi10) ) THEN
-            IF( ln_leadhfx ) THEN   ;   fhld(ji,jj) = rswitch * zqld * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90
-            ELSE                    ;   fhld(ji,jj) = 0._wp
-            ENDIF
+         !                              This implies: qsb_ice_bot(ji,jj) * at_i(ji,jj) * rtdice <= - zqfr_neg
+         !                              The following formulation is ok for both normal conditions and supercooling
+         qsb_ice_bot(ji,jj) = rswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) )
+
+         ! --- Energy Budget of the leads (qlead, J.m-2) --- !
+         !     qlead is the energy received from the atm. in the leads.
+         !     If warming (zqld >= 0), then the energy in the leads is used to melt ice (bottom melting) => fhld  (W/m2)
+         !     If cooling (zqld <  0), then the energy in the leads is used to grow ice in open water    => qlead (J.m-2)
+         IF( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) THEN
+            ! upper bound for fhld: fhld should be equal to zqld
+            !                        but we have to make sure that this heat will not make the sst drop below the freezing point
+            !                        so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr_pos
+            !                        The following formulation is ok for both normal conditions and supercooling
+            fhld (ji,jj) = rswitch * MAX( 0._wp, ( zqld - zqfr_pos ) * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) &  ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90
+               &                                 - qsb_ice_bot(ji,jj) )
             qlead(ji,jj) = 0._wp
          ELSE
             fhld (ji,jj) = 0._wp
+            ! upper bound for qlead: qlead should be equal to zqld
+            !                        but before using this heat for ice formation, we suppose that the ocean cools down till the freezing point.
+            !                        The energy for this cooling down is zqfr. Also some heat will be removed from the ocean from turbulent fluxes (qsb)
+            !                        and freezing point is reached if zqfr = zqld - qsb*a/dt
+            !                        so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr
+            !                        The following formulation is ok for both normal conditions and supercooling
+            qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rDt_ice ) - zqfr )
          ENDIF
          !
-         ! Net heat flux on top of the ice-ocean [W.m-2]
-         ! ---------------------------------------------
-         qt_atm_oi(ji,jj) = qns_tot(ji,jj) + qsr_tot(ji,jj) 
+         ! If ice is landfast and ice concentration reaches its max
+         ! => stop ice formation in open water
+         IF(  zvel(ji,jj) <= 5.e-04_wp .AND. at_i(ji,jj) >= rn_amax_2d(ji,jj)-epsi06 )   qlead(ji,jj) = 0._wp
+         !
+         ! If the grid cell is almost fully covered by ice (no leads)
+         ! => stop ice formation in open water
+         IF( at_i(ji,jj) >= (1._wp - epsi10) )   qlead(ji,jj) = 0._wp
+         !
+         ! If ln_leadhfx is false
+         ! => do not use energy of the leads to melt sea-ice
+         IF( .NOT.ln_leadhfx )   fhld(ji,jj) = 0._wp
+         !
       END_2D
       
@@ -195,15 +215,4 @@
       ENDIF
 
-      ! ---------------------------------------------------------------------
-      ! Net heat flux on top of the ocean after ice thermo (1st step) [W.m-2]
-      ! ---------------------------------------------------------------------
-      !     First  step here              :  non solar + precip - qlead - qsensible
-      !     Second step in icethd_dh      :  heat remaining if total melt (zq_rema) 
-      !     Third  step in iceupdate.F90  :  heat from ice-ocean mass exchange (zf_mass) + solar
-      qt_oce_ai(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + qemp_oce(:,:)  &  ! Non solar heat flux received by the ocean               
-         &             - qlead(:,:) * r1_Dt_ice                                &  ! heat flux taken from the ocean where there is open water ice formation
-         &             - at_i (:,:) * qsb_ice_bot(:,:)                         &  ! heat flux taken by sensible flux
-         &             - at_i (:,:) * fhld       (:,:)                            ! heat flux taken during bottom growth/melt 
-      !                                                                           !    (fhld should be 0 while bott growth)
       !-------------------------------------------------------------------------------------------!
       ! Thermodynamic computation (only on grid points covered by ice) => loop over ice categories
@@ -259,4 +268,8 @@
       IF( ln_icedO )          CALL ice_thd_do                       ! --- Frazil ice growth in leads --- !
       !
+                              CALL ice_cor( kt , 2 )                ! --- Corrections --- !
+      !
+      oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rdt_ice              ! ice natural aging incrementation     
+      !
       ! convergence tests
       IF( ln_zdf_chkcvg ) THEN
@@ -422,9 +435,9 @@
          CALL tab_2d_1d( npti, nptidx(1:npti), hfx_res_1d    (1:npti), hfx_res       )
          CALL tab_2d_1d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif   )
-         CALL tab_2d_1d( npti, nptidx(1:npti), qt_oce_ai_1d  (1:npti), qt_oce_ai     )
          !
          ! ocean surface fields
          CALL tab_2d_1d( npti, nptidx(1:npti), sst_1d(1:npti), sst_m )
          CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d(1:npti), sss_m )
+         CALL tab_2d_1d( npti, nptidx(1:npti), frq_m_1d(1:npti), frq_m )
          !
          ! to update ice age
@@ -514,5 +527,4 @@
          CALL tab_1d_2d( npti, nptidx(1:npti), hfx_res_1d    (1:npti), hfx_res     )
          CALL tab_1d_2d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif )
-         CALL tab_1d_2d( npti, nptidx(1:npti), qt_oce_ai_1d  (1:npti), qt_oce_ai   )
          !
          CALL tab_1d_2d( npti, nptidx(1:npti), qns_ice_1d    (1:npti), qns_ice    (:,:,kl) )
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icethd_dh.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icethd_dh.F90	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icethd_dh.F90	(revision 13877)
@@ -139,5 +139,5 @@
       !
       DO ji = 1, npti
-         zf_tt(ji)         = qcn_ice_bot_1d(ji) + qsb_ice_bot_1d(ji) + fhld_1d(ji) 
+         zf_tt(ji)         = qcn_ice_bot_1d(ji) + qsb_ice_bot_1d(ji) + fhld_1d(ji) + qtr_ice_bot_1d(ji) * frq_m_1d(ji) 
          zq_bot(ji)        = MAX( 0._wp, zf_tt(ji) * rDt_ice )
       END DO
@@ -556,5 +556,5 @@
          !    
          ! Remaining heat flux (W.m-2) is sent to the ocean heat budget
-         qt_oce_ai_1d(ji) = qt_oce_ai_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_Dt_ice
+         !!hfx_res_1d(ji) = hfx_res_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_Dt_ice
 
          IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji)
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icethd_do.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icethd_do.F90	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/icethd_do.F90	(revision 13877)
@@ -131,6 +131,8 @@
 
       ! Default new ice thickness
-      WHERE( qlead(:,:) < 0._wp  .AND. tau_icebfr(:,:) == 0._wp )   ;   ht_i_new(:,:) = rn_hinew ! if cooling and no landfast
-      ELSEWHERE                                                     ;   ht_i_new(:,:) = 0._wp
+      WHERE( qlead(:,:) < 0._wp ) ! cooling
+         ht_i_new(:,:) = rn_hinew
+      ELSEWHERE
+         ht_i_new(:,:) = 0._wp
       END WHERE
 
@@ -146,5 +148,5 @@
          !
          DO_2D( 0, 0, 0, 0 )
-            IF ( qlead(ji,jj) < 0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN ! activated if cooling and no landfast
+            IF ( qlead(ji,jj) < 0._wp ) THEN ! cooling
                ! -- Wind stress -- !
                ztaux         = ( utau_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)   &
@@ -202,10 +204,10 @@
       ! 2) Compute thickness, salinity, enthalpy, age, area and volume of new ice
       !------------------------------------------------------------------------------!
-      ! This occurs if open water energy budget is negative (cooling) and there is no landfast ice
+      ! it occurs if cooling
 
       ! Identify grid points where new ice forms
       npti = 0   ;   nptidx(:) = 0
       DO_2D( 1, 1, 1, 1 )
-         IF ( qlead(ji,jj)  <  0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN
+         IF ( qlead(ji,jj)  <  0._wp ) THEN
             npti = npti + 1
             nptidx( npti ) = (jj - 1) * jpi + ji
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/iceupdate.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/iceupdate.F90	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/ICE/iceupdate.F90	(revision 13877)
@@ -24,5 +24,4 @@
    USE traqsr         ! add penetration of solar flux in the calculation of heat budget
    USE icectl         ! sea-ice: control prints
-   USE bdy_oce , ONLY : ln_bdy
    USE zdfdrg  , ONLY : ln_drgice_imp
    !
@@ -92,5 +91,4 @@
       !
       INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
-      REAL(wp) ::   zqmass           ! Heat flux associated with mass exchange ice->ocean (W.m-2)
       REAL(wp) ::   zqsr             ! New solar flux received by the ocean
       REAL(wp), DIMENSION(jpi,jpj) ::   z2d                  ! 2D workspace
@@ -103,4 +101,8 @@
          WRITE(numout,*)'~~~~~~~~~~~~~~'
       ENDIF
+
+      ! Net heat flux on top of the ice-ocean (W.m-2)
+      !----------------------------------------------
+      qt_atm_oi(:,:) = qns_tot(:,:) + qsr_tot(:,:) 
 
       ! --- case we bypass ice thermodynamics --- !
@@ -115,5 +117,5 @@
       DO_2D( 1, 1, 1, 1 )
 
-         ! Solar heat flux reaching the ocean = zqsr (W.m-2) 
+         ! Solar heat flux reaching the ocean (max) = zqsr (W.m-2) 
          !---------------------------------------------------
          zqsr = qsr_tot(ji,jj) - SUM( a_i_b(ji,jj,:) * ( qsr_ice(ji,jj,:) - qtr_ice_bot(ji,jj,:) ) )
@@ -121,17 +123,27 @@
          ! Total heat flux reaching the ocean = qt_oce_ai (W.m-2) 
          !---------------------------------------------------
-         zqmass           = hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) ! heat flux from snow is 0 (T=0 degC)
-         qt_oce_ai(ji,jj) = qt_oce_ai(ji,jj) + zqmass + zqsr
-
-         ! Add the residual from heat diffusion equation and sublimation (W.m-2)
-         !----------------------------------------------------------------------
-         qt_oce_ai(ji,jj) = qt_oce_ai(ji,jj) + hfx_err_dif(ji,jj) +   &
-            &             ( hfx_sub(ji,jj) - SUM( qevap_ice(ji,jj,:) * a_i_b(ji,jj,:) ) )
-
+         qt_oce_ai(ji,jj) = qt_atm_oi(ji,jj) - hfx_sum(ji,jj) - hfx_bom(ji,jj) - hfx_bog(ji,jj) &
+            &                                - hfx_dif(ji,jj) - hfx_opw(ji,jj) - hfx_snw(ji,jj) &
+            &                                + hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_res(ji,jj) &
+            &                                + hfx_sub(ji,jj) - SUM( qevap_ice(ji,jj,:) * a_i_b(ji,jj,:) ) + hfx_spr(ji,jj)                 
+         
          ! New qsr and qns used to compute the oceanic heat flux at the next time step
          !----------------------------------------------------------------------------
-         qsr(ji,jj) = zqsr                                      
+         ! if warming and some ice remains, then we suppose that the whole solar flux has been consumed to melt the ice
+         ! else ( cooling or no ice left ), then we suppose that     no    solar flux has been consumed
+         !
+         IF( fhld(ji,jj) > 0._wp .AND. at_i(ji,jj) > 0._wp ) THEN   !-- warming and some ice remains
+            !                                        solar flux transmitted thru the 1st level of the ocean (i.e. not used by sea-ice)
+            qsr(ji,jj) = ( 1._wp - at_i_b(ji,jj) ) * qsr_oce(ji,jj) * ( 1._wp - frq_m(ji,jj) ) &
+               !                                   + solar flux transmitted thru ice and the 1st ocean level (also not used by sea-ice)
+               &             + SUM( a_i_b(ji,jj,:) * qtr_ice_bot(ji,jj,:) ) * ( 1._wp - frq_m(ji,jj) )
+            !
+         ELSE                                                       !-- cooling or no ice left
+            qsr(ji,jj) = zqsr
+         ENDIF
+         !
+         ! the non-solar is simply derived from the solar flux
          qns(ji,jj) = qt_oce_ai(ji,jj) - zqsr              
-
+         
          ! Mass flux at the atm. surface       
          !-----------------------------------
@@ -140,23 +152,14 @@
          ! Mass flux at the ocean surface      
          !------------------------------------
-         !  case of realistic freshwater flux (Tartinville et al., 2001) (presently ACTIVATED)
-         !  ------------------------------------------------------------------------------------- 
-         !  The idea of this approach is that the system that we consider is the ICE-OCEAN system
-         !  Thus  FW  flux  =  External ( E-P+snow melt)
-         !       Salt flux  =  Exchanges in the ice-ocean system then converted into FW
-         !                     Associated to Ice formation AND Ice melting
-         !                     Even if i see Ice melting as a FW and SALT flux
-         !        
-         ! mass flux from ice/ocean
+         ! ice-ocean  mass flux
          wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj)   &
             &           + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) + wfx_pnd(ji,jj)
-
-         ! add the snow melt water to snow mass flux to the ocean
+         
+         ! snw-ocean mass flux
          wfx_snw(ji,jj) = wfx_snw_sni(ji,jj) + wfx_snw_dyn(ji,jj) + wfx_snw_sum(ji,jj)
-
-         ! mass flux at the ocean/ice interface
-         fmmflx(ji,jj) = - ( wfx_ice(ji,jj) + wfx_snw(ji,jj) + wfx_err_sub(ji,jj) )              ! F/M mass flux save at least for biogeochemical model
-         emp(ji,jj)    = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_err_sub(ji,jj)   ! mass flux + F/M mass flux (always ice/ocean mass exchange)
-
+         
+         ! total mass flux at the ocean/ice interface
+         fmmflx(ji,jj) =                - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_err_sub(ji,jj)   ! ice-ocean mass flux saved at least for biogeochemical model
+         emp   (ji,jj) = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_err_sub(ji,jj)   ! atm-ocean + ice-ocean mass flux
 
          ! Salt flux at the ocean surface      
@@ -262,5 +265,5 @@
       CALL iom_put ('hfxdif'     , hfx_dif     )   ! heat flux used for ice temperature change
       CALL iom_put ('hfxsnw'     , hfx_snw     )   ! heat flux used for snow melt 
-      CALL iom_put ('hfxerr'     , hfx_err_dif )   ! heat flux error after heat diffusion (included in qt_oce_ai)
+      CALL iom_put ('hfxerr'     , hfx_err_dif )   ! heat flux error after heat diffusion
 
       ! heat fluxes associated with mass exchange (freeze/melt/precip...)
@@ -279,9 +282,9 @@
       !---------
 #if ! defined key_agrif
-      IF( ln_icediachk .AND. .NOT. ln_bdy)   CALL ice_cons_final('iceupdate')                                       ! conservation
+      IF( ln_icediachk      )   CALL ice_cons_final('iceupdate')                                       ! conservation
 #endif
-      IF( ln_icectl                      )   CALL ice_prt       (kt, iiceprt, jiceprt, 3, 'Final state ice_update') ! prints
-      IF( sn_cfctl%l_prtctl              )   CALL ice_prt3D     ('iceupdate')                                       ! prints
-      IF( ln_timing                      )   CALL timing_stop   ('ice_update')                                      ! timing
+      IF( ln_icectl         )   CALL ice_prt       (kt, iiceprt, jiceprt, 3, 'Final state ice_update') ! prints
+      IF( sn_cfctl%l_prtctl )   CALL ice_prt3D     ('iceupdate')                                       ! prints
+      IF( ln_timing         )   CALL timing_stop   ('ice_update')                                      ! timing
       !
    END SUBROUTINE ice_update_flx
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/BDY/bdyice.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/BDY/bdyice.F90	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/BDY/bdyice.F90	(revision 13877)
@@ -61,7 +61,5 @@
       !!----------------------------------------------------------------------
       ! controls
-      IF( ln_timing    )   CALL timing_start('bdy_ice_thd')                                                            ! timing
-      IF( ln_icediachk )   CALL ice_cons_hsm(0,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
-      IF( ln_icediachk )   CALL ice_cons2D  (0,'bdy_ice_thd',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation
+      IF( ln_timing )   CALL timing_start('bdy_ice_thd')   ! timing
       !
       CALL ice_var_glo2eqv
@@ -120,8 +118,6 @@
       !
       ! controls
-      IF( ln_icectl    )   CALL ice_prt     ( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' )                        ! prints
-      IF( ln_icediachk )   CALL ice_cons_hsm(1,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
-      IF( ln_icediachk )   CALL ice_cons2D  (1,'bdy_ice_thd',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation
-      IF( ln_timing    )   CALL timing_stop ('bdy_ice_thd')                                                            ! timing
+      IF( ln_icectl )   CALL ice_prt     ( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' )   ! prints
+      IF( ln_timing )   CALL timing_stop ('bdy_ice_thd')                                       ! timing
       !
    END SUBROUTINE bdy_ice
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/DIU/diu_coolskin.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/DIU/diu_coolskin.F90	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/DIU/diu_coolskin.F90	(revision 13877)
@@ -95,5 +95,5 @@
       !!----------------------------------------------------------------------
       !
-      IF( .NOT. ln_blk )   CALL ctl_stop("diu_coolskin.f90: diurnal flux processing only implemented for bulk forcing")
+      IF( .NOT. (ln_blk .OR. ln_abl) )   CALL ctl_stop("diu_coolskin.f90: diurnal flux processing only implemented for bulk forcing")
       !
       DO_2D( 1, 1, 1, 1 )
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/IOM/iom.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/IOM/iom.F90	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/IOM/iom.F90	(revision 13877)
@@ -1935,9 +1935,5 @@
       IF( iom_use(cdname) ) THEN
 #if defined key_iomput
-         IF( SIZE(pfield2d, dim=1) == jpi .AND. SIZE(pfield2d, dim=2) == jpj ) THEN
-            CALL xios_send_field( cdname, pfield2d(Nis0:Nie0, Njs0:Nje0) )       ! this extraction will create a copy of pfield2d
-         ELSE
-            CALL xios_send_field( cdname, pfield2d )
-         ENDIF
+         CALL xios_send_field( cdname, pfield2d )
 #else
          WRITE(numout,*) pfield2d   ! iom_use(cdname) = .F. -> useless test to avoid compilation warnings
@@ -1951,9 +1947,5 @@
       IF( iom_use(cdname) ) THEN
 #if defined key_iomput
-         IF( SIZE(pfield2d, dim=1) == jpi .AND. SIZE(pfield2d, dim=2) == jpj ) THEN
-            CALL xios_send_field( cdname, pfield2d(Nis0:Nie0, Njs0:Nje0) )       ! this extraction will create a copy of pfield2d
-         ELSE
-            CALL xios_send_field( cdname, pfield2d )
-         ENDIF
+         CALL xios_send_field( cdname, pfield2d )
 #else
          WRITE(numout,*) pfield2d   ! iom_use(cdname) = .F. -> useless test to avoid compilation warnings
@@ -1967,9 +1959,5 @@
       IF( iom_use(cdname) ) THEN
 #if defined key_iomput
-         IF( SIZE(pfield3d, dim=1) == jpi .AND. SIZE(pfield3d, dim=2) == jpj ) THEN
-            CALL xios_send_field( cdname, pfield3d(Nis0:Nie0, Njs0:Nje0,:) )     ! this extraction will create a copy of pfield3d
-         ELSE
-            CALL xios_send_field( cdname, pfield3d )
-         ENDIF
+         CALL xios_send_field( cdname, pfield3d )
 #else
          WRITE(numout,*) pfield3d   ! iom_use(cdname) = .F. -> useless test to avoid compilation warnings
@@ -1983,9 +1971,5 @@
       IF( iom_use(cdname) ) THEN
 #if defined key_iomput
-         IF( SIZE(pfield3d, dim=1) == jpi .AND. SIZE(pfield3d, dim=2) == jpj ) THEN
-            CALL xios_send_field( cdname, pfield3d(Nis0:Nie0, Njs0:Nje0,:) )     ! this extraction will create a copy of pfield3d
-         ELSE
-            CALL xios_send_field( cdname, pfield3d )
-         ENDIF
+         CALL xios_send_field( cdname, pfield3d )
 #else
          WRITE(numout,*) pfield3d   ! iom_use(cdname) = .F. -> useless test to avoid compilation warnings
@@ -1999,9 +1983,5 @@
       IF( iom_use(cdname) ) THEN
 #if defined key_iomput
-         IF( SIZE(pfield4d, dim=1) == jpi .AND. SIZE(pfield4d, dim=2) == jpj ) THEN
-            CALL xios_send_field( cdname, pfield4d(Nis0:Nie0, Njs0:Nje0,:,:) )   ! this extraction will create a copy of pfield4d
-         ELSE
-            CALL xios_send_field (cdname, pfield4d )
-         ENDIF
+         CALL xios_send_field (cdname, pfield4d )
 #else
          WRITE(numout,*) pfield4d   ! iom_use(cdname) = .F. -> useless test to avoid compilation warnings
@@ -2015,9 +1995,5 @@
       IF( iom_use(cdname) ) THEN
 #if defined key_iomput
-         IF( SIZE(pfield4d, dim=1) == jpi .AND. SIZE(pfield4d, dim=2) == jpj ) THEN
-            CALL xios_send_field( cdname, pfield4d(Nis0:Nie0, Njs0:Nje0,:,:) )   ! this extraction will create a copy of pfield4d
-         ELSE
-            CALL xios_send_field (cdname, pfield4d )
-         ENDIF
+         CALL xios_send_field (cdname, pfield4d )
 #else
          WRITE(numout,*) pfield4d   ! iom_use(cdname) = .F. -> useless test to avoid compilation warnings
@@ -2225,5 +2201,5 @@
       !
       CALL iom_set_domain_attr("grid_"//cdgrd, ni_glo=Ni0glo,nj_glo=Nj0glo,ibegin=mig0(Nis0)-1,jbegin=mjg0(Njs0)-1,ni=Ni_0,nj=Nj_0)
-      CALL iom_set_domain_attr("grid_"//cdgrd, data_dim=2, data_ibegin = 0, data_ni = Ni_0, data_jbegin = 0, data_nj = Nj_0)
+      CALL iom_set_domain_attr("grid_"//cdgrd, data_dim=2, data_ibegin = -nn_hls, data_ni = jpi, data_jbegin = -nn_hls, data_nj = jpj)
 !don't define lon and lat for restart reading context. 
       IF ( .NOT.ldrxios ) &
@@ -2328,8 +2304,8 @@
       CALL dom_ngb( 180.0_wp, 90.0_wp, ix, iy, 'T' ) !  i-line that passes near the North Pole : Reference latitude (used in plots)
       CALL iom_set_domain_attr("gznl", ni_glo=Ni0glo, nj_glo=Nj0glo, ibegin=mig0(Nis0)-1, jbegin=mjg0(Njs0)-1, ni=Ni_0, nj=Nj_0)
-      CALL iom_set_domain_attr("gznl", data_dim=2, data_ibegin = 0, data_ni = Ni_0, data_jbegin = 0, data_nj = Nj_0)
+      CALL iom_set_domain_attr("gznl", data_dim=2, data_ibegin = -nn_hls, data_ni = jpi, data_jbegin = -nn_hls, data_nj = jpj)
       CALL iom_set_domain_attr("gznl", lonvalue = real(zlon, dp),   &
          &                             latvalue = real(RESHAPE(plat(Nis0:Nie0, Njs0:Nje0),(/ Ni_0*Nj_0 /)),dp))  
-      CALL iom_set_zoom_domain_attr("ptr", ibegin=ix-1, jbegin=0, ni=1, nj=Nj_0)
+      CALL iom_set_zoom_domain_attr("ptr", ibegin=ix-1, jbegin=0, ni=1, nj=Nj0glo)
       !
       CALL iom_update_file_name('ptr')
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/LBC/lib_mpp.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/LBC/lib_mpp.F90	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/LBC/lib_mpp.F90	(revision 13877)
@@ -516,4 +516,5 @@
             ALLOCATE(todelay(idvar)%y1d(isz))
             todelay(idvar)%y1d(:) = CMPLX(todelay(idvar)%z1d(:), 0., wp)   ! create %y1d, complex variable needed by mpi_sumdd
+            ndelayid(idvar) = MPI_REQUEST_NULL                             ! initialised request to a valid value
          END IF
       ENDIF
@@ -523,8 +524,8 @@
          ALLOCATE(todelay(idvar)%z1d(isz), todelay(idvar)%y1d(isz))   ! allocate also %z1d as used for the restart
          CALL mpi_allreduce( y_in(:), todelay(idvar)%y1d(:), isz, MPI_DOUBLE_COMPLEX, mpi_sumdd, ilocalcomm, ierr )   ! get %y1d
-         todelay(idvar)%z1d(:) = REAL(todelay(idvar)%y1d(:), wp)      ! define %z1d from %y1d
-      ENDIF
-
-      IF( ndelayid(idvar) > 0 )   CALL mpp_delay_rcv( idvar )         ! make sure %z1d is received
+         ndelayid(idvar) = MPI_REQUEST_NULL
+      ENDIF
+
+      CALL mpp_delay_rcv( idvar )         ! make sure %z1d is received
 
       ! send back pout from todelay(idvar)%z1d defined at previous call
@@ -535,5 +536,5 @@
       IF( ln_timing ) CALL tic_tac( .TRUE., ld_global = .TRUE.)
       CALL  mpi_allreduce( y_in(:), todelay(idvar)%y1d(:), isz, MPI_DOUBLE_COMPLEX, mpi_sumdd, ilocalcomm, ierr )
-      ndelayid(idvar) = 1
+      ndelayid(idvar) = MPI_REQUEST_NULL
       IF( ln_timing ) CALL tic_tac(.FALSE., ld_global = .TRUE.)
 # else
@@ -596,4 +597,6 @@
             DEALLOCATE(todelay(idvar)%z1d)
             ndelayid(idvar) = -1                                      ! do as if we had no restart
+         ELSE
+            ndelayid(idvar) = MPI_REQUEST_NULL
          END IF
       ENDIF
@@ -603,7 +606,8 @@
          ALLOCATE(todelay(idvar)%z1d(isz))
          CALL mpi_allreduce( p_in(:), todelay(idvar)%z1d(:), isz, MPI_DOUBLE_PRECISION, mpi_max, ilocalcomm, ierr )   ! get %z1d
-      ENDIF
-
-      IF( ndelayid(idvar) > 0 )   CALL mpp_delay_rcv( idvar )         ! make sure %z1d is received
+         ndelayid(idvar) = MPI_REQUEST_NULL
+      ENDIF
+
+      CALL mpp_delay_rcv( idvar )         ! make sure %z1d is received
 
       ! send back pout from todelay(idvar)%z1d defined at previous call
@@ -611,7 +615,8 @@
 
       ! send p_in into todelay(idvar)%z1d with a non-blocking communication
+      ! (PM) Should we get rid of MPI2 option ? MPI3 was release in 2013. Who is still using MPI2 ?
 # if defined key_mpi2
       IF( ln_timing ) CALL tic_tac( .TRUE., ld_global = .TRUE.)
-      CALL  mpi_allreduce( p_in(:), todelay(idvar)%z1d(:), isz, MPI_TYPE, mpi_max, ilocalcomm, ndelayid(idvar), ierr )
+      CALL  mpi_allreduce( p_in(:), todelay(idvar)%z1d(:), isz, MPI_TYPE, mpi_max, ilocalcomm, ierr )
       IF( ln_timing ) CALL tic_tac(.FALSE., ld_global = .TRUE.)
 # else
@@ -636,13 +641,9 @@
       !!----------------------------------------------------------------------
 #if defined key_mpp_mpi
-      IF( ndelayid(kid) /= -2 ) THEN  
-#if ! defined key_mpi2
-         IF( ln_timing ) CALL tic_tac( .TRUE., ld_global = .TRUE.)
-         CALL mpi_wait( ndelayid(kid), MPI_STATUS_IGNORE, ierr )                        ! make sure todelay(kid) is received
-         IF( ln_timing ) CALL tic_tac(.FALSE., ld_global = .TRUE.)
-#endif
-         IF( ASSOCIATED(todelay(kid)%y1d) )   todelay(kid)%z1d(:) = REAL(todelay(kid)%y1d(:), wp)  ! define %z1d from %y1d
-         ndelayid(kid) = -2   ! add flag to know that mpi_wait was already called on kid
-      ENDIF
+      IF( ln_timing ) CALL tic_tac( .TRUE., ld_global = .TRUE.)
+      ! test on ndelayid(kid) useless as mpi_wait return immediatly if the request handle is MPI_REQUEST_NULL
+      CALL mpi_wait( ndelayid(kid), MPI_STATUS_IGNORE, ierr ) ! after this ndelayid(kid) = MPI_REQUEST_NULL
+      IF( ln_timing ) CALL tic_tac( .FALSE., ld_global = .TRUE.)
+      IF( ASSOCIATED(todelay(kid)%y1d) )   todelay(kid)%z1d(:) = REAL(todelay(kid)%y1d(:), wp)  ! define %z1d from %y1d
 #endif
    END SUBROUTINE mpp_delay_rcv
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcfwb.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcfwb.F90	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcfwb.F90	(revision 13877)
@@ -94,4 +94,5 @@
          snwice_mass_b(:,:) = 0.e0               ! no sea-ice model is being used : no snow+ice mass
          snwice_mass  (:,:) = 0.e0
+         snwice_fmass (:,:) = 0.e0
 #endif
          !
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcmod.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcmod.F90	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcmod.F90	(revision 13877)
@@ -229,5 +229,5 @@
       CASE DEFAULT                     !- not supported
       END SELECT
-      IF( ln_diurnal .AND. .NOT. ln_blk  )   CALL ctl_stop( "sbc_init: diurnal flux processing only implemented for bulk forcing" )
+      IF( ln_diurnal .AND. .NOT. (ln_blk.OR.ln_abl) )   CALL ctl_stop( "sbc_init: diurnal flux processing only implemented for bulk forcing" )
       !
       !                       !**  allocate and set required variables
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/stpctl.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/stpctl.F90	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/stpctl.F90	(revision 13877)
@@ -67,5 +67,5 @@
       REAL(wp)                        ::   zzz                                   ! local real 
       REAL(wp), DIMENSION(9)          ::   zmax, zmaxlocal
-      LOGICAL                         ::   ll_wrtstp, ll_colruns, ll_wrtruns
+      LOGICAL                         ::   ll_wrtstp, ll_colruns, ll_wrtruns, ll_0oce
       LOGICAL, DIMENSION(jpi,jpj,jpk) ::   llmsk
       CHARACTER(len=20)               ::   clname
@@ -125,4 +125,7 @@
       !
       llmsk(Nis0:Nie0,Njs0:Nje0,1) = ssmask(Nis0:Nie0,Njs0:Nje0) == 1._wp         ! define only the inner domain
+      !
+      ll_0oce = .NOT. ANY( llmsk(:,:,1) )                                         ! no ocean point in the inner domain?
+      !
       IF( ll_wd ) THEN
          zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm) + ssh_ref ), mask = llmsk(:,:,1) )   ! ssh max
@@ -149,25 +152,30 @@
       ENDIF
       zmax(9) = REAL( nstop, wp )                                                 ! stop indicator
+      !
       !                                   !==               get global extrema             ==!
       !                                   !==  done by all processes if writting run.stat  ==!
       IF( ll_colruns ) THEN
          zmaxlocal(:) = zmax(:)
-         CALL mpp_max( "stpctl", zmax )          ! max over the global domain
+         CALL mpp_max( "stpctl", zmax )          ! max over the global domain: ok even of ll_0oce = .true. 
          nstop = NINT( zmax(9) )                 ! update nstop indicator (now sheared among all local domains)
-      ENDIF
+      ELSE
+         ! if no ocean point: MAXVAL returns -HUGE => we must overwrite this value to avoid error handling bellow.
+         IF( ll_0oce )   zmax(1:4) = (/ 0._wp, 0._wp, -1._wp, 1._wp /)   ! default "valid" values...
+      ENDIF
+      !
+      zmax(3) = -zmax(3)                         ! move back from max(-zz) to min(zz) : easier to manage! 
+      zmax(5) = -zmax(5)                         ! move back from max(-zz) to min(zz) : easier to manage!
+      IF( ll_colruns ) THEN
+         zmaxlocal(3) = -zmaxlocal(3)            ! move back from max(-zz) to min(zz) : easier to manage! 
+         zmaxlocal(5) = -zmaxlocal(5)            ! move back from max(-zz) to min(zz) : easier to manage!
+      ENDIF
+      !
       !                                   !==              write "run.stat" files              ==!
       !                                   !==  done only by 1st subdomain at writting timestep  ==!
       IF( ll_wrtruns ) THEN
-         WRITE(numrun,9500) kt, zmax(1), zmax(2), -zmax(3), zmax(4)
-         istatus = NF90_PUT_VAR( nrunid, nvarid(1), (/ zmax(1)/), (/kt/), (/1/) )
-         istatus = NF90_PUT_VAR( nrunid, nvarid(2), (/ zmax(2)/), (/kt/), (/1/) )
-         istatus = NF90_PUT_VAR( nrunid, nvarid(3), (/-zmax(3)/), (/kt/), (/1/) )
-         istatus = NF90_PUT_VAR( nrunid, nvarid(4), (/ zmax(4)/), (/kt/), (/1/) )
-         istatus = NF90_PUT_VAR( nrunid, nvarid(5), (/-zmax(5)/), (/kt/), (/1/) )
-         istatus = NF90_PUT_VAR( nrunid, nvarid(6), (/ zmax(6)/), (/kt/), (/1/) )
-         IF( ln_zad_Aimp ) THEN
-            istatus = NF90_PUT_VAR( nrunid, nvarid(7), (/ zmax(7)/), (/kt/), (/1/) )
-            istatus = NF90_PUT_VAR( nrunid, nvarid(8), (/ zmax(8)/), (/kt/), (/1/) )
-         ENDIF
+         WRITE(numrun,9500) kt, zmax(1), zmax(2), zmax(3), zmax(4)
+         DO ji = 1, 6 + 2 * COUNT( (/ln_zad_Aimp/) )
+            istatus = NF90_PUT_VAR( nrunid, nvarid(ji), (/zmax(ji)/), (/kt/), (/1/) )
+         END DO
          IF( kt == nitend )   istatus = NF90_CLOSE(nrunid)
       END IF
@@ -175,11 +183,11 @@
       !                                   !==  done by all processes at every time step  ==!
       !
-      IF(   zmax(1) >   20._wp .OR.   &                   ! too large sea surface height ( > 20 m )
-         &  zmax(2) >   10._wp .OR.   &                   ! too large velocity ( > 10 m/s)
-         &  zmax(3) >=   0._wp .OR.   &                   ! negative or zero sea surface salinity
-         &  zmax(4) >= 100._wp .OR.   &                   ! too large sea surface salinity ( > 100 )
-         &  zmax(4) <    0._wp .OR.   &                   ! too large sea surface salinity (keep this line for sea-ice)
-         &  ISNAN( zmax(1) + zmax(2) + zmax(3) ) .OR.   &               ! NaN encounter in the tests
-         &  ABS(   zmax(1) + zmax(2) + zmax(3) ) > HUGE(1._wp) ) THEN   ! Infinity encounter in the tests
+      IF(  zmax(1) >   20._wp .OR.   &                   ! too large sea surface height ( > 20 m )
+         & zmax(2) >   10._wp .OR.   &                   ! too large velocity ( > 10 m/s)
+         & zmax(3) <=   0._wp .OR.   &                   ! negative or zero sea surface salinity
+         & zmax(4) >= 100._wp .OR.   &                   ! too large sea surface salinity ( > 100 )
+         & zmax(4) <    0._wp .OR.   &                   ! too large sea surface salinity (keep this line for sea-ice)
+         & ISNAN( zmax(1) + zmax(2) + zmax(3) ) .OR.   &               ! NaN encounter in the tests
+         & ABS(   zmax(1) + zmax(2) + zmax(3) ) > HUGE(1._wp) ) THEN   ! Infinity encounter in the tests
          !
          iloc(:,:) = 0
@@ -221,8 +229,8 @@
          !
          WRITE(ctmp1,*) ' stp_ctl: |ssh| > 20 m  or  |U| > 10 m/s  or  S <= 0  or  S >= 100  or  NaN encounter in the tests'
-         CALL wrt_line( ctmp2, kt, '|ssh| max',  zmax(1), iloc(:,1), iareasum(1), iareamin(1), iareamax(1) )
-         CALL wrt_line( ctmp3, kt, '|U|   max',  zmax(2), iloc(:,2), iareasum(2), iareamin(2), iareamax(2) )
-         CALL wrt_line( ctmp4, kt, 'Sal   min', -zmax(3), iloc(:,3), iareasum(3), iareamin(3), iareamax(3) )
-         CALL wrt_line( ctmp5, kt, 'Sal   max',  zmax(4), iloc(:,4), iareasum(4), iareamin(4), iareamax(4) )
+         CALL wrt_line( ctmp2, kt, '|ssh| max', zmax(1), iloc(:,1), iareasum(1), iareamin(1), iareamax(1) )
+         CALL wrt_line( ctmp3, kt, '|U|   max', zmax(2), iloc(:,2), iareasum(2), iareamin(2), iareamax(2) )
+         CALL wrt_line( ctmp4, kt, 'Sal   min', zmax(3), iloc(:,3), iareasum(3), iareamin(3), iareamax(3) )
+         CALL wrt_line( ctmp5, kt, 'Sal   max', zmax(4), iloc(:,4), iareasum(4), iareamin(4), iareamax(4) )
          IF( Agrif_Root() ) THEN
             WRITE(ctmp6,*) '      ===> output of last computed fields in output.abort* files'
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/SAS/stpctl.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/SAS/stpctl.F90	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/SAS/stpctl.F90	(revision 13877)
@@ -68,5 +68,5 @@
       REAL(wp)                        ::   zzz                                   ! local real 
       REAL(wp), DIMENSION(4)          ::   zmax, zmaxlocal
-      LOGICAL                         ::   ll_wrtstp, ll_colruns, ll_wrtruns
+      LOGICAL                         ::   ll_wrtstp, ll_colruns, ll_wrtruns, ll_0oce
       LOGICAL, DIMENSION(jpi,jpj)     ::   llmsk
       CHARACTER(len=20)               ::   clname
@@ -124,31 +124,34 @@
       llmsk(:,Nje1: jpj) = .FALSE.
       !
-      llmsk(Nis0:Nie0,Njs0:Nje0) = tmask(Nis0:Nie0,Njs0:Nje0,1) == 1._wp   ! test only the inner domain
-      IF( COUNT( llmsk(:,:) ) > 0 ) THEN   ! avoid huge values sent back for land processors...
-         zmax(1) = MAXVAL(      vt_i (:,:)      , mask = llmsk )   ! max ice thickness
-         zmax(2) = MAXVAL( ABS( u_ice(:,:) )    , mask = llmsk )   ! max ice velocity (zonal only)
-         zmax(3) = MAXVAL(     -tm_i (:,:) + rt0, mask = llmsk )   ! min ice temperature (in degC)
-      ELSE
-         IF( ll_colruns ) THEN    ! default value: must not be kept when calling mpp_max -> must be as small as possible
-            zmax(1:3) = -HUGE(1._wp)
-         ELSE                     ! default value: must not give true for any of the tests bellow (-> avoid manipulating HUGE...)
-            zmax(1:3) = 0._wp
-         ENDIF
-      ENDIF
-      zmax(4) = REAL( nstop, wp )                                     ! stop indicator
+      llmsk(Nis0:Nie0,Njs0:Nje0) = tmask(Nis0:Nie0,Njs0:Nje0,1) == 1._wp        ! test only the inner domain
+      !
+      ll_0oce = .NOT. ANY( llmsk(:,:) )                                         ! no ocean point in the inner domain?
+      !
+      zmax(1) = MAXVAL(      vt_i (:,:)      , mask = llmsk )                   ! max ice thickness
+      zmax(2) = MAXVAL( ABS( u_ice(:,:) )    , mask = llmsk )                   ! max ice velocity (zonal only)
+      zmax(3) = MAXVAL(     -tm_i (:,:) + rt0, mask = llmsk )                   ! min ice temperature (in degC)
+      zmax(4) = REAL( nstop, wp )                                               ! stop indicator
+      !
       !                                   !==               get global extrema             ==!
       !                                   !==  done by all processes if writting run.stat  ==!
       IF( ll_colruns ) THEN
          zmaxlocal(:) = zmax(:)
-         CALL mpp_max( "stpctl", zmax )          ! max over the global domain
+         CALL mpp_max( "stpctl", zmax )          ! max over the global domain: ok even of ll_0oce = .true.
          nstop = NINT( zmax(4) )                 ! update nstop indicator (now sheared among all local domains)
-      ENDIF
+      ELSE
+         ! if no ocean point: MAXVAL returns -HUGE => we must overwrite this value to avoid error handling bellow.
+         IF( ll_0oce )   zmax(1:3) = 0._wp       ! default "valid" values...
+      ENDIF
+      !
+      zmax(3) = -zmax(3)                              ! move back from max(-zz) to min(zz) : easier to manage!
+      IF( ll_colruns ) zmaxlocal(3) = -zmaxlocal(3)   ! move back from max(-zz) to min(zz) : easier to manage!
+      !
       !                                   !==              write "run.stat" files              ==!
       !                                   !==  done only by 1st subdomain at writting timestep  ==!
       IF( ll_wrtruns ) THEN
-         WRITE(numrun,9500) kt, zmax(1), zmax(2), -zmax(3)
-         istatus = NF90_PUT_VAR( nrunid, nvarid(1), (/ zmax(1)/), (/kt/), (/1/) )
-         istatus = NF90_PUT_VAR( nrunid, nvarid(2), (/ zmax(2)/), (/kt/), (/1/) )
-         istatus = NF90_PUT_VAR( nrunid, nvarid(3), (/-zmax(3)/), (/kt/), (/1/) )
+         WRITE(numrun,9500) kt, zmax(1), zmax(2), zmax(3)
+         DO ji = 1, 3
+            istatus = NF90_PUT_VAR( nrunid, nvarid(ji), (/zmax(ji)/), (/kt/), (/1/) )
+         END DO
          IF( kt == nitend )   istatus = NF90_CLOSE(nrunid)
       END IF
@@ -158,5 +161,5 @@
       IF(   zmax(1) >  100._wp .OR.   &                   ! too large ice thickness maximum ( > 100 m)
          &  zmax(2) >   10._wp .OR.   &                   ! too large ice velocity ( > 10 m/s)
-         &  zmax(3) >  101._wp .OR.   &                   ! too cold ice temperature ( < -100 degC)
+         &  zmax(3) < -101._wp .OR.   &                   ! too cold ice temperature ( < -100 degC)
          &  ISNAN( zmax(1) + zmax(2) + zmax(3) ) .OR.   &               ! NaN encounter in the tests
          &  ABS(   zmax(1) + zmax(2) + zmax(3) ) > HUGE(1._wp) ) THEN   ! Infinity encounter in the tests
@@ -192,7 +195,7 @@
          !
          WRITE(ctmp1,*) ' stp_ctl: ice_thick > 100 m or |ice_vel| > 10 m/s or ice_temp < -100 degC or NaN encounter in the tests'
-         CALL wrt_line( ctmp2, kt, 'ice_thick max',  zmax(1), iloc(:,1), iareasum(1), iareamin(1), iareamax(1) )
-         CALL wrt_line( ctmp3, kt, '|ice_vel| max',  zmax(2), iloc(:,2), iareasum(2), iareamin(2), iareamax(2) )
-         CALL wrt_line( ctmp4, kt, 'ice_temp  min', -zmax(3), iloc(:,3), iareasum(3), iareamin(3), iareamax(3) )
+         CALL wrt_line( ctmp2, kt, 'ice_thick max', zmax(1), iloc(:,1), iareasum(1), iareamin(1), iareamax(1) )
+         CALL wrt_line( ctmp3, kt, '|ice_vel| max', zmax(2), iloc(:,2), iareasum(2), iareamin(2), iareamax(2) )
+         CALL wrt_line( ctmp4, kt, 'ice_temp  min', zmax(3), iloc(:,3), iareasum(3), iareamin(3), iareamax(3) )
          IF( Agrif_Root() ) THEN
             WRITE(ctmp6,*) '      ===> output of last computed fields in output.abort* files'
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/CANAL/MY_SRC/stpctl.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/CANAL/MY_SRC/stpctl.F90	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/CANAL/MY_SRC/stpctl.F90	(revision 13877)
@@ -67,5 +67,5 @@
       REAL(wp)                        ::   zzz                                   ! local real 
       REAL(wp), DIMENSION(9)          ::   zmax, zmaxlocal
-      LOGICAL                         ::   ll_wrtstp, ll_colruns, ll_wrtruns
+      LOGICAL                         ::   ll_wrtstp, ll_colruns, ll_wrtruns, ll_0oce
       LOGICAL, DIMENSION(jpi,jpj,jpk) ::   llmsk
       CHARACTER(len=20)               ::   clname
@@ -125,4 +125,7 @@
       !
       llmsk(Nis0:Nie0,Njs0:Nje0,1) = ssmask(Nis0:Nie0,Njs0:Nje0) == 1._wp         ! define only the inner domain
+      !
+      ll_0oce = .NOT. ANY( llmsk(:,:,1) )                                         ! no ocean point in the inner domain?
+      !
       IF( ll_wd ) THEN
          zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm) + ssh_ref ), mask = llmsk(:,:,1) )   ! ssh max
@@ -149,25 +152,30 @@
       ENDIF
       zmax(9) = REAL( nstop, wp )                                                 ! stop indicator
+      !
       !                                   !==               get global extrema             ==!
       !                                   !==  done by all processes if writting run.stat  ==!
       IF( ll_colruns ) THEN
          zmaxlocal(:) = zmax(:)
-         CALL mpp_max( "stpctl", zmax )          ! max over the global domain
+         CALL mpp_max( "stpctl", zmax )          ! max over the global domain: ok even of ll_0oce = .true. 
          nstop = NINT( zmax(9) )                 ! update nstop indicator (now sheared among all local domains)
-      ENDIF
+      ELSE
+         ! if no ocean point: MAXVAL returns -HUGE => we must overwrite this value to avoid error handling bellow.
+         IF( ll_0oce )   zmax(1:4) = (/ 0._wp, 0._wp, -1._wp, 1._wp /)   ! default "valid" values...
+      ENDIF
+      !
+      zmax(3) = -zmax(3)                         ! move back from max(-zz) to min(zz) : easier to manage! 
+      zmax(5) = -zmax(5)                         ! move back from max(-zz) to min(zz) : easier to manage!
+      IF( ll_colruns ) THEN
+         zmaxlocal(3) = -zmaxlocal(3)            ! move back from max(-zz) to min(zz) : easier to manage! 
+         zmaxlocal(5) = -zmaxlocal(5)            ! move back from max(-zz) to min(zz) : easier to manage!
+      ENDIF
+      !
       !                                   !==              write "run.stat" files              ==!
       !                                   !==  done only by 1st subdomain at writting timestep  ==!
       IF( ll_wrtruns ) THEN
-         WRITE(numrun,9500) kt, zmax(1), zmax(2), -zmax(3), zmax(4)
-         istatus = NF90_PUT_VAR( nrunid, nvarid(1), (/ zmax(1)/), (/kt/), (/1/) )
-         istatus = NF90_PUT_VAR( nrunid, nvarid(2), (/ zmax(2)/), (/kt/), (/1/) )
-         istatus = NF90_PUT_VAR( nrunid, nvarid(3), (/-zmax(3)/), (/kt/), (/1/) )
-         istatus = NF90_PUT_VAR( nrunid, nvarid(4), (/ zmax(4)/), (/kt/), (/1/) )
-         istatus = NF90_PUT_VAR( nrunid, nvarid(5), (/-zmax(5)/), (/kt/), (/1/) )
-         istatus = NF90_PUT_VAR( nrunid, nvarid(6), (/ zmax(6)/), (/kt/), (/1/) )
-         IF( ln_zad_Aimp ) THEN
-            istatus = NF90_PUT_VAR( nrunid, nvarid(7), (/ zmax(7)/), (/kt/), (/1/) )
-            istatus = NF90_PUT_VAR( nrunid, nvarid(8), (/ zmax(8)/), (/kt/), (/1/) )
-         ENDIF
+         WRITE(numrun,9500) kt, zmax(1), zmax(2), zmax(3), zmax(4)
+         DO ji = 1, 6 + 2 * COUNT( (/ln_zad_Aimp/) )
+            istatus = NF90_PUT_VAR( nrunid, nvarid(ji), (/zmax(ji)/), (/kt/), (/1/) )
+         END DO
          IF( kt == nitend )   istatus = NF90_CLOSE(nrunid)
       END IF
@@ -175,11 +183,11 @@
       !                                   !==  done by all processes at every time step  ==!
       !
-      IF(   zmax(1) >   20._wp .OR.   &                   ! too large sea surface height ( > 20 m )
-         &  zmax(2) >   10._wp .OR.   &                   ! too large velocity ( > 10 m/s)
-!!$         &  zmax(3) >=   0._wp .OR.   &                   ! negative or zero sea surface salinity
-!!$         &  zmax(4) >= 100._wp .OR.   &                   ! too large sea surface salinity ( > 100 )
-!!$         &  zmax(4) <    0._wp .OR.   &                   ! too large sea surface salinity (keep this line for sea-ice)
-         &  ISNAN( zmax(1) + zmax(2) + zmax(3) ) .OR.   &               ! NaN encounter in the tests
-         &  ABS(   zmax(1) + zmax(2) + zmax(3) ) > HUGE(1._wp) ) THEN   ! Infinity encounter in the tests
+      IF(  zmax(1) >   20._wp .OR.   &                   ! too large sea surface height ( > 20 m )
+         & zmax(2) >   10._wp .OR.   &                   ! too large velocity ( > 10 m/s)
+!!$         & zmax(3) <=   0._wp .OR.   &                   ! negative or zero sea surface salinity
+!!$         & zmax(4) >= 100._wp .OR.   &                   ! too large sea surface salinity ( > 100 )
+!!$         & zmax(4) <    0._wp .OR.   &                   ! too large sea surface salinity (keep this line for sea-ice)
+         & ISNAN( zmax(1) + zmax(2) + zmax(3) ) .OR.   &               ! NaN encounter in the tests
+         & ABS(   zmax(1) + zmax(2) + zmax(3) ) > HUGE(1._wp) ) THEN   ! Infinity encounter in the tests
          !
          iloc(:,:) = 0
@@ -207,5 +215,5 @@
          ELSE                    ! find local min and max locations:
             ! if we are here, this means that the subdomain contains some oce points -> no need to test the mask used in maxloc
-            llmsk(Nis0:Nie0,Njs0:Nje0,1) = ssmask(Nis0:Nie0,Njs0:Nje0 ) == 1._wp         ! define only the inner domain
+            llmsk(Nis0:Nie0,Njs0:Nje0,1) = ssmask(Nis0:Nie0,Njs0:Nje0 ) == 1._wp        ! define only the inner domain
             iloc(1:2,1) = MAXLOC( ABS( ssh(:,:,         Kmm)), mask = llmsk(:,:,1) )
             llmsk(Nis0:Nie0,Njs0:Nje0,:) = umask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp        ! define only the inner domain
@@ -221,8 +229,8 @@
          !
          WRITE(ctmp1,*) ' stp_ctl: |ssh| > 20 m  or  |U| > 10 m/s  or  S <= 0  or  S >= 100  or  NaN encounter in the tests'
-         CALL wrt_line( ctmp2, kt, '|ssh| max',  zmax(1), iloc(:,1), iareasum(1), iareamin(1), iareamax(1) )
-         CALL wrt_line( ctmp3, kt, '|U|   max',  zmax(2), iloc(:,2), iareasum(2), iareamin(2), iareamax(2) )
-         CALL wrt_line( ctmp4, kt, 'Sal   min', -zmax(3), iloc(:,3), iareasum(3), iareamin(3), iareamax(3) )
-         CALL wrt_line( ctmp5, kt, 'Sal   max',  zmax(4), iloc(:,4), iareasum(4), iareamin(4), iareamax(4) )
+         CALL wrt_line( ctmp2, kt, '|ssh| max', zmax(1), iloc(:,1), iareasum(1), iareamin(1), iareamax(1) )
+         CALL wrt_line( ctmp3, kt, '|U|   max', zmax(2), iloc(:,2), iareasum(2), iareamin(2), iareamax(2) )
+         CALL wrt_line( ctmp4, kt, 'Sal   min', zmax(3), iloc(:,3), iareasum(3), iareamin(3), iareamax(3) )
+         CALL wrt_line( ctmp5, kt, 'Sal   max', zmax(4), iloc(:,4), iareasum(4), iareamin(4), iareamax(4) )
          IF( Agrif_Root() ) THEN
             WRITE(ctmp6,*) '      ===> output of last computed fields in output.abort* files'
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/CPL_OASIS/EXPREF/file_def_nemo-ice.xml
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/CPL_OASIS/EXPREF/file_def_nemo-ice.xml	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/CPL_OASIS/EXPREF/file_def_nemo-ice.xml	(revision 13877)
@@ -53,7 +53,6 @@
        <field field_ref="normstr"          name="normstr" />
        <field field_ref="sheastr"          name="sheastr" />
-       <field field_ref="isig1"            name="isig1"   />
-       <field field_ref="isig2"            name="isig2"   />
-       <field field_ref="isig3"            name="isig3"   />
+       <field field_ref="sig1_pnorm"       name="sig1_pnorm"/>
+       <field field_ref="sig2_pnorm"       name="sig2_pnorm"/>
        
        <!-- heat fluxes -->
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ICE_ADV2D/EXPREF/file_def_nemo-ice.xml
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ICE_ADV2D/EXPREF/file_def_nemo-ice.xml	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ICE_ADV2D/EXPREF/file_def_nemo-ice.xml	(revision 13877)
@@ -55,7 +55,6 @@
         <field field_ref="normstr"          name="normstr" />
         <field field_ref="sheastr"          name="sheastr" />
-        <field field_ref="isig1"            name="isig1"   />
-        <field field_ref="isig2"            name="isig2"   />
-        <field field_ref="isig3"            name="isig3"   />
+       <field field_ref="sig1_pnorm"       name="sig1_pnorm"/>
+       <field field_ref="sig2_pnorm"       name="sig2_pnorm"/>
 
         <!-- heat fluxes -->
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ICE_AGRIF/EXPREF/file_def_nemo-ice.xml
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ICE_AGRIF/EXPREF/file_def_nemo-ice.xml	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ICE_AGRIF/EXPREF/file_def_nemo-ice.xml	(revision 13877)
@@ -53,7 +53,6 @@
         <field field_ref="normstr"          name="normstr" />
         <field field_ref="sheastr"          name="sheastr" />
-        <field field_ref="isig1"            name="isig1"   />
-        <field field_ref="isig2"            name="isig2"   />
-        <field field_ref="isig3"            name="isig3"   />
+       <field field_ref="sig1_pnorm"       name="sig1_pnorm"/>
+       <field field_ref="sig2_pnorm"       name="sig2_pnorm"/>
 
         <!-- heat fluxes -->
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/EXPREF/namelist_cfg
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/EXPREF/namelist_cfg	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/EXPREF/namelist_cfg	(revision 13877)
@@ -477,7 +477,4 @@
 /
 !-----------------------------------------------------------------------
-!-----------------------------------------------------------------------
-/
-!-----------------------------------------------------------------------
 &namhsb        !  Heat and salt budgets                                 (default: OFF)
 !-----------------------------------------------------------------------
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/MY_SRC/dtatsd.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/MY_SRC/dtatsd.F90	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/MY_SRC/dtatsd.F90	(revision 13877)
@@ -191,5 +191,5 @@
          ENDIF
          !
-         DO_2D( 1, 1, 1, 1 )
+         DO_2D( 1, 1, 1, 1 )                  ! vertical interpolation of T & S
             DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points
                zl = gdept_0(ji,jj,jk)
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/MY_SRC/eosbn2.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/MY_SRC/eosbn2.F90	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/MY_SRC/eosbn2.F90	(revision 13877)
@@ -182,4 +182,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -965,5 +966,5 @@
       IF( ln_timing )   CALL timing_start('bn2')
       !
-      DO_3D( 1, 1, 1, 1, 2, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 2, jpkm1 )      ! interior points only (2=< jk =< jpkm1 ); surface and bottom value set to zero one for all in istate.F90
          zrw =   ( gdepw(ji,jj,jk  ,Kmm) - gdept(ji,jj,jk,Kmm) )   &
             &  / ( gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm) ) 
@@ -1723,4 +1724,7 @@
          !
       CASE( np_leos )                        !==  Linear ISOMIP EOS     ==!
+
+         r1_S0  = 0.875_wp/35.16504_wp   ! Used to convert CT in potential temperature when using bulk formulae (eos_pt_from_ct)
+         
          IF(lwp) THEN
             WRITE(numout,*)
@@ -1731,4 +1735,5 @@
             WRITE(numout,*) '             saline  cont. coef.   rn_b0      = ', rn_b0
          ENDIF
+         l_useCT = .TRUE.          ! Use conservative temperature
          !
       CASE DEFAULT                     !==  ERROR in neos  ==!
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/MY_SRC/isf_oce.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/MY_SRC/isf_oce.F90	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/MY_SRC/isf_oce.F90	(revision 13877)
@@ -75,5 +75,5 @@
    !
    ! 2.1 -------- ice shelf cavity parameter --------------
-   LOGICAL , PUBLIC            :: l_isfoasis
+   LOGICAL , PUBLIC            :: l_isfoasis = .FALSE.
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)    ::   risfload                    !: ice shelf load
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)    ::   fwfisf_oasis
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/MY_SRC/isfcavgam.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/MY_SRC/isfcavgam.F90	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/MY_SRC/isfcavgam.F90	(revision 13877)
@@ -30,4 +30,5 @@
    PUBLIC   isfcav_gammats
 
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/MY_SRC/isfstp.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/MY_SRC/isfstp.F90	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/MY_SRC/isfstp.F90	(revision 13877)
@@ -13,5 +13,4 @@
    !!   isfstp       : compute iceshelf melt and heat flux
    !!----------------------------------------------------------------------
-   !
    USE isf_oce                                      ! isf variables
    USE isfload, ONLY: isf_load                      ! ice shelf load
@@ -21,5 +20,6 @@
    USE isfcpl , ONLY: isfcpl_rst_write, isfcpl_init ! isf variables
 
-   USE dom_oce, ONLY: ht, e3t, ln_isfcav, ln_linssh     ! ocean space and time domain
+   USE dom_oce        ! ocean space and time domain
+   USE oce      , ONLY: ssh                           ! sea surface height
    USE domvvl,  ONLY: ln_vvl_zstar                      ! zstar logical
    USE zdfdrg,  ONLY: r_Cdmin_top, r_ke0_top            ! vertical physics: top/bottom drag coef.
@@ -31,9 +31,10 @@
 
    IMPLICIT NONE
-
    PRIVATE
 
    PUBLIC   isf_stp, isf_init, isf_nam  ! routine called in sbcmod and divhor
 
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -41,4 +42,5 @@
    !! Software governed by the CeCILL license (see ./LICENSE)
    !!----------------------------------------------------------------------
+
 CONTAINS
  
@@ -60,4 +62,7 @@
       INTEGER, INTENT(in) ::   kt   ! ocean time step
       INTEGER, INTENT(in) ::   Kmm  ! ocean time level index
+      !!----------------------------------------------------------------------
+      INTEGER :: jk                               ! loop index
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t    ! e3t 
       !!---------------------------------------------------------------------
       !
@@ -78,5 +83,8 @@
          ! 1.2: compute misfkb, rhisf_tbl, rfrac (deepest level, thickness, fraction of deepest cell affected by tbl)
          rhisf_tbl_cav(:,:) = rn_htbl * mskisf_cav(:,:)
-         CALL isf_tbl_lvl(ht, e3t(:,:,:,Kmm), misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav)
+         DO jk = 1, jpk
+            ze3t(:,:,jk) = e3t(:,:,jk,Kmm)
+         END DO 
+         CALL isf_tbl_lvl(ht(:,:), ze3t, misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav)
          !
          ! 1.3: compute ice shelf melt
@@ -100,5 +108,8 @@
          ! by simplicity, we assume the top level where param applied do not change with time (done in init part)
          rhisf_tbl_par(:,:) = rhisf0_tbl_par(:,:)
-         CALL isf_tbl_lvl(ht, e3t(:,:,:,Kmm), misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par)
+         DO jk = 1, jpk
+            ze3t(:,:,jk) = e3t(:,:,jk,Kmm)
+         END DO
+         CALL isf_tbl_lvl(ht(:,:), ze3t, misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par)
          !
          ! 2.3: compute ice shelf melt
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/MY_SRC/istate.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/MY_SRC/istate.F90	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/MY_SRC/istate.F90	(revision 13877)
@@ -24,5 +24,4 @@
    USE dom_oce        ! ocean space and time domain 
    USE daymod         ! calendar
-   USE divhor         ! horizontal divergence            (div_hor routine)
    USE dtatsd         ! data temperature and salinity   (dta_tsd routine)
    USE dtauvd         ! data: U & V current             (dta_uvd routine)
@@ -35,4 +34,8 @@
    USE lib_mpp         ! MPP library
    USE restart         ! restart
+#if defined key_agrif
+   USE agrif_oce_interp
+   USE agrif_oce
+#endif   
 
    IMPLICIT NONE
@@ -43,4 +46,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -59,4 +63,5 @@
       !
       INTEGER ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zgdept     ! 3D table  !!st patch to use gdept subtitute
 !!gm see comment further down
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   zuvd    ! U & V data workspace
@@ -70,5 +75,5 @@
 !!gm  Why not include in the first call of dta_tsd ?  
 !!gm  probably associated with the use of internal damping...
-                     CALL dta_tsd_init        ! Initialisation of T & S input data
+       CALL dta_tsd_init        ! Initialisation of T & S input data
 !!gm to be moved in usrdef of C1D case
 !      IF( lk_c1d )   CALL dta_uvd_init        ! Initialization of U & V input data
@@ -84,4 +89,17 @@
 #endif
 
+#if defined key_agrif
+      IF ( (.NOT.Agrif_root()).AND.ln_init_chfrpar ) THEN
+         numror = 0                           ! define numror = 0 -> no restart file to read
+         ln_1st_euler = .true.                ! Set time-step indicator at nit000 (euler forward)
+         CALL day_init 
+         CALL agrif_istate( Kbb, Kmm, Kaa )   ! Interp from parent
+         !
+         ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb) 
+         ssh (:,:,Kmm)     = ssh(:,:,Kbb)
+         uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb)
+         vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb)
+      ELSE
+#endif
       IF( ln_rstart ) THEN                    ! Restart from a file
          !                                    ! -------------------
@@ -100,4 +118,7 @@
             !
             ssh(:,:,Kbb)   = 0._wp               ! set the ocean at rest
+            uu  (:,:,:,Kbb) = 0._wp
+            vv  (:,:,:,Kbb) = 0._wp  
+            !
             IF( ll_wd ) THEN
                ssh(:,:,Kbb) =  -ssh_ref  ! Added in 30 here for bathy that adds 30 as Iterative test CEOD 
@@ -111,9 +132,10 @@
                END_2D
             ENDIF 
-            uu  (:,:,:,Kbb) = 0._wp
-            vv  (:,:,:,Kbb) = 0._wp  
-            !
+             !
          ELSE                                 ! user defined initial T and S
-            CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )         
+            DO jk = 1, jpk
+               zgdept(:,:,jk) = gdept(:,:,jk,Kbb)
+            END DO
+            CALL usr_def_istate( zgdept, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )         
          ENDIF
          ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones
@@ -121,11 +143,8 @@
          uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb)
          vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb)
-         hdiv(:,:,jpk) = 0._wp               ! bottom divergence set one for 0 to zero at jpk level
-         CALL div_hor( 0, Kbb, Kmm )         ! compute interior hdiv value  
-!!gm                                    hdiv(:,:,:) = 0._wp
 
 !!gm POTENTIAL BUG :
 !!gm  ISSUE :  if ssh(:,:,Kbb) /= 0  then, in non linear free surface, the e3._n, e3._b should be recomputed
-!!             as well as gdept and gdepw....   !!!!! 
+!!             as well as gdept_ and gdepw_....   !!!!! 
 !!      ===>>>>   probably a call to domvvl initialisation here....
 
@@ -151,4 +170,7 @@
          ! 
       ENDIF 
+#if defined key_agrif
+      ENDIF
+#endif
       ! 
       ! Initialize "now" and "before" barotropic velocities:
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/MY_SRC/sbcfwb.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/MY_SRC/sbcfwb.F90	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP+/MY_SRC/sbcfwb.F90	(revision 13877)
@@ -17,5 +17,5 @@
    USE dom_oce        ! ocean space and time domain
    USE sbc_oce        ! surface ocean boundary condition
-   USE isf_oce        ! ice shelf melting contribution
+   USE isf_oce , ONLY : fwfisf_cav, fwfisf_par, ln_isfcpl, ln_isfcpl_cons, risfcpl_cons_ssh ! ice shelf melting contribution
    USE sbc_ice , ONLY : snwice_mass, snwice_mass_b, snwice_fmass
    USE phycst         ! physical constants
@@ -71,5 +71,5 @@
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   ztmsk_tospread, zerp_cor    !   -      -
       REAL(wp)   ,DIMENSION(1) ::   z_fwfprv  
-      COMPLEX(wp),DIMENSION(1) ::   y_fwfnow  
+      COMPLEX(dp),DIMENSION(1) ::   y_fwfnow  
       !!----------------------------------------------------------------------
       !
@@ -207,7 +207,7 @@
 !!gm   ===>>>>  lbc_lnk should be useless as all the computation is done over the whole domain !
 #if defined key_mpi3
-            CALL lbc_lnk_nc_multi( 'sbcfwb', zerp_cor, 'T', 1. )
+            CALL lbc_lnk_nc_lnk( 'sbcfwb', zerp_cor, 'T', 1.0_wp )
 #else
-            CALL lbc_lnk( 'sbcfwb', zerp_cor, 'T', 1. )
+            CALL lbc_lnk( 'sbcfwb', zerp_cor, 'T', 1.0_wp )
 #endif
             !
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP/EXPREF/namelist_cfg
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP/EXPREF/namelist_cfg	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/ISOMIP/EXPREF/namelist_cfg	(revision 13877)
@@ -437,7 +437,4 @@
 /
 !-----------------------------------------------------------------------
-!-----------------------------------------------------------------------
-/
-!-----------------------------------------------------------------------
 &namhsb        !  Heat and salt budgets                                 (default: OFF)
 !-----------------------------------------------------------------------
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/STATION_ASF/MY_SRC/stpctl.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/STATION_ASF/MY_SRC/stpctl.F90	(revision 13807)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/tests/STATION_ASF/MY_SRC/stpctl.F90	(revision 13877)
@@ -65,5 +65,5 @@
       REAL(wp)                        ::   zzz                                   ! local real 
       REAL(wp), DIMENSION(4)          ::   zmax, zmaxlocal
-      LOGICAL                         ::   ll_wrtstp, ll_colruns, ll_wrtruns
+      LOGICAL                         ::   ll_wrtstp, ll_colruns, ll_wrtruns, ll_0oce
       LOGICAL, DIMENSION(jpi,jpj)     ::   llmsk
       CHARACTER(len=20)               ::   clname
@@ -115,17 +115,13 @@
       llmsk(:,Nje1: jpj) = .FALSE.
       !
-      llmsk(Nis0:Nie0,Njs0:Nje0) = tmask(Nis0:Nie0,Njs0:Nje0,1) == 1._wp   ! test only the inner domain
-      IF( COUNT( llmsk(:,:) ) > 0 ) THEN   ! avoid huge values sent back for land processors...
-         zmax(1) = MAXVAL(     taum(:,:)   , mask = llmsk )   ! max wind stress module
-         zmax(2) = MAXVAL( ABS( qns(:,:) ) , mask = llmsk )   ! max non-solar heat flux
-         zmax(3) = MAXVAL( ABS( emp(:,:) ) , mask = llmsk )   ! max E-P
-      ELSE
-         IF( ll_colruns ) THEN    ! default value: must not be kept when calling mpp_max -> must be as small as possible
-            zmax(1:3) = -HUGE(1._wp)
-         ELSE                     ! default value: must not give true for any of the tests bellow (-> avoid manipulating HUGE...)
-            zmax(1:3) = 0._wp
-         ENDIF
-      ENDIF
-      zmax(4) = REAL( nstop, wp )                                     ! stop indicator
+      llmsk(Nis0:Nie0,Njs0:Nje0) = tmask(Nis0:Nie0,Njs0:Nje0,1) == 1._wp        ! test only the inner domain
+      !
+      ll_0oce = .NOT. ANY( llmsk(:,:) )                                         ! no ocean point in the inner domain?
+      !
+      zmax(1) = MAXVAL(     taum(:,:)  , mask = llmsk )                         ! max wind stress module
+      zmax(2) = MAXVAL( ABS( qns(:,:) ), mask = llmsk )                         ! max non-solar heat flux
+      zmax(3) = MAXVAL( ABS( emp(:,:) ), mask = llmsk )                         ! max E-P
+      zmax(4) = REAL( nstop, wp )                                               ! stop indicator
+      !
       !                                   !==               get global extrema             ==!
       !                                   !==  done by all processes if writting run.stat  ==!
@@ -134,12 +130,16 @@
          CALL mpp_max( "stpctl", zmax )          ! max over the global domain
          nstop = NINT( zmax(4) )                 ! update nstop indicator (now sheared among all local domains)
-      ENDIF
+      ELSE
+         ! if no ocean point: MAXVAL returns -HUGE => we must overwrite this value to avoid error handling bellow.
+         IF( ll_0oce )   zmax(1:3) = 0._wp       ! default "valid" values...
+      ENDIF
+      !                                   !==               error handling               ==!
       !                                   !==              write "run.stat" files              ==!
       !                                   !==  done only by 1st subdomain at writting timestep  ==!
       IF( ll_wrtruns ) THEN
          WRITE(numrun,9500) kt, zmax(1), zmax(2), zmax(3)
-         istatus = NF90_PUT_VAR( nrunid, nvarid(1), (/ zmax(1)/), (/kt/), (/1/) )
-         istatus = NF90_PUT_VAR( nrunid, nvarid(2), (/ zmax(2)/), (/kt/), (/1/) )
-         istatus = NF90_PUT_VAR( nrunid, nvarid(3), (/ zmax(3)/), (/kt/), (/1/) )
+         DO ji = 1, 3
+            istatus = NF90_PUT_VAR( nrunid, nvarid(ji), (/zmax(ji)/), (/kt/), (/1/) )
+         END DO
          IF( kt == nitend ) istatus = NF90_CLOSE(nrunid)
       END IF
