Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ABL/ablmod.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ABL/ablmod.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ABL/ablmod.F90	(revision 13193)
@@ -592,8 +592,8 @@
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !                            !  8 *** Swap time indices for the next timestep
-      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
-      nt_n = 1 + MOD( kt  , 2)
-      nt_a = 1 + MOD( kt+1, 2)
-      !	  
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      nt_n = 1 + MOD( nt_n, 2)
+      nt_a = 1 + MOD( nt_a, 2)
+      !
 !---------------------------------------------------------------------------------------------------
    END SUBROUTINE abl_stp
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ABL/par_abl.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ABL/par_abl.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ABL/par_abl.F90	(revision 13193)
@@ -29,4 +29,5 @@
    LOGICAL , PUBLIC            ::   ln_smth_pblh   !: smoothing of atmospheric PBL height 
 
+   LOGICAL           , PUBLIC ::   ln_rstart_abl    !: (de)activate abl restart
    CHARACTER(len=256), PUBLIC ::   cn_ablrst_in     !: suffix of abl restart name (input)
    CHARACTER(len=256), PUBLIC ::   cn_ablrst_out    !: suffix of abl restart name (output)
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ABL/sbcabl.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ABL/sbcabl.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ABL/sbcabl.F90	(revision 13193)
@@ -68,5 +68,5 @@
       LOGICAL            ::   lluldl
       NAMELIST/namsbc_abl/ cn_dir, cn_dom, cn_ablrst_in, cn_ablrst_out,           &
-         &                 cn_ablrst_indir, cn_ablrst_outdir,                     &
+         &                 cn_ablrst_indir, cn_ablrst_outdir, ln_rstart_abl,      &
          &                 ln_hpgls_frc, ln_geos_winds, nn_dyn_restore,           &
          &                 rn_ldyn_min , rn_ldyn_max, rn_ltra_min, rn_ltra_max,   &
@@ -263,9 +263,8 @@
 
       ! Initialize the time index for now time (nt_n) and after time (nt_a)
-      nt_n = 1 + MOD( nit000  , 2)
-      nt_a = 1 + MOD( nit000+1, 2)
+      nt_n = 1; nt_a = 2
 
       ! initialize ABL from data or restart
-      IF( ln_rstart ) THEN
+      IF( ln_rstart_abl ) THEN
          CALL abl_rst_read
       ELSE
@@ -288,6 +287,4 @@
       ENDIF
 
-      rhoa(:,:) = rho_air( tq_abl(:,:,2,nt_n,jp_ta), tq_abl(:,:,2,nt_n,jp_qa), sf(jp_slp)%fnow(:,:,1) ) !!GS: rhoa must be (re)computed here here to avoid division by zero in blk_ice_1 (TBI)
-
    END SUBROUTINE sbc_abl_init
 
@@ -329,55 +326,58 @@
       CALL fld_read( kt, nn_fsbc, sf )             ! input fields provided at the current time-step
 
-      !!-------------------------------------------------------------------------------------------
-      !! 2 - Compute Cd x ||U||, Ch x ||U||, Ce x ||U||, and SSQ using now fields
-      !!-------------------------------------------------------------------------------------------
-
-      CALL blk_oce_1( kt,  u_abl(:,:,2,nt_n      ),  v_abl(:,:,2,nt_n      ),   &   !   <<= in
-         &                tq_abl(:,:,2,nt_n,jp_ta), tq_abl(:,:,2,nt_n,jp_qa),   &   !   <<= in
-         &                sf(jp_slp )%fnow(:,:,1) , sst_m, ssu_m, ssv_m     ,   &   !   <<= in
-         &                sf(jp_qsr )%fnow(:,:,1) , sf(jp_qlw )%fnow(:,:,1) ,   &   !   <<= in
-         &                tsk_m, zssq, zcd_du, zsen, zevp                       )   !   =>> out
-
-#if defined key_si3
-      CALL blk_ice_1(  u_abl(:,:,2,nt_n      ),  v_abl(:,:,2,nt_n      ),    &   !   <<= in
-         &            tq_abl(:,:,2,nt_n,jp_ta), tq_abl(:,:,2,nt_n,jp_qa),    &   !   <<= in
-         &            sf(jp_slp)%fnow(:,:,1)  ,  u_ice, v_ice, tm_su    ,    &   !   <<= in
-         &            pseni=zseni, pevpi=zevpi, pssqi=zssqi, pcd_dui=zcd_dui )   !   <<= out
-#endif
-
-      !!-------------------------------------------------------------------------------------------
-      !! 3 - Advance ABL variables from now (n) to after (n+1)
-      !!-------------------------------------------------------------------------------------------
-
-      CALL abl_stp( kt, tsk_m, ssu_m, ssv_m, zssq,                          &   !   <<= in
-         &              sf(jp_wndi)%fnow(:,:,:), sf(jp_wndj)%fnow(:,:,:),   &   !   <<= in
-         &              sf(jp_tair)%fnow(:,:,:), sf(jp_humi)%fnow(:,:,:),   &   !   <<= in
-         &              sf(jp_slp )%fnow(:,:,1),                            &   !   <<= in
-         &              sf(jp_hpgi)%fnow(:,:,:), sf(jp_hpgj)%fnow(:,:,:),   &   !   <<= in
-         &              zcd_du, zsen, zevp,                                 &   !   <=> in/out
-         &              wndm, utau, vtau, taum                              &   !   =>> out
-#if defined key_si3
-         &            , tm_su, u_ice, v_ice, zssqi, zcd_dui                 &   !   <<= in
-         &            , zseni, zevpi, wndm_ice, ato_i                       &   !   <<= in
-         &            , utau_ice, vtau_ice                                  &   !   =>> out
-#endif
-         &                                                                  )
-      !!-------------------------------------------------------------------------------------------
-      !! 4 - Finalize flux computation using ABL variables at (n+1), nt_n corresponds to (n+1) since
-      !!                                                                time swap is done in abl_stp
-      !!-------------------------------------------------------------------------------------------
-
-      CALL blk_oce_2( tq_abl(:,:,2,nt_n,jp_ta),                            &
-         &            sf(jp_qsr )%fnow(:,:,1) , sf(jp_qlw )%fnow(:,:,1),   &
-         &            sf(jp_prec)%fnow(:,:,1) , sf(jp_snow)%fnow(:,:,1),   &
-         &            tsk_m, zsen, zevp                                )
-
-      CALL abl_rst_opn( kt )                       ! Open abl restart file (if necessary)
-      IF( lrst_abl ) CALL abl_rst_write( kt )      ! -- abl restart file
-
-#if defined key_si3
-      ! Avoid a USE abl in icesbc module
-      sf(jp_tair)%fnow(:,:,1) = tq_abl(:,:,2,nt_n,jp_ta);  sf(jp_humi)%fnow(:,:,1) = tq_abl(:,:,2,nt_n,jp_qa)
-#endif
+      IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
+
+         !!-------------------------------------------------------------------------------------------
+         !! 2 - Compute Cd x ||U||, Ch x ||U||, Ce x ||U||, and SSQ using now fields
+         !!-------------------------------------------------------------------------------------------
+
+         CALL blk_oce_1( kt,  u_abl(:,:,2,nt_n      ),  v_abl(:,:,2,nt_n      ),   &   !   <<= in
+            &                tq_abl(:,:,2,nt_n,jp_ta), tq_abl(:,:,2,nt_n,jp_qa),   &   !   <<= in
+            &                sf(jp_slp )%fnow(:,:,1) , sst_m, ssu_m, ssv_m     ,   &   !   <<= in
+            &                sf(jp_qsr )%fnow(:,:,1) , sf(jp_qlw )%fnow(:,:,1) ,   &   !   <<= in
+            &                tsk_m, zssq, zcd_du, zsen, zevp                       )   !   =>> out
+
+#if defined key_si3
+         CALL blk_ice_1(  u_abl(:,:,2,nt_n      ),  v_abl(:,:,2,nt_n      ),    &   !   <<= in
+            &            tq_abl(:,:,2,nt_n,jp_ta), tq_abl(:,:,2,nt_n,jp_qa),    &   !   <<= in
+            &            sf(jp_slp)%fnow(:,:,1)  ,  u_ice, v_ice, tm_su    ,    &   !   <<= in
+            &            pseni=zseni, pevpi=zevpi, pssqi=zssqi, pcd_dui=zcd_dui )   !   <<= out
+#endif
+
+         !!-------------------------------------------------------------------------------------------
+         !! 3 - Advance ABL variables from now (n) to after (n+1)
+         !!-------------------------------------------------------------------------------------------
+   
+         CALL abl_stp( kt, tsk_m, ssu_m, ssv_m, zssq,                          &   !   <<= in
+            &              sf(jp_wndi)%fnow(:,:,:), sf(jp_wndj)%fnow(:,:,:),   &   !   <<= in
+            &              sf(jp_tair)%fnow(:,:,:), sf(jp_humi)%fnow(:,:,:),   &   !   <<= in
+            &              sf(jp_slp )%fnow(:,:,1),                            &   !   <<= in
+            &              sf(jp_hpgi)%fnow(:,:,:), sf(jp_hpgj)%fnow(:,:,:),   &   !   <<= in
+            &              zcd_du, zsen, zevp,                                 &   !   <=> in/out
+            &              wndm, utau, vtau, taum                              &   !   =>> out
+#if defined key_si3
+            &            , tm_su, u_ice, v_ice, zssqi, zcd_dui                 &   !   <<= in
+            &            , zseni, zevpi, wndm_ice, ato_i                       &   !   <<= in
+            &            , utau_ice, vtau_ice                                  &   !   =>> out
+#endif
+            &                                                                  )
+         !!-------------------------------------------------------------------------------------------
+         !! 4 - Finalize flux computation using ABL variables at (n+1), nt_n corresponds to (n+1) since
+         !!                                                                time swap is done in abl_stp
+         !!-------------------------------------------------------------------------------------------
+
+         CALL blk_oce_2( tq_abl(:,:,2,nt_n,jp_ta),                            &
+            &            sf(jp_qsr )%fnow(:,:,1) , sf(jp_qlw )%fnow(:,:,1),   &
+            &            sf(jp_prec)%fnow(:,:,1) , sf(jp_snow)%fnow(:,:,1),   &
+            &            tsk_m, zsen, zevp                                )
+   
+         CALL abl_rst_opn( kt )                       ! Open abl restart file (if necessary)
+         IF( lrst_abl ) CALL abl_rst_write( kt )      ! -- abl restart file
+
+#if defined key_si3
+         ! Avoid a USE abl in icesbc module
+         sf(jp_tair)%fnow(:,:,1) = tq_abl(:,:,2,nt_n,jp_ta);  sf(jp_humi)%fnow(:,:,1) = tq_abl(:,:,2,nt_n,jp_qa)
+#endif
+      END IF
 
    END SUBROUTINE sbc_abl
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/iceistate.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/iceistate.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/iceistate.F90	(revision 13193)
@@ -179,7 +179,7 @@
             !
             ! -- mandatory fields -- !
-            zht_i_ini(:,:) = si(jp_hti)%fnow(:,:,1)
-            zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1)
-            zat_i_ini(:,:) = si(jp_ati)%fnow(:,:,1)
+            zht_i_ini(:,:) = si(jp_hti)%fnow(:,:,1) * tmask(:,:,1)
+            zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1) * tmask(:,:,1)
+            zat_i_ini(:,:) = si(jp_ati)%fnow(:,:,1) * tmask(:,:,1)
 
             ! -- optional fields -- !
@@ -219,10 +219,10 @@
                &     si(jp_hpd)%fnow(:,:,1) = ( rn_hpd_ini_n * zswitch + rn_hpd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
             !
-            zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1)
-            ztm_i_ini(:,:) = si(jp_tmi)%fnow(:,:,1)
-            zt_su_ini(:,:) = si(jp_tsu)%fnow(:,:,1)
-            ztm_s_ini(:,:) = si(jp_tms)%fnow(:,:,1)
-            zapnd_ini(:,:) = si(jp_apd)%fnow(:,:,1)
-            zhpnd_ini(:,:) = si(jp_hpd)%fnow(:,:,1)
+            zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1) * tmask(:,:,1)
+            ztm_i_ini(:,:) = si(jp_tmi)%fnow(:,:,1) * tmask(:,:,1)
+            zt_su_ini(:,:) = si(jp_tsu)%fnow(:,:,1) * tmask(:,:,1)
+            ztm_s_ini(:,:) = si(jp_tms)%fnow(:,:,1) * tmask(:,:,1)
+            zapnd_ini(:,:) = si(jp_apd)%fnow(:,:,1) * tmask(:,:,1)
+            zhpnd_ini(:,:) = si(jp_hpd)%fnow(:,:,1) * tmask(:,:,1)
             !
             ! change the switch for the following
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/BDY/bdydta.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/BDY/bdydta.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/BDY/bdydta.F90	(revision 13193)
@@ -92,5 +92,4 @@
       INTEGER ::  jbdy, jfld, jstart, jend, ib, jl    ! dummy loop indices
       INTEGER ::  ii, ij, ik, igrd, ipl               ! local integers
-      INTEGER,   DIMENSION(jpbgrd)     ::   ilen1 
       TYPE(OBC_DATA)         , POINTER ::   dta_alias        ! short cut
       TYPE(FLD), DIMENSION(:), POINTER ::   bf_alias
@@ -117,13 +116,15 @@
                   END DO
                ENDIF
-               IF( dta_bdy(jbdy)%lneed_dyn2d .AND. ASSOCIATED(dta_bdy(jbdy)%u2d) ) THEN   ! no SIZE with a unassociated pointer
+               IF( ASSOCIATED(dta_bdy(jbdy)%u2d) ) THEN   ! no SIZE with a unassociated pointer. v2d and u2d can differ on subdomain
                   igrd = 2
-                  DO ib = 1, SIZE(dta_bdy(jbdy)%u2d)   ! u2d is used only on the rim except if ln_full_vel = T, see bdy_dta_init
+                  DO ib = 1, SIZE(dta_bdy(jbdy)%u2d)      ! u2d is used either over the whole bdy or only on the rim
                      ii = idx_bdy(jbdy)%nbi(ib,igrd)
                      ij = idx_bdy(jbdy)%nbj(ib,igrd)
                      dta_bdy(jbdy)%u2d(ib) = uu_b(ii,ij,Kmm) * umask(ii,ij,1)         
                   END DO
+               ENDIF
+               IF( ASSOCIATED(dta_bdy(jbdy)%v2d) ) THEN   ! no SIZE with a unassociated pointer. v2d and u2d can differ on subdomain
                   igrd = 3
-                  DO ib = 1, SIZE(dta_bdy(jbdy)%v2d)   ! v2d is used only on the rim except if ln_full_vel = T, see bdy_dta_init
+                  DO ib = 1, SIZE(dta_bdy(jbdy)%v2d)      ! v2d is used either over the whole bdy or only on the rim
                      ii = idx_bdy(jbdy)%nbi(ib,igrd)
                      ij = idx_bdy(jbdy)%nbj(ib,igrd)
@@ -211,25 +212,29 @@
          !
          ! if runoff condition: change river flow we read (in m3/s) into barotropic velocity (m/s)
-         IF( cn_tra(jbdy) == 'runoff' .AND. TRIM(bf_alias(jp_bdyu2d)%clrootname) /= 'NOT USED' ) THEN   ! runoff and we read u/v2d
+         IF( cn_tra(jbdy) == 'runoff' ) THEN   ! runoff
             !
-            igrd = 2                      ! zonal flow (m3/s) to barotropic zonal velocity (m/s)
-            DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
-               ii   = idx_bdy(jbdy)%nbi(ib,igrd)
-               ij   = idx_bdy(jbdy)%nbj(ib,igrd)
-               dta_alias%u2d(ib) = dta_alias%u2d(ib) / ( e2u(ii,ij) * hu_0(ii,ij) )
-            END DO
-            igrd = 3                      ! meridional flow (m3/s) to barotropic meridional velocity (m/s)
-            DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
-               ii   = idx_bdy(jbdy)%nbi(ib,igrd)
-               ij   = idx_bdy(jbdy)%nbj(ib,igrd)
-               dta_alias%v2d(ib) = dta_alias%v2d(ib) / ( e1v(ii,ij) * hv_0(ii,ij) )
-            END DO
+            IF( ASSOCIATED(dta_bdy(jbdy)%u2d) ) THEN   ! no SIZE with a unassociated pointer. v2d and u2d can differ on subdomain
+               igrd = 2                         ! zonal flow (m3/s) to barotropic zonal velocity (m/s)
+               DO ib = 1, SIZE(dta_alias%u2d)   ! u2d is used either over the whole bdy or only on the rim
+                  ii   = idx_bdy(jbdy)%nbi(ib,igrd)
+                  ij   = idx_bdy(jbdy)%nbj(ib,igrd)
+                  dta_alias%u2d(ib) = dta_alias%u2d(ib) / ( e2u(ii,ij) * hu_0(ii,ij) )
+               END DO
+            ENDIF
+            IF( ASSOCIATED(dta_bdy(jbdy)%v2d) ) THEN   ! no SIZE with a unassociated pointer. v2d and u2d can differ on subdomain
+               igrd = 3                         ! meridional flow (m3/s) to barotropic meridional velocity (m/s)
+               DO ib = 1, SIZE(dta_alias%v2d)   ! v2d is used either over the whole bdy or only on the rim
+                  ii   = idx_bdy(jbdy)%nbi(ib,igrd)
+                  ij   = idx_bdy(jbdy)%nbj(ib,igrd)
+                  dta_alias%v2d(ib) = dta_alias%v2d(ib) / ( e1v(ii,ij) * hv_0(ii,ij) )
+               END DO
+            ENDIF
          ENDIF
 
          ! tidal harmonic forcing ONLY: initialise arrays
          IF( nn_dyn2d_dta(jbdy) == 2 ) THEN   ! we did not read ssh, u/v2d 
-            IF( dta_alias%lneed_ssh   .AND. ASSOCIATED(dta_alias%ssh) ) dta_alias%ssh(:) = 0._wp
-            IF( dta_alias%lneed_dyn2d .AND. ASSOCIATED(dta_alias%u2d) ) dta_alias%u2d(:) = 0._wp
-            IF( dta_alias%lneed_dyn2d .AND. ASSOCIATED(dta_alias%v2d) ) dta_alias%v2d(:) = 0._wp
+            IF( ASSOCIATED(dta_alias%ssh) ) dta_alias%ssh(:) = 0._wp
+            IF( ASSOCIATED(dta_alias%u2d) ) dta_alias%u2d(:) = 0._wp
+            IF( ASSOCIATED(dta_alias%v2d) ) dta_alias%v2d(:) = 0._wp
          ENDIF
 
@@ -333,15 +338,11 @@
             DO jbdy = 1, nb_bdy      ! Tidal component added in ts loop
                IF ( nn_dyn2d_dta(jbdy) .GE. 2 ) THEN
-                  IF( cn_dyn2d(jbdy) == 'frs' ) THEN   ;   ilen1(:)=idx_bdy(jbdy)%nblen(:)
-                  ELSE                                 ;   ilen1(:)=idx_bdy(jbdy)%nblenrim(:)
-                  ENDIF
-                  IF ( dta_bdy(jbdy)%lneed_ssh   ) dta_bdy_s(jbdy)%ssh(1:ilen1(1)) = dta_bdy(jbdy)%ssh(1:ilen1(1))
-                  IF ( dta_bdy(jbdy)%lneed_dyn2d ) dta_bdy_s(jbdy)%u2d(1:ilen1(2)) = dta_bdy(jbdy)%u2d(1:ilen1(2))
-                  IF ( dta_bdy(jbdy)%lneed_dyn2d ) dta_bdy_s(jbdy)%v2d(1:ilen1(3)) = dta_bdy(jbdy)%v2d(1:ilen1(3))
+                  IF( ASSOCIATED(dta_bdy(jbdy)%ssh) ) dta_bdy_s(jbdy)%ssh(:) = dta_bdy(jbdy)%ssh(:)
+                  IF( ASSOCIATED(dta_bdy(jbdy)%u2d) ) dta_bdy_s(jbdy)%u2d(:) = dta_bdy(jbdy)%u2d(:)
+                  IF( ASSOCIATED(dta_bdy(jbdy)%v2d) ) dta_bdy_s(jbdy)%v2d(:) = dta_bdy(jbdy)%v2d(:)
                ENDIF
             END DO
          ELSE ! Add tides if not split-explicit free surface else this is done in ts loop
             !
-            ! BDY: use pt_offset=1.0 as applied at the end of the step and bdy_dta_tides is referenced at the middle of the step
             CALL bdy_dta_tides( kt=kt, pt_offset = 1._wp )
          ENDIF
@@ -351,5 +352,5 @@
       !
    END SUBROUTINE bdy_dta
-
+   
 
    SUBROUTINE bdy_dta_init
@@ -383,4 +384,5 @@
       LOGICAL                                ::   llneed        !
       LOGICAL                                ::   llread        !
+      LOGICAL                                ::   llfullbdy     !
       TYPE(FLD_N), DIMENSION(1), TARGET  ::   bn_tem, bn_sal, bn_u3d, bn_v3d   ! must be an array to be used with fld_fill
       TYPE(FLD_N), DIMENSION(1), TARGET  ::   bn_ssh, bn_u2d, bn_v2d           ! informations about the fields to be read
@@ -497,10 +499,11 @@
                igrd = 2                                                    ! U point
                ipk = 1                                                     ! surface data
-               llneed = dta_bdy(jbdy)%lneed_dyn2d                          ! dta_bdy(jbdy)%ssh will be needed
+               llneed = dta_bdy(jbdy)%lneed_dyn2d                          ! dta_bdy(jbdy)%u2d will be needed
                llread = .NOT. ln_full_vel .AND. MOD(nn_dyn2d_dta(jbdy),2) == 1   ! don't get u2d from u3d and read NetCDF file
                bf_alias => bf(jp_bdyu2d,jbdy:jbdy)                         ! alias for u2d structure of bdy number jbdy
                bn_alias => bn_u2d                                          ! alias for u2d structure of nambdy_dta
-               IF( ln_full_vel ) THEN  ;   iszdim = idx_bdy(jbdy)%nblen(igrd)      ! will be computed from u3d -> need on the full bdy
-               ELSE                    ;   iszdim = idx_bdy(jbdy)%nblenrim(igrd)   ! used only on the rim
+               llfullbdy = ln_full_vel .OR. cn_dyn2d(jbdy) == 'frs'        ! need u2d over the whole bdy or only over the rim?
+               IF( llfullbdy ) THEN  ;   iszdim = idx_bdy(jbdy)%nblen(igrd)
+               ELSE                  ;   iszdim = idx_bdy(jbdy)%nblenrim(igrd)
                ENDIF
             ENDIF
@@ -509,10 +512,11 @@
                igrd = 3                                                    ! V point
                ipk = 1                                                     ! surface data
-               llneed = dta_bdy(jbdy)%lneed_dyn2d                          ! dta_bdy(jbdy)%ssh will be needed
+               llneed = dta_bdy(jbdy)%lneed_dyn2d                          ! dta_bdy(jbdy)%v2d will be needed
                llread = .NOT. ln_full_vel .AND. MOD(nn_dyn2d_dta(jbdy),2) == 1   ! don't get v2d from v3d and read NetCDF file
                bf_alias => bf(jp_bdyv2d,jbdy:jbdy)                         ! alias for v2d structure of bdy number jbdy
                bn_alias => bn_v2d                                          ! alias for v2d structure of nambdy_dta 
-               IF( ln_full_vel ) THEN  ;   iszdim = idx_bdy(jbdy)%nblen(igrd)      ! will be computed from v3d -> need on the full bdy
-               ELSE                    ;   iszdim = idx_bdy(jbdy)%nblenrim(igrd)   ! used only on the rim
+               llfullbdy = ln_full_vel .OR. cn_dyn2d(jbdy) == 'frs'        ! need v2d over the whole bdy or only over the rim?
+               IF( llfullbdy ) THEN  ;   iszdim = idx_bdy(jbdy)%nblen(igrd)
+               ELSE                  ;   iszdim = idx_bdy(jbdy)%nblenrim(igrd)
                ENDIF
             ENDIF
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/BDY/bdyini.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/BDY/bdyini.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/BDY/bdyini.F90	(revision 13193)
@@ -19,9 +19,10 @@
    USE oce            ! ocean dynamics and tracers variables
    USE dom_oce        ! ocean space and time domain
+   USE sbc_oce , ONLY: nn_ice
    USE bdy_oce        ! unstructured open boundary conditions
    USE bdydta         ! open boundary cond. setting   (bdy_dta_init routine)
    USE bdytides       ! open boundary cond. setting   (bdytide_init routine)
    USE tide_mod, ONLY: ln_tide ! tidal forcing
-   USE phycst   , ONLY: rday
+   USE phycst  , ONLY: rday
    !
    USE in_out_manager ! I/O units
@@ -315,4 +316,9 @@
 
          dta_bdy(ib_bdy)%lneed_ice = cn_ice(ib_bdy) /= 'none'
+
+         IF( dta_bdy(ib_bdy)%lneed_ice .AND. nn_ice /= 2 ) THEN
+            WRITE(ctmp1,*) 'bdy number ', ib_bdy,', needs ice model but nn_ice = ', nn_ice
+            CALL ctl_stop( ctmp1 )
+         ENDIF
 
          IF( lwp .AND. dta_bdy(ib_bdy)%lneed_ice ) THEN 
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/BDY/bdytides.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/BDY/bdytides.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/BDY/bdytides.F90	(revision 13193)
@@ -65,19 +65,20 @@
       !! namelist variables
       !!-------------------
-      CHARACTER(len=80)                         ::   filtide             !: Filename root for tidal input files
-      LOGICAL                                   ::   ln_bdytide_2ddta    !: If true, read 2d harmonic data
+      CHARACTER(len=80)                         ::   filtide             ! Filename root for tidal input files
+      LOGICAL                                   ::   ln_bdytide_2ddta    ! If true, read 2d harmonic data
       !!
-      INTEGER                                   ::   ib_bdy, itide, ib   !: dummy loop indices
-      INTEGER                                   ::   ii, ij              !: dummy loop indices
+      INTEGER                                   ::   ib_bdy, itide, ib   ! dummy loop indices
+      INTEGER                                   ::   ii, ij              ! dummy loop indices
       INTEGER                                   ::   inum, igrd
-      INTEGER, DIMENSION(3)                     ::   ilen0       !: length of boundary data (from OBC arrays)
+      INTEGER                                   ::   isz                 ! bdy data size
       INTEGER                                   ::   ios                 ! Local integer output status for namelist read
       INTEGER                                   ::   nbdy_rdstart, nbdy_loc
-      CHARACTER(LEN=50)                         ::   cerrmsg             !: error string
-      CHARACTER(len=80)                         ::   clfile              !: full file name for tidal input file 
-      REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)    ::   dta_read            !: work space to read in tidal harmonics data
-      REAL(wp),ALLOCATABLE, DIMENSION(:,:)      ::   ztr, zti            !:  "     "    "   "   "   "        "      " 
+      CHARACTER(LEN=50)                         ::   cerrmsg             ! error string
+      CHARACTER(len=80)                         ::   clfile              ! full file name for tidal input file 
+      REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)    ::   dta_read            ! work space to read in tidal harmonics data
+      REAL(wp),ALLOCATABLE, DIMENSION(:,:)      ::   ztr, zti            !  "     "    "   "   "   "        "      " 
       !!
-      TYPE(TIDES_DATA),  POINTER                ::   td                  !: local short cut   
+      TYPE(TIDES_DATA), POINTER                 ::   td                  ! local short cut   
+      TYPE(  OBC_DATA), POINTER                 ::   dta                 ! local short cut
       !!
       NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta
@@ -93,6 +94,7 @@
          IF( nn_dyn2d_dta(ib_bdy) >= 2 ) THEN
             !
-            td => tides(ib_bdy)
-
+            td  => tides(ib_bdy)
+            dta => dta_bdy(ib_bdy)
+         
             ! Namelist nambdy_tide : tidal harmonic forcing at open boundaries
             filtide(:) = ''
@@ -130,30 +132,28 @@
             IF(lwp) WRITE(numout,*) ' '
 
-            ! Allocate space for tidal harmonics data - get size from OBC data arrays
+            ! Allocate space for tidal harmonics data - get size from BDY data arrays
+            ! Allocate also slow varying data in the case of time splitting:
+            ! Do it anyway because at this stage knowledge of free surface scheme is unknown
             ! -----------------------------------------------------------------------
-
-            ! JC: If FRS scheme is used, we assume that tidal is needed over the whole
-            ! relaxation area      
-            IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN   ;   ilen0(:) = idx_bdy(ib_bdy)%nblen   (:)
-            ELSE                                   ;   ilen0(:) = idx_bdy(ib_bdy)%nblenrim(:)
-            ENDIF
-
-            ALLOCATE( td%ssh0( ilen0(1), nb_harmo, 2 ) )
-            ALLOCATE( td%ssh ( ilen0(1), nb_harmo, 2 ) )
-
-            ALLOCATE( td%u0( ilen0(2), nb_harmo, 2 ) )
-            ALLOCATE( td%u ( ilen0(2), nb_harmo, 2 ) )
-
-            ALLOCATE( td%v0( ilen0(3), nb_harmo, 2 ) )
-            ALLOCATE( td%v ( ilen0(3), nb_harmo, 2 ) )
-
-            td%ssh0(:,:,:) = 0._wp
-            td%ssh (:,:,:) = 0._wp
-            td%u0  (:,:,:) = 0._wp
-            td%u   (:,:,:) = 0._wp
-            td%v0  (:,:,:) = 0._wp
-            td%v   (:,:,:) = 0._wp
-
+            IF( ASSOCIATED(dta%ssh) ) THEN   ! we use bdy ssh on this mpi subdomain
+               isz = SIZE(dta%ssh)
+               ALLOCATE( td%ssh0( isz, nb_harmo, 2 ), td%ssh( isz, nb_harmo, 2 ), dta_bdy_s(ib_bdy)%ssh( isz ) )
+               dta_bdy_s(ib_bdy)%ssh(:) = 0._wp   ! needed?
+            ENDIF
+            IF( ASSOCIATED(dta%u2d) ) THEN   ! we use bdy u2d on this mpi subdomain
+               isz = SIZE(dta%u2d)
+               ALLOCATE( td%u0  ( isz, nb_harmo, 2 ), td%u  ( isz, nb_harmo, 2 ), dta_bdy_s(ib_bdy)%u2d( isz ) )
+               dta_bdy_s(ib_bdy)%u2d(:) = 0._wp   ! needed?
+            ENDIF
+            IF( ASSOCIATED(dta%v2d) ) THEN   ! we use bdy v2d on this mpi subdomain
+               isz = SIZE(dta%v2d)
+               ALLOCATE( td%v0  ( isz, nb_harmo, 2 ), td%v  ( isz, nb_harmo, 2 ), dta_bdy_s(ib_bdy)%v2d( isz ) )
+               dta_bdy_s(ib_bdy)%v2d(:) = 0._wp   ! needed?
+            ENDIF
+
+            ! fill td%ssh0, td%u0, td%v0
+            ! -----------------------------------------------------------------------
             IF( ln_bdytide_2ddta ) THEN
+               !
                ! It is assumed that each data file contains all complex harmonic amplitudes
                ! given on the global domain (ie global, jpiglo x jpjglo)
@@ -162,53 +162,56 @@
                !
                ! SSH fields
-               clfile = TRIM(filtide)//'_grid_T.nc'
-               CALL iom_open( clfile , inum ) 
-               igrd = 1                       ! Everything is at T-points here
-               DO itide = 1, nb_harmo
-                  CALL iom_get( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_z1', ztr(:,:) )
-                  CALL iom_get( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_z2', zti(:,:) ) 
-                  DO ib = 1, ilen0(igrd)
-                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                     IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove?
-                     td%ssh0(ib,itide,1) = ztr(ii,ij)
-                     td%ssh0(ib,itide,2) = zti(ii,ij)
-                  END DO
-               END DO 
-               CALL iom_close( inum )
+               IF( ASSOCIATED(dta%ssh) ) THEN   ! we use bdy ssh on this mpi subdomain
+                  clfile = TRIM(filtide)//'_grid_T.nc'
+                  CALL iom_open( clfile , inum ) 
+                  igrd = 1                       ! Everything is at T-points here
+                  DO itide = 1, nb_harmo
+                     CALL iom_get( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_z1', ztr(:,:) )
+                     CALL iom_get( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_z2', zti(:,:) ) 
+                     DO ib = 1, SIZE(dta%ssh)
+                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                        td%ssh0(ib,itide,1) = ztr(ii,ij)
+                        td%ssh0(ib,itide,2) = zti(ii,ij)
+                     END DO
+                  END DO
+                  CALL iom_close( inum )
+               ENDIF
                !
                ! U fields
-               clfile = TRIM(filtide)//'_grid_U.nc'
-               CALL iom_open( clfile , inum ) 
-               igrd = 2                       ! Everything is at U-points here
-               DO itide = 1, nb_harmo
-                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_u1', ztr(:,:) )
-                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_u2', zti(:,:) )
-                  DO ib = 1, ilen0(igrd)
-                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                     IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove?
-                     td%u0(ib,itide,1) = ztr(ii,ij)
-                     td%u0(ib,itide,2) = zti(ii,ij)
-                  END DO
-               END DO
-               CALL iom_close( inum )
+               IF( ASSOCIATED(dta%u2d) ) THEN   ! we use bdy u2d on this mpi subdomain
+                  clfile = TRIM(filtide)//'_grid_U.nc'
+                  CALL iom_open( clfile , inum ) 
+                  igrd = 2                       ! Everything is at U-points here
+                  DO itide = 1, nb_harmo
+                     CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_u1', ztr(:,:) )
+                     CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_u2', zti(:,:) )
+                     DO ib = 1, SIZE(dta%u2d)
+                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                        td%u0(ib,itide,1) = ztr(ii,ij)
+                        td%u0(ib,itide,2) = zti(ii,ij)
+                     END DO
+                  END DO
+                  CALL iom_close( inum )
+               ENDIF
                !
                ! V fields
-               clfile = TRIM(filtide)//'_grid_V.nc'
-               CALL iom_open( clfile , inum ) 
-               igrd = 3                       ! Everything is at V-points here
-               DO itide = 1, nb_harmo
-                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_v1', ztr(:,:) )
-                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_v2', zti(:,:) )
-                  DO ib = 1, ilen0(igrd)
-                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                     IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove?
-                     td%v0(ib,itide,1) = ztr(ii,ij)
-                     td%v0(ib,itide,2) = zti(ii,ij)
-                  END DO
-               END DO  
-               CALL iom_close( inum )
+               IF( ASSOCIATED(dta%v2d) ) THEN   ! we use bdy v2d on this mpi subdomain
+                  clfile = TRIM(filtide)//'_grid_V.nc'
+                  CALL iom_open( clfile , inum ) 
+                  igrd = 3                       ! Everything is at V-points here
+                  DO itide = 1, nb_harmo
+                     CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_v1', ztr(:,:) )
+                     CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_v2', zti(:,:) )
+                     DO ib = 1, SIZE(dta%v2d)
+                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                        td%v0(ib,itide,1) = ztr(ii,ij)
+                        td%v0(ib,itide,2) = zti(ii,ij)
+                     END DO
+                  END DO
+                  CALL iom_close( inum )
+               ENDIF
                !
                DEALLOCATE( ztr, zti ) 
@@ -218,5 +221,5 @@
                ! Read tidal data only on bdy segments
                ! 
-               ALLOCATE( dta_read( MAXVAL(ilen0(1:3)), 1, 1 ) )
+               ALLOCATE( dta_read( MAXVAL( idx_bdy(ib_bdy)%nblen(:) ), 1, 1 ) )
                !
                ! Open files and read in tidal forcing data
@@ -225,27 +228,36 @@
                DO itide = 1, nb_harmo
                   !                                                              ! SSH fields
-                  clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_T.nc'
-                  CALL iom_open( clfile, inum )
-                  CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
-                  td%ssh0(:,itide,1) = dta_read(1:ilen0(1),1,1)
-                  CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
-                  td%ssh0(:,itide,2) = dta_read(1:ilen0(1),1,1)
-                  CALL iom_close( inum )
+                  IF( ASSOCIATED(dta%ssh) ) THEN   ! we use bdy ssh on this mpi subdomain
+                     isz = SIZE(dta%ssh)
+                     clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_T.nc'
+                     CALL iom_open( clfile, inum )
+                     CALL fld_map( inum, 'z1', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
+                     td%ssh0(:,itide,1) = dta_read(1:isz,1,1)
+                     CALL fld_map( inum, 'z2', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
+                     td%ssh0(:,itide,2) = dta_read(1:isz,1,1)
+                     CALL iom_close( inum )
+                  ENDIF
                   !                                                              ! U fields
-                  clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_U.nc'
-                  CALL iom_open( clfile, inum )
-                  CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
-                  td%u0(:,itide,1) = dta_read(1:ilen0(2),1,1)
-                  CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
-                  td%u0(:,itide,2) = dta_read(1:ilen0(2),1,1)
-                  CALL iom_close( inum )
+                  IF( ASSOCIATED(dta%u2d) ) THEN   ! we use bdy u2d on this mpi subdomain
+                     isz = SIZE(dta%u2d)
+                     clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_U.nc'
+                     CALL iom_open( clfile, inum )
+                     CALL fld_map( inum, 'u1', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
+                     td%u0(:,itide,1) = dta_read(1:isz,1,1)
+                     CALL fld_map( inum, 'u2', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
+                     td%u0(:,itide,2) = dta_read(1:isz,1,1)
+                     CALL iom_close( inum )
+                  ENDIF
                   !                                                              ! V fields
-                  clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_V.nc'
-                  CALL iom_open( clfile, inum )
-                  CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
-                  td%v0(:,itide,1) = dta_read(1:ilen0(3),1,1)
-                  CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
-                  td%v0(:,itide,2) = dta_read(1:ilen0(3),1,1)
-                  CALL iom_close( inum )
+                  IF( ASSOCIATED(dta%v2d) ) THEN   ! we use bdy v2d on this mpi subdomain
+                     isz = SIZE(dta%v2d)
+                     clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_V.nc'
+                     CALL iom_open( clfile, inum )
+                     CALL fld_map( inum, 'v1', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
+                     td%v0(:,itide,1) = dta_read(1:isz,1,1)
+                     CALL fld_map( inum, 'v2', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
+                     td%v0(:,itide,2) = dta_read(1:isz,1,1)
+                     CALL iom_close( inum )
+                  ENDIF
                   !
                END DO ! end loop on tidal components
@@ -254,13 +266,4 @@
                !
             ENDIF ! ln_bdytide_2ddta=.true.
-            !
-            ! Allocate slow varying data in the case of time splitting:
-            ! Do it anyway because at this stage knowledge of free surface scheme is unknown
-            ALLOCATE( dta_bdy_s(ib_bdy)%ssh ( ilen0(1) ) )
-            ALLOCATE( dta_bdy_s(ib_bdy)%u2d ( ilen0(2) ) )
-            ALLOCATE( dta_bdy_s(ib_bdy)%v2d ( ilen0(3) ) )
-            dta_bdy_s(ib_bdy)%ssh(:) = 0._wp
-            dta_bdy_s(ib_bdy)%u2d(:) = 0._wp
-            dta_bdy_s(ib_bdy)%v2d(:) = 0._wp
             !
          ENDIF ! nn_dyn2d_dta(ib_bdy) >= 2
@@ -283,7 +286,5 @@
       !
       LOGICAL  ::   lk_first_btstp            ! =.TRUE. if time splitting and first barotropic step
-      INTEGER  ::   itide, ib_bdy, ib, igrd   ! loop indices
-      INTEGER, DIMENSION(jpbgrd)   ::   ilen0 
-      INTEGER, DIMENSION(1:jpbgrd) ::   nblen, nblenrim  ! short cuts
+      INTEGER  ::   itide, ib_bdy, ib         ! loop indices
       REAL(wp) ::   z_arg, z_sarg, zramp, zoff, z_cost, z_sist, zt_offset   
       !!----------------------------------------------------------------------
@@ -310,11 +311,4 @@
          IF( nn_dyn2d_dta(ib_bdy) >= 2 ) THEN
             !
-            nblen(1:jpbgrd) = idx_bdy(ib_bdy)%nblen(1:jpbgrd)
-            nblenrim(1:jpbgrd) = idx_bdy(ib_bdy)%nblenrim(1:jpbgrd)
-            !
-            IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN   ;   ilen0(:) = nblen   (:)
-            ELSE                                   ;   ilen0(:) = nblenrim(:)
-            ENDIF     
-            !
             ! We refresh nodal factors every day below
             ! This should be done somewhere else
@@ -337,7 +331,7 @@
             ! If time splitting, initialize arrays from slow varying open boundary data:
             IF ( PRESENT(kit) ) THEN           
-               IF ( dta_bdy(ib_bdy)%lneed_ssh   ) dta_bdy(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy_s(ib_bdy)%ssh(1:ilen0(1))
-               IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) dta_bdy(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy_s(ib_bdy)%u2d(1:ilen0(2))
-               IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) dta_bdy(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy_s(ib_bdy)%v2d(1:ilen0(3))
+               IF ( ASSOCIATED(dta_bdy(ib_bdy)%ssh) ) dta_bdy(ib_bdy)%ssh(:) = dta_bdy_s(ib_bdy)%ssh(:)
+               IF ( ASSOCIATED(dta_bdy(ib_bdy)%u2d) ) dta_bdy(ib_bdy)%u2d(:) = dta_bdy_s(ib_bdy)%u2d(:)
+               IF ( ASSOCIATED(dta_bdy(ib_bdy)%v2d) ) dta_bdy(ib_bdy)%v2d(:) = dta_bdy_s(ib_bdy)%v2d(:)
             ENDIF
             !
@@ -349,7 +343,6 @@
                z_sist = zramp * SIN( z_sarg )
                !
-               IF ( dta_bdy(ib_bdy)%lneed_ssh ) THEN
-                  igrd=1                              ! SSH on tracer grid
-                  DO ib = 1, ilen0(igrd)
+               IF ( ASSOCIATED(dta_bdy(ib_bdy)%ssh) ) THEN   ! SSH on tracer grid
+                  DO ib = 1, SIZE(dta_bdy(ib_bdy)%ssh)
                      dta_bdy(ib_bdy)%ssh(ib) = dta_bdy(ib_bdy)%ssh(ib) + &
                         &                      ( tides(ib_bdy)%ssh(ib,itide,1)*z_cost + &
@@ -358,13 +351,14 @@
                ENDIF
                !
-               IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) THEN
-                  igrd=2                              ! U grid
-                  DO ib = 1, ilen0(igrd)
+               IF ( ASSOCIATED(dta_bdy(ib_bdy)%u2d) ) THEN  ! U grid
+                  DO ib = 1, SIZE(dta_bdy(ib_bdy)%u2d)
                      dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) + &
                         &                      ( tides(ib_bdy)%u(ib,itide,1)*z_cost + &
                         &                        tides(ib_bdy)%u(ib,itide,2)*z_sist )
                   END DO
-                  igrd=3                              ! V grid
-                  DO ib = 1, ilen0(igrd) 
+               ENDIF
+               !
+               IF ( ASSOCIATED(dta_bdy(ib_bdy)%v2d) ) THEN   ! V grid
+                  DO ib = 1, SIZE(dta_bdy(ib_bdy)%v2d)
                      dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) + &
                         &                      ( tides(ib_bdy)%v(ib,itide,1)*z_cost + &
@@ -372,6 +366,7 @@
                   END DO
                ENDIF
+               !
             END DO             
-         END IF
+         ENDIF
       END DO
       !
@@ -386,31 +381,31 @@
       TYPE(TIDES_DATA), INTENT(inout) ::   td    ! tidal harmonics data
       !
-      INTEGER ::   itide, igrd, ib       ! dummy loop indices
-      INTEGER, DIMENSION(1) ::   ilen0   ! length of boundary data (from OBC arrays)
+      INTEGER ::   itide, isz, ib       ! dummy loop indices
       REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide
       !!----------------------------------------------------------------------
       !
-      igrd=1   
-                              ! SSH on tracer grid.
-      ilen0(1) =  SIZE(td%ssh0(:,1,1))
-      !
-      ALLOCATE( mod_tide(ilen0(igrd)), phi_tide(ilen0(igrd)) )
-      !
-      DO itide = 1, nb_harmo
-         DO ib = 1, ilen0(igrd)
-            mod_tide(ib)=SQRT(td%ssh0(ib,itide,1)**2.+td%ssh0(ib,itide,2)**2.)
-            phi_tide(ib)=ATAN2(-td%ssh0(ib,itide,2),td%ssh0(ib,itide,1))
+      IF( ASSOCIATED(td%ssh0) ) THEN   ! SSH on tracer grid.
+         !
+         isz = SIZE( td%ssh0, dim = 1 )
+         ALLOCATE( mod_tide(isz), phi_tide(isz) )
+         !
+         DO itide = 1, nb_harmo
+            DO ib = 1, isz
+               mod_tide(ib)=SQRT( td%ssh0(ib,itide,1)*td%ssh0(ib,itide,1) + td%ssh0(ib,itide,2)*td%ssh0(ib,itide,2) )
+               phi_tide(ib)=ATAN2(-td%ssh0(ib,itide,2),td%ssh0(ib,itide,1))
+            END DO
+            DO ib = 1, isz
+               mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f
+               phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0+tide_harmonics(itide)%u
+            END DO
+            DO ib = 1, isz
+               td%ssh(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
+               td%ssh(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
+            END DO
          END DO
-         DO ib = 1 , ilen0(igrd)
-            mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f
-            phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0+tide_harmonics(itide)%u
-         ENDDO
-         DO ib = 1 , ilen0(igrd)
-            td%ssh(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
-            td%ssh(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
-         ENDDO
-      END DO
-      !
-      DEALLOCATE( mod_tide, phi_tide )
+         !
+         DEALLOCATE( mod_tide, phi_tide )
+         !
+      ENDIF
       !
    END SUBROUTINE tide_init_elevation
@@ -424,55 +419,57 @@
       TYPE(TIDES_DATA), INTENT(inout) ::   td    ! tidal harmonics data
       !
-      INTEGER ::   itide, igrd, ib       ! dummy loop indices
-      INTEGER, DIMENSION(3) ::   ilen0   ! length of boundary data (from OBC arrays)
+      INTEGER ::   itide, isz, ib        ! dummy loop indices
       REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide
       !!----------------------------------------------------------------------
       !
-      ilen0(2) =  SIZE(td%u0(:,1,1))
-      ilen0(3) =  SIZE(td%v0(:,1,1))
-      !
-      igrd=2                                 ! U grid.
-      !
-      ALLOCATE( mod_tide(ilen0(igrd)) , phi_tide(ilen0(igrd)) )
-      !
-      DO itide = 1, nb_harmo
-         DO ib = 1, ilen0(igrd)
-            mod_tide(ib)=SQRT(td%u0(ib,itide,1)**2.+td%u0(ib,itide,2)**2.)
-            phi_tide(ib)=ATAN2(-td%u0(ib,itide,2),td%u0(ib,itide,1))
+      IF( ASSOCIATED(td%u0) ) THEN   ! U grid. we use bdy u2d on this mpi subdomain
+         !
+         isz = SIZE( td%u0, dim = 1 )
+         ALLOCATE( mod_tide(isz), phi_tide(isz) )
+         !
+         DO itide = 1, nb_harmo
+            DO ib = 1, isz
+               mod_tide(ib)=SQRT( td%u0(ib,itide,1)*td%u0(ib,itide,1) + td%u0(ib,itide,2)*td%u0(ib,itide,2) )
+               phi_tide(ib)=ATAN2(-td%u0(ib,itide,2),td%u0(ib,itide,1))
+            END DO
+            DO ib = 1, isz
+               mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f
+               phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0 + tide_harmonics(itide)%u
+            END DO
+            DO ib = 1, isz
+               td%u(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
+               td%u(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
+            END DO
          END DO
-         DO ib = 1, ilen0(igrd)
-            mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f
-            phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0 + tide_harmonics(itide)%u
-         ENDDO
-         DO ib = 1, ilen0(igrd)
-            td%u(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
-            td%u(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
-         ENDDO
-      END DO
-      !
-      DEALLOCATE( mod_tide , phi_tide )
-      !
-      igrd=3                                 ! V grid.
-      !
-      ALLOCATE( mod_tide(ilen0(igrd)) , phi_tide(ilen0(igrd)) )
-
-      DO itide = 1, nb_harmo
-         DO ib = 1, ilen0(igrd)
-            mod_tide(ib)=SQRT(td%v0(ib,itide,1)**2.+td%v0(ib,itide,2)**2.)
-            phi_tide(ib)=ATAN2(-td%v0(ib,itide,2),td%v0(ib,itide,1))
+         !
+         DEALLOCATE( mod_tide, phi_tide )
+         !
+      ENDIF
+      !
+      IF( ASSOCIATED(td%v0) ) THEN   ! V grid. we use bdy u2d on this mpi subdomain
+         !
+         isz = SIZE( td%v0, dim = 1 )
+         ALLOCATE( mod_tide(isz), phi_tide(isz) )
+         !
+         DO itide = 1, nb_harmo
+            DO ib = 1, isz
+               mod_tide(ib)=SQRT( td%v0(ib,itide,1)*td%v0(ib,itide,1) + td%v0(ib,itide,2)*td%v0(ib,itide,2) )
+               phi_tide(ib)=ATAN2(-td%v0(ib,itide,2),td%v0(ib,itide,1))
+            END DO
+            DO ib = 1, isz
+               mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f
+               phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0 + tide_harmonics(itide)%u
+            END DO
+            DO ib = 1, isz
+               td%v(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
+               td%v(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
+            END DO
          END DO
-         DO ib = 1, ilen0(igrd)
-            mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f
-            phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0 + tide_harmonics(itide)%u
-         ENDDO
-         DO ib = 1, ilen0(igrd)
-            td%v(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
-            td%v(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
-         ENDDO
-      END DO
-      !
-      DEALLOCATE( mod_tide, phi_tide )
-      !
-  END SUBROUTINE tide_init_velocities
+         !
+         DEALLOCATE( mod_tide, phi_tide )
+         !
+      ENDIF
+      !
+   END SUBROUTINE tide_init_velocities
 
    !!======================================================================
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/C1D/step_c1d.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/C1D/step_c1d.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/C1D/step_c1d.F90	(revision 13193)
@@ -27,5 +27,5 @@
    PRIVATE
 
-   PUBLIC stp_c1d      ! called by opa.F90
+   PUBLIC stp_c1d      ! called by nemogcm.F90
 
    !!----------------------------------------------------------------------
@@ -56,8 +56,5 @@
       !
       INTEGER ::   jk       ! dummy loop indice
-      INTEGER ::   indic    ! error indicator if < 0
       !! ---------------------------------------------------------------------
-
-                             indic = 0                ! reset to no error condition
       IF( kstp == nit000 )   CALL iom_init( "nemo")   ! iom_put initialization (must be done after nemo_init for AGRIF+XIOS+OASIS)
       IF( kstp /= nit000 )   CALL day( kstp )         ! Calendar (day was already called at nit000 in day_init)
@@ -88,5 +85,5 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
                          CALL dia_wri( kstp, Nnn )  ! ocean model: outputs
-      IF( lk_diahth  )   CALL dia_hth( kstp, Nnn )  ! Thermocline depth (20°C)
+                         CALL dia_hth( kstp, Nnn )  ! Thermocline depth (20°C)
 
 
@@ -111,7 +108,5 @@
                         CALL eos( ts(:,:,:,:,Nnn), rhd, rhop, gdept_0(:,:,:) )  ! now potential density for zdfmxl
       IF( ln_zdfnpc )   CALL tra_npc( kstp,      Nnn, Nrhs, ts, Naa   )         ! applied non penetrative convective adjustment on (t,s)
-                        CALL tra_atf( kstp, Nbb, Nnn, Nrhs,     Naa, ts   )     ! time filtering of "now" tracer fields
-
-
+                        CALL tra_atf( kstp, Nbb, Nnn, Naa, ts )                 ! time filtering of "now" tracer arrays
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -139,10 +134,10 @@
       ! Control and restarts
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                             CALL stp_ctl( kstp, Nnn, indic )
+                             CALL stp_ctl( kstp, Nnn )
       IF( kstp == nit000 )   CALL iom_close( numror )          ! close input  ocean restart file
       IF( lrst_oce       )   CALL rst_write( kstp, Nbb, Nnn )  ! write output ocean restart file
       !
 #if defined key_iomput
-      IF( kstp == nitend .OR. indic < 0 )   CALL xios_context_finalize()   ! needed for XIOS
+      IF( kstp == nitend .OR. nstop > 0 )   CALL xios_context_finalize()   ! needed for XIOS
       !
 #endif
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DIA/diaar5.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DIA/diaar5.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DIA/diaar5.F90	(revision 13193)
@@ -76,7 +76,7 @@
       REAL(wp) ::   zaw, zbw, zrw
       !
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     :: zarea_ssh , zbotpres       ! 2D workspace 
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     :: zarea_ssh, zbotpres        ! 2D workspace 
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)     :: z2d, zpe                   ! 2D workspace 
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   :: z3d, zrhd , zrhop, ztpot, zgdept   ! 3D workspace (zgdept: needed to use the substitute)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   :: z3d, zrhd, ztpot, zgdept   ! 3D workspace (zgdept: needed to use the substitute)
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ztsn                       ! 4D workspace
 
@@ -88,5 +88,5 @@
       IF( l_ar5 ) THEN 
          ALLOCATE( zarea_ssh(jpi,jpj), zbotpres(jpi,jpj), z2d(jpi,jpj) )
-         ALLOCATE( zrhd(jpi,jpj,jpk) , zrhop(jpi,jpj,jpk) )
+         ALLOCATE( zrhd(jpi,jpj,jpk) )
          ALLOCATE( ztsn(jpi,jpj,jpk,jpts) )
          zarea_ssh(:,:) = e1e2t(:,:) * ssh(:,:,Kmm)
@@ -163,11 +163,7 @@
       
          !                                         ! steric sea surface height
-         CALL eos( ts(:,:,:,:,Kmm), zrhd, zrhop, zgdept )                 ! now in situ and potential density
-         zrhop(:,:,jpk) = 0._wp
-         CALL iom_put( 'rhop', zrhop )
-         !
          zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice
          DO jk = 1, jpkm1
-            zbotpres(:,:) = zbotpres(:,:) + e3t(:,:,jk,Kmm) * zrhd(:,:,jk)
+            zbotpres(:,:) = zbotpres(:,:) + e3t(:,:,jk,Kmm) * rhd(:,:,jk)
          END DO
          IF( ln_linssh ) THEN
@@ -176,9 +172,9 @@
                   DO jj = 1,jpj
                      iks = mikt(ji,jj)
-                     zbotpres(ji,jj) = zbotpres(ji,jj) + ssh(ji,jj,Kmm) * zrhd(ji,jj,iks) + riceload(ji,jj)
+                     zbotpres(ji,jj) = zbotpres(ji,jj) + ssh(ji,jj,Kmm) * rhd(ji,jj,iks) + riceload(ji,jj)
                   END DO
                END DO
             ELSE
-               zbotpres(:,:) = zbotpres(:,:) + ssh(:,:,Kmm) * zrhd(:,:,1)
+               zbotpres(:,:) = zbotpres(:,:) + ssh(:,:,Kmm) * rhd(:,:,1)
             END IF
          END IF
@@ -303,5 +299,4 @@
       IF( l_ar5 ) THEN
         DEALLOCATE( zarea_ssh , zbotpres, z2d )
-        DEALLOCATE( zrhd      , zrhop    )
         DEALLOCATE( ztsn                 )
       ENDIF
@@ -377,5 +372,6 @@
       IF(   iom_use( 'voltot'  ) .OR. iom_use( 'sshtot'    )  .OR. iom_use( 'sshdyn' )  .OR.  & 
          &  iom_use( 'masstot' ) .OR. iom_use( 'temptot'   )  .OR. iom_use( 'saltot' ) .OR.  &    
-         &  iom_use( 'botpres' ) .OR. iom_use( 'sshthster' )  .OR. iom_use( 'sshsteric' )  ) L_ar5 = .TRUE.
+         &  iom_use( 'botpres' ) .OR. iom_use( 'sshthster' )  .OR. iom_use( 'sshsteric' ) .OR. &
+         &  iom_use( 'rhop' )  ) L_ar5 = .TRUE.
   
       IF( l_ar5 ) THEN
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DIA/diawri.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DIA/diawri.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DIA/diawri.F90	(revision 13193)
@@ -191,4 +191,6 @@
          CALL iom_put( "sbs", z2d )                ! bottom salinity
       ENDIF
+
+      CALL iom_put( "rhop", rhop(:,:,:) )          ! 3D potential density (sigma0)
 
       IF ( iom_use("taubot") ) THEN                ! bottom stress
@@ -1033,6 +1035,6 @@
          CALL iom_close( inum )
       ENDIF
+      !
 #endif
-
    END SUBROUTINE dia_wri_state
 
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/dom_oce.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/dom_oce.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/dom_oce.F90	(revision 13193)
@@ -17,4 +17,5 @@
    !!----------------------------------------------------------------------
    !!   Agrif_Root    : dummy function used when lk_agrif=F
+   !!   Agrif_Fixed   : dummy function used when lk_agrif=F
    !!   Agrif_CFixed  : dummy function used when lk_agrif=F
    !!   dom_oce_alloc : dynamical allocation of dom_oce arrays
@@ -243,4 +244,8 @@
    END FUNCTION Agrif_Root
 
+   INTEGER FUNCTION Agrif_Fixed()
+      Agrif_Fixed = 0
+   END FUNCTION Agrif_Fixed
+
    CHARACTER(len=3) FUNCTION Agrif_CFixed()
       Agrif_CFixed = '0'
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/dommsk.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/dommsk.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/dommsk.F90	(revision 13193)
@@ -261,12 +261,4 @@
                ENDIF
             END DO
-#if defined key_agrif
-            IF( .NOT. AGRIF_Root() ) THEN
-               IF ((nbondi ==  1).OR.(nbondi == 2)) fmask(nlci-1 , :     ,jk) = 0.e0      ! east
-               IF ((nbondi == -1).OR.(nbondi == 2)) fmask(1      , :     ,jk) = 0.e0      ! west
-               IF ((nbondj ==  1).OR.(nbondj == 2)) fmask(:      ,nlcj-1 ,jk) = 0.e0      ! north
-               IF ((nbondj == -1).OR.(nbondj == 2)) fmask(:      ,1      ,jk) = 0.e0      ! south
-            ENDIF
-#endif
          END DO
          !
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domvvl.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domvvl.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domvvl.F90	(revision 13193)
@@ -915,11 +915,9 @@
                e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
 
-               DO ji = 1, jpi
-                  DO jj = 1, jpj
-                     IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN
-                       CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' )
-                     ENDIF
-                  END DO
-               END DO
+               DO_2D_11_11
+                  IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN
+                     CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' )
+                  ENDIF
+               END_2D
                !
             ELSE
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/istate.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/istate.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/istate.F90	(revision 13193)
@@ -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)
@@ -126,7 +125,4 @@
          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 :
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/divhor.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/divhor.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/divhor.F90	(revision 13193)
@@ -85,13 +85,4 @@
       END_3D
       !
-#if defined key_agrif
-      IF( .NOT. Agrif_Root() ) THEN
-         IF( nbondi == -1 .OR. nbondi == 2 )   hdiv(   2   ,  :   ,:) = 0._wp      ! west
-         IF( nbondi ==  1 .OR. nbondi == 2 )   hdiv( nlci-1,  :   ,:) = 0._wp      ! east
-         IF( nbondj == -1 .OR. nbondj == 2 )   hdiv(   :   ,  2   ,:) = 0._wp      ! south
-         IF( nbondj ==  1 .OR. nbondj == 2 )   hdiv(   :   ,nlcj-1,:) = 0._wp      ! north
-      ENDIF
-#endif
-      !
       IF( ln_rnf )   CALL sbc_rnf_div( hdiv, Kmm )                     !==  runoffs    ==!   (update hdiv field)
       !
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynldf_lap_blp.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynldf_lap_blp.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynldf_lap_blp.F90	(revision 13193)
@@ -75,30 +75,21 @@
          DO_2D_01_01
             !                                      ! ahm * e3 * curl  (computed from 1 to jpim1/jpjm1)
-!!gm open question here : e3f  at before or now ?    probably now...
-!!gm note that ahmf has already been multiplied by fmask
-            zcur(ji-1,jj-1) =  &
-               &      ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)      &
-               &  * (  e2v(ji  ,jj-1) * pv(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk)  &
-               &     - e1u(ji-1,jj  ) * pu(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk)  )
+            zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)       &   ! ahmf already * by fmask
+               &     * (  e2v(ji  ,jj-1) * pv(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk)  &
+               &        - e1u(ji-1,jj  ) * pu(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk)  )
             !                                      ! ahm * div        (computed from 2 to jpi/jpj)
-!!gm note that ahmt has already been multiplied by tmask
-            zdiv(ji,jj)     =   &
-               &   ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb)      &
-               &     * (  e2u(ji,jj)*e3u(ji,jj,jk,Kbb) * pu(ji,jj,jk)        &
-               &        - e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb) * pu(ji-1,jj,jk)  &
-               &        + e1v(ji,jj)*e3v(ji,jj,jk,Kbb) * pv(ji,jj,jk)        &
-               &        - e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kbb) * pv(ji,jj-1,jk)  )
+            zdiv(ji,jj)     = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb)               &   ! ahmt already * by tmask
+               &     * (  e2u(ji,jj)*e3u(ji,jj,jk,Kbb) * pu(ji,jj,jk) - e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb) * pu(ji-1,jj,jk)  &
+               &        + e1v(ji,jj)*e3v(ji,jj,jk,Kbb) * pv(ji,jj,jk) - e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kbb) * pv(ji,jj-1,jk)  )
          END_2D
          !
          DO_2D_00_00
-            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * (                             &
-               &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj)   &
-               &              / e3u(ji,jj,jk,Kmm)   &
-               &              + ( zdiv(ji+1,jj) - zdiv(ji,jj  ) ) * r1_e1u(ji,jj)        )
+            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * umask(ji,jj,jk) * (    &    ! * by umask is mandatory for dyn_ldf_blp use
+               &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u(ji,jj,jk,Kmm)   &
+               &              + ( zdiv(ji+1,jj) - zdiv(ji,jj  ) ) * r1_e1u(ji,jj)                      )
                !
-            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * (                              &
-               &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj)   &
-               &              / e3v(ji,jj,jk,Kmm)   &
-               &              + ( zdiv(ji,jj+1) - zdiv(ji  ,jj) ) * r1_e2v(ji,jj)       )
+            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * vmask(ji,jj,jk) * (    &    ! * by vmask is mandatory for dyn_ldf_blp use
+               &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v(ji,jj,jk,Kmm)   &
+               &              + ( zdiv(ji,jj+1) - zdiv(ji  ,jj) ) * r1_e2v(ji,jj)                      )
          END_2D
          !                                             ! ===============
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynvor.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynvor.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/dynvor.F90	(revision 13193)
@@ -820,5 +820,5 @@
          DO_3D_10_10( 1, jpk )
             IF(    tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk)              &
-               & + tmask(ji,jj  ,jk) + tmask(ji+1,jj+1,jk) == 3._wp )   fmask(ji,jj,jk) = 1._wp
+               & + tmask(ji,jj  ,jk) + tmask(ji+1,jj  ,jk) == 3._wp )   fmask(ji,jj,jk) = 1._wp
          END_3D
          !
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/sshwzv.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/sshwzv.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/sshwzv.F90	(revision 13193)
@@ -215,12 +215,37 @@
       ENDIF
       !
-#if defined key_agrif
-      IF( .NOT. AGRIF_Root() ) THEN
-         IF ((nbondi ==  1).OR.(nbondi == 2)) pww(nlci-1 , :     ,:) = 0.e0      ! east
-         IF ((nbondi == -1).OR.(nbondi == 2)) pww(2      , :     ,:) = 0.e0      ! west
-         IF ((nbondj ==  1).OR.(nbondj == 2)) pww(:      ,nlcj-1 ,:) = 0.e0      ! north
-         IF ((nbondj == -1).OR.(nbondj == 2)) pww(:      ,2      ,:) = 0.e0      ! south
-      ENDIF
-#endif
+#if defined key_agrif 
+      IF( .NOT. AGRIF_Root() ) THEN 
+         ! Mask vertical velocity at first/last columns/row 
+         ! inside computational domain (cosmetic) 
+         ! --- West --- !
+         DO ji = mi0(2), mi1(2)
+            DO jj = 1, jpj
+               pww(ji,jj,:) = 0._wp 
+            ENDDO
+         ENDDO
+         !
+         ! --- East --- !
+         DO ji = mi0(jpiglo-1), mi1(jpiglo-1)
+            DO jj = 1, jpj
+               pww(ji,jj,:) = 0._wp
+            ENDDO
+         ENDDO
+         !
+         ! --- South --- !
+         DO jj = mj0(2), mj1(2)
+            DO ji = 1, jpi
+               pww(ji,jj,:) = 0._wp
+            ENDDO
+         ENDDO
+         !
+         ! --- North --- !
+         DO jj = mj0(jpjglo-1), mj1(jpjglo-1)
+            DO ji = 1, jpi
+               pww(ji,jj,:) = 0._wp
+            ENDDO
+         ENDDO
+      ENDIF 
+#endif 
       !
       IF( ln_timing )   CALL timing_stop('wzv')
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ICB/icbrst.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ICB/icbrst.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ICB/icbrst.F90	(revision 13193)
@@ -188,8 +188,10 @@
       !
       INTEGER ::   jn   ! dummy loop index
+      INTEGER ::   idg  ! number of digits
       INTEGER ::   ix_dim, iy_dim, ik_dim, in_dim
       CHARACTER(len=256)     :: cl_path
       CHARACTER(len=256)     :: cl_filename
-      CHARACTER(len=256)     :: cl_kt
+      CHARACTER(len=8  )     :: cl_kt
+      CHARACTER(LEN=12 )     :: clfmt            ! writing format
       TYPE(iceberg), POINTER :: this
       TYPE(point)  , POINTER :: pt
@@ -211,9 +213,11 @@
          ! file name
          WRITE(cl_kt, '(i8.8)') kt
-         cl_filename = TRIM(cexper)//"_"//TRIM(ADJUSTL(cl_kt))//"_"//TRIM(cn_icbrst_out)
+         cl_filename = TRIM(cexper)//"_"//cl_kt//"_"//TRIM(cn_icbrst_out)
          IF( lk_mpp ) THEN
-            WRITE(cl_filename,'(A,"_",I4.4,".nc")') TRIM(cl_filename), narea-1
+            idg = MAX( INT(LOG10(REAL(MAX(1,jpnij-1),wp))) + 1, 4 )          ! how many digits to we need to write? min=4, max=9
+            WRITE(clfmt, "('(a,a,i', i1, '.', i1, ',a)')") idg, idg          ! '(a,a,ix.x,a)'
+            WRITE(cl_filename,  clfmt) TRIM(cl_filename), '_', narea-1, '.nc'
          ELSE
-            WRITE(cl_filename,'(A,".nc")') TRIM(cl_filename)
+            WRITE(cl_filename,'(a,a)') TRIM(cl_filename),               '.nc'
          ENDIF
 
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ICB/icbtrj.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ICB/icbtrj.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ICB/icbtrj.F90	(revision 13193)
@@ -62,7 +62,9 @@
       !
       INTEGER                ::   iret, iyear, imonth, iday
+      INTEGER                ::   idg  ! number of digits
       REAL(wp)               ::   zfjulday, zsec
       CHARACTER(len=80)      ::   cl_filename
-      CHARACTER(LEN=20)      ::   cldate_ini, cldate_end
+      CHARACTER(LEN=12)      ::   clfmt            ! writing format
+      CHARACTER(LEN=8 )      ::   cldate_ini, cldate_end
       TYPE(iceberg), POINTER ::   this
       TYPE(point)  , POINTER ::   pt
@@ -80,8 +82,11 @@
 
       ! define trajectory output name
-      IF ( lk_mpp ) THEN   ;   WRITE(cl_filename,'("trajectory_icebergs_",A,"-",A,"_",I4.4,".nc")')   &
-         &                        TRIM(ADJUSTL(cldate_ini)), TRIM(ADJUSTL(cldate_end)), narea-1
-      ELSE                 ;   WRITE(cl_filename,'("trajectory_icebergs_",A,"-",A         ,".nc")')   &
-         &                        TRIM(ADJUSTL(cldate_ini)), TRIM(ADJUSTL(cldate_end))
+      cl_filename = 'trajectory_icebergs_'//cldate_ini//'-'//cldate_end
+      IF ( lk_mpp ) THEN
+         idg = MAX( INT(LOG10(REAL(MAX(1,jpnij-1),wp))) + 1, 4 )          ! how many digits to we need to write? min=4, max=9
+         WRITE(clfmt, "('(a,a,i', i1, '.', i1, ',a)')") idg, idg          ! '(a,a,ix.x,a)'
+         WRITE(cl_filename,  clfmt) TRIM(cl_filename), '_', narea-1, '.nc'
+      ELSE
+         WRITE(cl_filename,'(a,a)') TRIM(cl_filename),               '.nc'
       ENDIF
       IF( lwp .AND. nn_verbose_level >= 0 )   WRITE(numout,'(2a)') 'icebergs, icb_trj_init: creating ',TRIM(cl_filename)
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/IOM/in_out_manager.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/IOM/in_out_manager.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/IOM/in_out_manager.F90	(revision 13193)
@@ -100,14 +100,4 @@
    !!----------------------------------------------------------------------
    TYPE :: sn_ctl                !: structure for control over output selection
-      LOGICAL :: l_glochk  = .FALSE.  !: range sanity checks are local (F) or global (T)
-                                      !  Use global setting for debugging only;
-                                      !  local breaches will still be reported
-                                      !  and stop the code in most cases.
-      LOGICAL :: l_allon   = .FALSE.  !: overall control; activate all following output options
-      LOGICAL :: l_config  = .FALSE.  !: activate/deactivate finer control
-                                      !  Note if l_config is True then sn_cfctl%l_allon is ignored.
-                                      !  Otherwise setting sn_cfctl%l_allon T/F is equivalent to 
-                                      !  setting all the following logicals in this structure T/F
-                                      !  and disabling subsetting of processors
       LOGICAL :: l_runstat = .FALSE.  !: Produce/do not produce run.stat file (T/F)
       LOGICAL :: l_trcstat = .FALSE.  !: Produce/do not produce tracer.stat file (T/F)
@@ -169,4 +159,7 @@
    INTEGER       ::   no_print = 0          !: optional argument of fld_fill (if present, suppress some control print)
    INTEGER       ::   nstop = 0             !: error flag (=number of reason for a premature stop run)
+!$AGRIF_DO_NOT_TREAT
+   INTEGER       ::   ngrdstop = -1         !: grid number having nstop > 1
+!$AGRIF_END_DO_NOT_TREAT
    INTEGER       ::   nwarn = 0             !: warning flag (=number of warning found during the run)
    CHARACTER(lc) ::   ctmp1, ctmp2, ctmp3   !: temporary characters 1 to 3
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/IOM/iom_def.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/IOM/iom_def.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/IOM/iom_def.F90	(revision 13193)
@@ -33,5 +33,5 @@
    INTEGER, PARAMETER, PUBLIC ::   jpmax_vars   = 1200 !: maximum number of variables in one file
    INTEGER, PARAMETER, PUBLIC ::   jpmax_dims   =  4   !: maximum number of dimensions for one variable
-   INTEGER, PARAMETER, PUBLIC ::   jpmax_digits =  5   !: maximum number of digits for the cpu number in the file name
+   INTEGER, PARAMETER, PUBLIC ::   jpmax_digits =  9   !: maximum number of digits for the cpu number in the file name
 
 
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/IOM/iom_nf90.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/IOM/iom_nf90.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/IOM/iom_nf90.F90	(revision 13193)
@@ -62,5 +62,7 @@
       CHARACTER(LEN=256) ::   clinfo           ! info character
       CHARACTER(LEN=256) ::   cltmp            ! temporary character
+      CHARACTER(LEN=12 ) ::   clfmt            ! writing format
       CHARACTER(LEN=3  ) ::   clcomp           ! name of component calling iom_nf90_open
+      INTEGER            ::   idg              ! number of digits
       INTEGER            ::   iln              ! lengths of character
       INTEGER            ::   istop            ! temporary storage of nstop
@@ -109,5 +111,7 @@
          IF( ldwrt ) THEN              !* the file should be open in write mode so we create it...
             IF( jpnij > 1 ) THEN
-               WRITE(cltmp,'(a,a,i4.4,a)') cdname(1:iln-1), '_', narea-1, '.nc'
+               idg = MAX( INT(LOG10(REAL(MAX(1,jpnij-1),wp))) + 1, 4 )          ! how many digits to we need to write? min=4, max=9
+               WRITE(clfmt, "('(a,a,i', i1, '.', i1, ',a)')") idg, idg          ! '(a,a,ix.x,a)'
+               WRITE(cltmp,clfmt) cdname(1:iln-1), '_', narea-1, '.nc'
                cdname = TRIM(cltmp)
             ENDIF
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ISF/isfdiags.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ISF/isfdiags.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ISF/isfdiags.F90	(revision 13193)
@@ -89,5 +89,5 @@
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: phtbl, pfrac  ! thickness of the tbl and fraction of last cell affected by the tbl
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pvar2d        ! 2d var to map in 3d
-      CHARACTER(LEN=256), INTENT(in) :: cdvar
+      CHARACTER(LEN=*), INTENT(in) :: cdvar
       !!---------------------------------------------------------------------
       INTEGER  :: ji, jj, jk                       ! loop indices
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/LBC/lib_mpp.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/LBC/lib_mpp.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/LBC/lib_mpp.F90	(revision 13193)
@@ -1112,10 +1112,20 @@
       CHARACTER(len=*), INTENT(in   ), OPTIONAL ::        cd2, cd3, cd4, cd5
       CHARACTER(len=*), INTENT(in   ), OPTIONAL ::   cd6, cd7, cd8, cd9, cd10
+      !
+      CHARACTER(LEN=8) ::   clfmt            ! writing format
+      INTEGER          ::   inum
       !!----------------------------------------------------------------------
       !
       nstop = nstop + 1
       !
-      ! force to open ocean.output file if not already opened
-      IF( numout == 6 ) CALL ctl_opn( numout, 'ocean.output', 'APPEND', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE. )
+      IF( cd1 == 'STOP' .AND. narea /= 1 ) THEN    ! Immediate stop: add an arror message in 'ocean.output' file
+         CALL ctl_opn( inum, 'ocean.output', 'APPEND', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE. )
+         WRITE(inum,*)
+         WRITE(inum,*) ' ==>>>   Look for "E R R O R" messages in all existing *ocean.output* files'
+         CLOSE(inum)
+      ENDIF
+      IF( numout == 6 ) THEN                       ! force to open ocean.output file if not already opened
+         CALL ctl_opn( numout, 'ocean.output', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, -1, .FALSE., narea )
+      ENDIF
       !
                             WRITE(numout,*)
@@ -1145,4 +1155,6 @@
          WRITE(numout,*)  'huge E-R-R-O-R : immediate stop'
          WRITE(numout,*)  
+         CALL FLUSH(numout)
+         CALL SLEEP(60)   ! make sure that all output and abort files are written by all cores. 60s should be enough...
          CALL mppstop( ld_abort = .true. )
       ENDIF
@@ -1207,5 +1219,7 @@
       !
       CHARACTER(len=80) ::   clfile
+      CHARACTER(LEN=10) ::   clfmt            ! writing format
       INTEGER           ::   iost
+      INTEGER           ::   idg              ! number of digits
       !!----------------------------------------------------------------------
       !
@@ -1214,5 +1228,10 @@
       clfile = TRIM(cdfile)
       IF( PRESENT( karea ) ) THEN
-         IF( karea > 1 )   WRITE(clfile, "(a,'_',i4.4)") TRIM(clfile), karea-1
+         IF( karea > 1 ) THEN
+            ! Warning: jpnij is maybe not already defined when calling ctl_opn -> use mppsize instead of jpnij
+            idg = MAX( INT(LOG10(REAL(MAX(1,mppsize-1),wp))) + 1, 4 )      ! how many digits to we need to write? min=4, max=9
+            WRITE(clfmt, "('(a,a,i', i1, '.', i1, ')')") idg, idg          ! '(a,a,ix.x)'
+            WRITE(clfile, clfmt) TRIM(clfile), '_', karea-1
+         ENDIF
       ENDIF
 #if defined key_agrif
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/LBC/mpp_loc_generic.h90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/LBC/mpp_loc_generic.h90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/LBC/mpp_loc_generic.h90	(revision 13193)
@@ -32,5 +32,4 @@
       REAL(wp)        , INTENT(  out) ::   pmin    ! Global minimum of ptab
       INDEX_TYPE(:)                                ! index of minimum in global frame
-# if defined key_mpp_mpi
       !
       INTEGER  ::   ierror, ii, idim
@@ -56,20 +55,20 @@
          !
          kindex(1) = mig( ilocs(1) )
-#  if defined DIM_2d || defined DIM_3d    /* avoid warning when kindex has 1 element */
+#if defined DIM_2d || defined DIM_3d    /* avoid warning when kindex has 1 element */
          kindex(2) = mjg( ilocs(2) )
-#  endif
-#  if defined DIM_3d                      /* avoid warning when kindex has 2 elements */
+#endif
+#if defined DIM_3d                      /* avoid warning when kindex has 2 elements */
          kindex(3) = ilocs(3)
-#  endif
+#endif
          ! 
          DEALLOCATE (ilocs)
          !
          index0 = kindex(1)-1   ! 1d index starting at 0
-#  if defined DIM_2d || defined DIM_3d   /* avoid warning when kindex has 1 element */
+#if defined DIM_2d || defined DIM_3d   /* avoid warning when kindex has 1 element */
          index0 = index0 + jpiglo * (kindex(2)-1)
-#  endif
-#  if defined DIM_3d                     /* avoid warning when kindex has 2 elements */
+#endif
+#if defined DIM_3d                     /* avoid warning when kindex has 2 elements */
          index0 = index0 + jpiglo * jpjglo * (kindex(3)-1)
-#  endif
+#endif
       END IF
       zain(1,:) = zmin
@@ -77,23 +76,23 @@
       !
       IF( ln_timing ) CALL tic_tac(.TRUE., ld_global = .TRUE.)
+#if defined key_mpp_mpi
       CALL MPI_ALLREDUCE( zain, zaout, 1, MPI_2DOUBLE_PRECISION, MPI_OPERATION ,MPI_COMM_OCE, ierror)
+#else
+      zaout(:,:) = zain(:,:)
+#endif
       IF( ln_timing ) CALL tic_tac(.FALSE., ld_global = .TRUE.)
       !
       pmin      = zaout(1,1)
       index0    = NINT( zaout(2,1) )
-#  if defined DIM_3d                     /* avoid warning when kindex has 2 elements */
+#if defined DIM_3d                     /* avoid warning when kindex has 2 elements */
       kindex(3) = index0 / (jpiglo*jpjglo)
       index0    = index0 - kindex(3) * (jpiglo*jpjglo)
-#  endif
-#  if defined DIM_2d || defined DIM_3d   /* avoid warning when kindex has 1 element */
+#endif
+#if defined DIM_2d || defined DIM_3d   /* avoid warning when kindex has 1 element */
       kindex(2) = index0 / jpiglo
       index0 = index0 - kindex(2) * jpiglo
-#  endif
+#endif
       kindex(1) = index0
       kindex(:) = kindex(:) + 1   ! start indices at 1
-#else
-      kindex = 0 ; pmin = 0.
-      WRITE(*,*) 'ROUTINE_LOC: You should not have seen this print! error?'
-#endif
 
    END SUBROUTINE ROUTINE_LOC
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/OBS/obs_grid.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/OBS/obs_grid.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/OBS/obs_grid.F90	(revision 13193)
@@ -684,5 +684,7 @@
          & fhistx1, fhistx2, fhisty1, fhisty2
       REAL(wp) :: histtol
-      
+      CHARACTER(LEN=26) :: clfmt            ! writing format
+      INTEGER           :: idg              ! number of digits
+ 
       IF (ln_grid_search_lookup) THEN
          
@@ -709,9 +711,10 @@
 
          IF ( ln_grid_global ) THEN
-            WRITE(cfname, FMT="(A,'_',A)") &
-               &          TRIM(cn_gridsearchfile), 'global.nc'
+            WRITE(cfname, FMT="(A,'_',A)") TRIM(cn_gridsearchfile), 'global.nc'
          ELSE
-            WRITE(cfname, FMT="(A,'_',I4.4,'of',I4.4,'by',I4.4,'.nc')") &
-               &          TRIM(cn_gridsearchfile), nproc, jpni, jpnj
+            idg = MAX( INT(LOG10(REAL(jpnij,wp))) + 1, 4 )        ! how many digits to we need to write? min=4, max=9
+            ! define the following format: "(a,a,ix.x,a,ix.x,a,ix.x,a)"
+            WRITE(clfmt, "('(a,a,i', i1, '.', i1',a,i', i1, '.', i1',a,i', i1, '.', i1',a)')") idg, idg, idg, idg, idg, idg
+            WRITE(cfname,      clfmt     ) TRIM(cn_gridsearchfile),'_', nproc,'of', jpni,'by', jpnj,'.nc'
          ENDIF
 
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/OBS/obs_write.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/OBS/obs_write.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/OBS/obs_write.F90	(revision 13193)
@@ -86,4 +86,6 @@
       CHARACTER(LEN=40) :: clfname
       CHARACTER(LEN=10) :: clfiletype
+      CHARACTER(LEN=12) :: clfmt            ! writing format
+      INTEGER :: idg                        ! number of digits
       INTEGER :: ilevel
       INTEGER :: jvar
@@ -181,5 +183,7 @@
       fbdata%caddname(1)   = 'Hx'
 
-      WRITE(clfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(clfiletype), nproc
+      idg = MAX( INT(LOG10(REAL(jpnij,wp))) + 1, 4 )            ! how many digits to we need to write? min=4, max=9
+      WRITE(clfmt, "('(a,a,i', i1, '.', i1, ',a)')") idg, idg   ! '(a,a,ix.x,a)'
+      WRITE(clfname,clfmt) TRIM(clfiletype), '_fdbk_', nproc, '.nc'
 
       IF(lwp) THEN
@@ -326,4 +330,6 @@
       CHARACTER(LEN=10) :: clfiletype
       CHARACTER(LEN=12), PARAMETER :: cpname = 'obs_wri_surf'
+      CHARACTER(LEN=12) :: clfmt           ! writing format
+      INTEGER :: idg                       ! number of digits
       INTEGER :: jo
       INTEGER :: ja
@@ -453,5 +459,7 @@
       fbdata%caddname(1)   = 'Hx'
 
-      WRITE(clfname, FMT="(A,'_fdbk_',I4.4,'.nc')") TRIM(clfiletype), nproc
+      idg = MAX( INT(LOG10(REAL(jpnij,wp))) + 1, 4 )            ! how many digits to we need to write? min=4, max=9
+      WRITE(clfmt, "('(a,a,i', i1, '.', i1, ',a)')") idg, idg   ! '(a,a,ix.x,a)'
+      WRITE(clfname,clfmt) TRIM(clfiletype), '_fdbk_', nproc, '.nc'
 
       IF(lwp) THEN
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC/sbcblk.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC/sbcblk.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC/sbcblk.F90	(revision 13193)
@@ -627,5 +627,16 @@
 
       END SELECT
-
+      
+      IF( iom_use('Cd_oce') )   CALL iom_put("Cd_oce",   zcd_oce * tmask(:,:,1))
+      IF( iom_use('Ce_oce') )   CALL iom_put("Ce_oce",   zce_oce * tmask(:,:,1))
+      IF( iom_use('Ch_oce') )   CALL iom_put("Ch_oce",   zch_oce * tmask(:,:,1))
+      !! LB: mainly here for debugging purpose:
+      IF( iom_use('theta_zt') ) CALL iom_put("theta_zt", (ztpot-rt0) * tmask(:,:,1)) ! potential temperature at z=zt
+      IF( iom_use('q_zt') )     CALL iom_put("q_zt",     zqair       * tmask(:,:,1)) ! specific humidity       "
+      IF( iom_use('theta_zu') ) CALL iom_put("theta_zu", (t_zu -rt0) * tmask(:,:,1)) ! potential temperature at z=zu
+      IF( iom_use('q_zu') )     CALL iom_put("q_zu",     q_zu        * tmask(:,:,1)) ! specific humidity       "
+      IF( iom_use('ssq') )      CALL iom_put("ssq",      pssq        * tmask(:,:,1)) ! saturation specific humidity at z=0
+      IF( iom_use('wspd_blk') ) CALL iom_put("wspd_blk", zU_zu       * tmask(:,:,1)) ! bulk wind speed at z=zu
+      
       IF( ln_skin_cs .OR. ln_skin_wl ) THEN
          !! ptsk and pssq have been updated!!!
@@ -643,11 +654,11 @@
 
       IF( ln_abl ) THEN         !==  ABL formulation  ==!   multiplication by rho_air and turbulent fluxes computation done in ablstp
-         !! FL do we need this multiplication by tmask ... ???
          DO_2D_11_11
-            zztmp = zU_zu(ji,jj) !* tmask(ji,jj,1)
+            zztmp = zU_zu(ji,jj)
             wndm(ji,jj)   = zztmp                   ! Store zU_zu in wndm to compute ustar2 in ablmod
             pcd_du(ji,jj) = zztmp * zcd_oce(ji,jj)
             psen(ji,jj)   = zztmp * zch_oce(ji,jj)
             pevp(ji,jj)   = zztmp * zce_oce(ji,jj)
+            rhoa(ji,jj)   = rho_air( ptair(ji,jj), phumi(ji,jj), pslp(ji,jj) )
          END_2D
       ELSE                      !==  BLK formulation  ==!   turbulent fluxes computation
@@ -673,6 +684,6 @@
          ! ... utau, vtau at U- and V_points, resp.
          !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
-         !     Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
-         DO_2D_10_10
+         !     Note that coastal wind stress is not used in the code... so this extra care has no effect
+         DO_2D_00_00
             utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj  ) ) &
                &          * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1))
@@ -878,24 +889,28 @@
          Ce_ice(:,:) = Ch_ice(:,:)       ! sensible and latent heat transfer coef. are considered identical
       ENDIF
-
-      !! IF ( iom_use("Cd_ice") ) CALL iom_put("Cd_ice", Cd_ice)   ! output value of pure ice-atm. transfer coef.
-      !! IF ( iom_use("Ch_ice") ) CALL iom_put("Ch_ice", Ch_ice)   ! output value of pure ice-atm. transfer coef.
-
+      
+      IF( iom_use('Cd_ice') ) CALL iom_put("Cd_ice", Cd_ice)
+      IF( iom_use('Ce_ice') ) CALL iom_put("Ce_ice", Ce_ice)
+      IF( iom_use('Ch_ice') ) CALL iom_put("Ch_ice", Ch_ice)
+      
       ! local scalars ( place there for vector optimisation purposes)
-      !IF (ln_abl) rhoa  (:,:)  = rho_air( ptair(:,:), phumi(:,:), pslp(:,:) ) !!GS: rhoa must be (re)computed here with ABL to avoid division by zero after (TBI)
       zcd_dui(:,:) = wndm_ice(:,:) * Cd_ice(:,:)
 
       IF( ln_blk ) THEN
-         ! ------------------------------------------------------------ !
-         !    Wind stress relative to the moving ice ( U10m - U_ice )   !
-         ! ------------------------------------------------------------ !
-         ! C-grid ice dynamics :   U & V-points (same as ocean)
-         DO_2D_00_00
-            putaui(ji,jj) = 0.5_wp * (  rhoa(ji+1,jj) * zcd_dui(ji+1,jj)             &
-               &                      + rhoa(ji  ,jj) * zcd_dui(ji  ,jj)  )          &
-               &         * ( 0.5_wp * ( pwndi(ji+1,jj) + pwndi(ji,jj) ) - rn_vfac * puice(ji,jj) )
-            pvtaui(ji,jj) = 0.5_wp * (  rhoa(ji,jj+1) * zcd_dui(ji,jj+1)             &
-               &                      + rhoa(ji,jj  ) * zcd_dui(ji,jj  )  )          &
-               &         * ( 0.5_wp * ( pwndj(ji,jj+1) + pwndj(ji,jj) ) - rn_vfac * pvice(ji,jj) )
+         ! ------------------------------------------------------------- !
+         !    Wind stress relative to the moving ice ( U10m - U_ice )    !
+         ! ------------------------------------------------------------- !
+         zztmp1 = rn_vfac * 0.5_wp
+         DO_2D_01_01    ! at T point 
+            putaui(ji,jj) = rhoa(ji,jj) * zcd_dui(ji,jj) * ( pwndi(ji,jj) - zztmp1 * ( puice(ji-1,jj  ) + puice(ji,jj) ) )
+            pvtaui(ji,jj) = rhoa(ji,jj) * zcd_dui(ji,jj) * ( pwndj(ji,jj) - zztmp1 * ( pvice(ji  ,jj-1) + pvice(ji,jj) ) )
+         END_2D
+         !
+         DO_2D_00_00    ! U & V-points (same as ocean).
+            ! take care of the land-sea mask to avoid "pollution" of coastal stress. p[uv]taui used in frazil and  rheology 
+            zztmp1 = 0.5_wp * ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji+1,jj  ,1) )
+            zztmp2 = 0.5_wp * ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji  ,jj+1,1) )
+            putaui(ji,jj) = zztmp1 * ( putaui(ji,jj) + putaui(ji+1,jj  ) )
+            pvtaui(ji,jj) = zztmp2 * ( pvtaui(ji,jj) + pvtaui(ji  ,jj+1) )
          END_2D
          CALL lbc_lnk_multi( 'sbcblk', putaui, 'U', -1., pvtaui, 'V', -1. )
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC/sbcblk_phy.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC/sbcblk_phy.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC/sbcblk_phy.F90	(revision 13193)
@@ -31,5 +31,5 @@
    REAL(wp), PARAMETER, PUBLIC :: R_vap   = 461.495_wp  !: Specific gas constant for water vapor          [J/K/kg]
    REAL(wp), PARAMETER, PUBLIC :: reps0   = R_dry/R_vap !: ratio of gas constant for dry air and water vapor => ~ 0.622
-   REAL(wp), PARAMETER, PUBLIC :: rctv0   = R_vap/R_dry !: for virtual temperature (== (1-eps)/eps) => ~ 0.608
+   REAL(wp), PARAMETER, PUBLIC :: rctv0   = R_vap/R_dry - 1._wp  !: for virtual temperature (== (1-eps)/eps) => ~ 0.608
    REAL(wp), PARAMETER, PUBLIC :: rCp_air = 1000.5_wp   !: specific heat of air (only used for ice fluxes now...)
    REAL(wp), PARAMETER, PUBLIC :: rCd_ice = 1.4e-3_wp   !: transfer coefficient over ice
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC/sbccpl.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC/sbccpl.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC/sbccpl.F90	(revision 13193)
@@ -365,5 +365,7 @@
       ! 
       ! Vectors: change of sign at north fold ONLY if on the local grid
-      IF( TRIM( sn_rcv_tau%cldes ) == 'oce only' .OR. TRIM(sn_rcv_tau%cldes ) == 'oce and ice') THEN ! avoid working with the atmospheric fields if they are not coupled
+      IF(       TRIM( sn_rcv_tau%cldes ) == 'oce only' .OR. TRIM( sn_rcv_tau%cldes ) == 'oce and ice'  &
+           .OR. TRIM( sn_rcv_tau%cldes ) == 'mixed oce-ice' ) THEN ! avoid working with the atmospheric fields if they are not coupled
+
       IF( TRIM( sn_rcv_tau%clvor ) == 'local grid' )   srcv(jpr_otx1:jpr_itz2)%nsgn = -1.
       
@@ -1482,4 +1484,5 @@
       INTEGER ::   ji, jj   ! dummy loop indices
       INTEGER ::   itx      ! index of taux over ice
+      REAL(wp)                     ::   zztmp1, zztmp2
       REAL(wp), DIMENSION(jpi,jpj) ::   ztx, zty 
       !!----------------------------------------------------------------------
@@ -1545,23 +1548,14 @@
             p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1)                   ! (U,V) ==> (U,V)
             p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1)
-         CASE( 'F' )
-            DO_2D_00_00
-               p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji,jj,1) + frcv(jpr_itx1)%z3(ji  ,jj-1,1) )
-               p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji,jj,1) + frcv(jpr_ity1)%z3(ji-1,jj  ,1) )
-            END_2D
          CASE( 'T' )
             DO_2D_00_00
-               p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji+1,jj  ,1) + frcv(jpr_itx1)%z3(ji,jj,1) )
-               p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji  ,jj+1,1) + frcv(jpr_ity1)%z3(ji,jj,1) )
+               ! take care of the land-sea mask to avoid "pollution" of coastal stress. p[uv]taui used in frazil and  rheology 
+               zztmp1 = 0.5_wp * ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji+1,jj  ,1) )
+               zztmp2 = 0.5_wp * ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji  ,jj+1,1) )
+               p_taui(ji,jj) = zztmp1 * ( frcv(jpr_itx1)%z3(ji+1,jj  ,1) + frcv(jpr_itx1)%z3(ji,jj,1) )
+               p_tauj(ji,jj) = zztmp2 * ( frcv(jpr_ity1)%z3(ji  ,jj+1,1) + frcv(jpr_ity1)%z3(ji,jj,1) )
             END_2D
-         CASE( 'I' )
-            DO_2D_00_00
-               p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji+1,jj+1,1) + frcv(jpr_itx1)%z3(ji+1,jj  ,1) )
-               p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji+1,jj+1,1) + frcv(jpr_ity1)%z3(ji  ,jj+1,1) )
-            END_2D
+            CALL lbc_lnk_multi( 'sbccpl', p_taui, 'U',  -1., p_tauj, 'V',  -1. )
          END SELECT
-         IF( srcv(jpr_itx1)%clgrid /= 'U' ) THEN 
-            CALL lbc_lnk_multi( 'sbccpl', p_taui, 'U',  -1., p_tauj, 'V',  -1. )
-         ENDIF
          
       ENDIF
@@ -1792,7 +1786,6 @@
             ENDDO
          ELSE
-            qns_tot(:,:) = qns_tot(:,:) + picefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
+            zqns_tot(:,:) = zqns_tot(:,:) + picefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
             DO jl = 1, jpl
-               zqns_tot(:,:   ) = zqns_tot(:,:) + picefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
                zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1)
             END DO
@@ -1935,7 +1928,6 @@
             END DO
          ELSE
-            qsr_tot(:,:   ) = qsr_tot(:,:) + picefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1)
+            zqsr_tot(:,:) = zqsr_tot(:,:) + picefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1)
             DO jl = 1, jpl
-               zqsr_tot(:,:   ) = zqsr_tot(:,:) + picefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1)
                zqsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1)
             END DO
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC/sbcmod.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC/sbcmod.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC/sbcmod.F90	(revision 13193)
@@ -120,15 +120,4 @@
       ncom_fsbc = nn_fsbc    ! make nn_fsbc available for lib_mpp
 #endif
-      !                             !* overwrite namelist parameter using CPP key information
-#if defined key_agrif
-      IF( Agrif_Root() ) THEN                ! AGRIF zoom (cf r1242: possibility to run without ice in fine grid)
-         IF( lk_si3  )   nn_ice      = 2
-         IF( lk_cice )   nn_ice      = 3
-      ENDIF
-!!GS: TBD
-!#else
-!      IF( lk_si3  )   nn_ice      = 2
-!      IF( lk_cice )   nn_ice      = 3
-#endif
       !
       IF(lwp) THEN                  !* Control print
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC/sbcwave.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC/sbcwave.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC/sbcwave.F90	(revision 13193)
@@ -212,13 +212,4 @@
       END_3D
       !
-#if defined key_agrif
-      IF( .NOT. Agrif_Root() ) THEN
-         IF( nbondi == -1 .OR. nbondi == 2 )   ze3divh( 2:nbghostcells+1,:        ,:) = 0._wp      ! west
-         IF( nbondi ==  1 .OR. nbondi == 2 )   ze3divh( nlci-nbghostcells:nlci-1,:,:) = 0._wp      ! east
-         IF( nbondj == -1 .OR. nbondj == 2 )   ze3divh( :,2:nbghostcells+1        ,:) = 0._wp      ! south
-         IF( nbondj ==  1 .OR. nbondj == 2 )   ze3divh( :,nlcj-nbghostcells:nlcj-1,:) = 0._wp      ! north
-      ENDIF
-#endif
-      !
       CALL lbc_lnk( 'sbcwave', ze3divh, 'T', 1. )
       !
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/STO/stopar.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/STO/stopar.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/STO/stopar.F90	(revision 13193)
@@ -684,10 +684,12 @@
       !! ** Purpose :   read stochastic parameters from restart file
       !!----------------------------------------------------------------------
-      INTEGER  :: jsto, jseed
+      INTEGER             ::   jsto, jseed
+      INTEGER             ::   idg                 ! number of digits
       INTEGER(KIND=8)     ::   ziseed(4)           ! RNG seeds in integer type
       REAL(KIND=8)        ::   zrseed(4)           ! RNG seeds in real type (with same bits to save in restart)
       CHARACTER(LEN=9)    ::   clsto2d='sto2d_000' ! stochastic parameter variable name
       CHARACTER(LEN=9)    ::   clsto3d='sto3d_000' ! stochastic parameter variable name
-      CHARACTER(LEN=10)   ::   clseed='seed0_0000' ! seed variable name
+      CHARACTER(LEN=15)   ::   clseed='seed0_0000' ! seed variable name
+      CHARACTER(LEN=6)    ::   clfmt               ! writing format
       !!----------------------------------------------------------------------
 
@@ -717,8 +719,10 @@
          IF (ln_rstseed) THEN
             ! Get saved state of the random number generator
+            idg = MAX( INT(LOG10(REAL(jpnij,wp))) + 1, 4 )        ! how many digits to we need to write? min=4, max=9
+            WRITE(clfmt, "('(i', i1, '.', i1, ')')") idg, idg     ! "(ix.x)"
             DO jseed = 1 , 4
-               WRITE(clseed(5:5) ,'(i1.1)') jseed
-               WRITE(clseed(7:10),'(i4.4)') narea
-               CALL iom_get( numstor, clseed , zrseed(jseed) )
+               WRITE(clseed(5:5)      ,'(i1.1)') jseed
+               WRITE(clseed(7:7+idg-1),  clfmt ) narea
+               CALL iom_get( numstor, clseed(1:7+idg-1) , zrseed(jseed) )
             END DO
             ziseed = TRANSFER( zrseed , ziseed)
@@ -742,5 +746,6 @@
       INTEGER, INTENT(in) ::   kt     ! ocean time-step
       !!
-      INTEGER  :: jsto, jseed
+      INTEGER             ::   jsto, jseed
+      INTEGER             ::   idg                 ! number of digits
       INTEGER(KIND=8)     ::   ziseed(4)           ! RNG seeds in integer type
       REAL(KIND=8)        ::   zrseed(4)           ! RNG seeds in real type (with same bits to save in restart)
@@ -749,5 +754,6 @@
       CHARACTER(LEN=9)    ::   clsto2d='sto2d_000' ! stochastic parameter variable name
       CHARACTER(LEN=9)    ::   clsto3d='sto3d_000' ! stochastic parameter variable name
-      CHARACTER(LEN=10)   ::   clseed='seed0_0000' ! seed variable name
+      CHARACTER(LEN=15)   ::   clseed='seed0_0000' ! seed variable name
+      CHARACTER(LEN=6)    ::   clfmt               ! writing format
       !!----------------------------------------------------------------------
 
@@ -771,8 +777,10 @@
             CALL kiss_state( ziseed(1) , ziseed(2) , ziseed(3) , ziseed(4) )
             zrseed = TRANSFER( ziseed , zrseed)
+            idg = MAX( INT(LOG10(REAL(jpnij,wp))) + 1, 4 )        ! how many digits to we need to write? min=4, max=9
+            WRITE(clfmt, "('(i', i1, '.', i1, ')')") idg, idg     ! "(ix.x)"
             DO jseed = 1 , 4
-               WRITE(clseed(5:5) ,'(i1.1)') jseed
-               WRITE(clseed(7:10),'(i4.4)') narea
-               CALL iom_rstput( kt, nitrst, numstow, clseed , zrseed(jseed) )
+               WRITE(clseed(5:5)      ,'(i1.1)') jseed
+               WRITE(clseed(7:7+idg-1),  clfmt ) narea
+               CALL iom_rstput( kt, nitrst, numstow, clseed(1:7+idg-1), zrseed(jseed) )
             END DO
             ! 2D stochastic parameters
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRD/trdtra.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRD/trdtra.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/TRD/trdtra.F90	(revision 13193)
@@ -83,5 +83,6 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL ::   ptra    ! now tracer variable
       !
-      INTEGER ::   jk   ! loop indices
+      INTEGER ::   jk    ! loop indices
+      INTEGER ::   i01   ! 0 or 1
       REAL(wp),        DIMENSION(jpi,jpj,jpk) ::   ztrds             ! 3D workspace
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zwt, zws, ztrdt   ! 3D workspace
@@ -91,8 +92,10 @@
          IF( trd_tra_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'trd_tra : unable to allocate arrays' )
       ENDIF
-
+      !
+      i01 = COUNT( (/ PRESENT(pu) .OR. ( ktrd /= jptra_xad .AND. ktrd /= jptra_yad .AND. ktrd /= jptra_zad ) /) )
+      !
       IF( ctype == 'TRA' .AND. ktra == jp_tem ) THEN   !==  Temperature trend  ==!
          !
-         SELECT CASE( ktrd )
+         SELECT CASE( ktrd*i01 )
          !                            ! advection: transform the advective flux into a trend
          CASE( jptra_xad )   ;   CALL trd_tra_adv( ptrd, pu, ptra, 'X', trdtx, Kmm ) 
@@ -113,5 +116,5 @@
       IF( ctype == 'TRA' .AND. ktra == jp_sal ) THEN      !==  Salinity trends  ==!
          !
-         SELECT CASE( ktrd )
+         SELECT CASE( ktrd*i01 )
          !                            ! advection: transform the advective flux into a trend
          !                            !            and send T & S trends to trd_tra_mng
@@ -168,5 +171,5 @@
       IF( ctype == 'TRC' ) THEN                           !==  passive tracer trend  ==!
          !
-         SELECT CASE( ktrd )
+         SELECT CASE( ktrd*i01 )
          !                            ! advection: transform the advective flux into a masked trend
          CASE( jptra_xad )   ;   CALL trd_tra_adv( ptrd , pu , ptra, 'X', ztrds, Kmm ) 
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/USR/usrdef_zgr.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/USR/usrdef_zgr.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/USR/usrdef_zgr.F90	(revision 13193)
@@ -202,5 +202,5 @@
       CALL lbc_lnk( 'usrdef_zgr', z2d, 'T', 1. )           ! set surrounding land to zero (here jperio=0 ==>> closed)
       !
-      k_bot(:,:) = NINT( z2d(:,:) )           ! =jpkm1 over the ocean point, =0 elsewhere
+      k_bot(:,:) = NINT( z2d(:,:) )          ! =jpkm1 over the ocean point, =0 elsewhere
       !
       k_top(:,:) = MIN( 1 , k_bot(:,:) )     ! = 1    over the ocean point, =0 elsewhere
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ZDF/zdftke.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ZDF/zdftke.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/ZDF/zdftke.F90	(revision 13193)
@@ -45,4 +45,10 @@
    USE zdfdrg         ! vertical physics: top/bottom drag coef.
    USE zdfmxl         ! vertical physics: mixed layer
+#if defined key_si3
+   USE ice, ONLY: hm_i, h_i
+#endif
+#if defined key_cice
+   USE sbc_ice, ONLY: h_i
+#endif
    !
    USE in_out_manager ! I/O manager
@@ -64,4 +70,6 @@
    INTEGER  ::   nn_mxl    ! type of mixing length (=0/1/2/3)
    REAL(wp) ::   rn_mxl0   ! surface  min value of mixing length (kappa*z_o=0.4*0.1 m)  [m]
+   INTEGER  ::      nn_mxlice ! type of scaling under sea-ice
+   REAL(wp) ::      rn_mxlice ! max constant ice thickness value when scaling under sea-ice ( nn_mxlice=1)
    INTEGER  ::   nn_pdl    ! Prandtl number or not (ratio avt/avm) (=0/1)
    REAL(wp) ::   rn_ediff  ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e)
@@ -245,5 +253,7 @@
                zetop = - 0.001875_wp * rCdU_top(ji,jj) * SQRT(  ( zmsku*( uu(ji,jj,mikt(ji,jj),Kbb)+uu(ji-1,jj,mikt(ji,jj),Kbb) ) )**2  &
                   &                                           + ( zmskv*( vv(ji,jj,mikt(ji,jj),Kbb)+vv(ji,jj-1,mikt(ji,jj),Kbb) ) )**2  )
-               en(ji,jj,mikt(ji,jj)) = en(ji,jj,1) * tmask(ji,jj,1) + MAX( zetop, rn_emin ) * (1._wp - tmask(ji,jj,1)) * ssmask(ji,jj)   ! masked at ocean surface
+               ! (1._wp - tmask(ji,jj,1)) * ssmask(ji,jj) = 1 where ice shelves are present
+               en(ji,jj,mikt(ji,jj)) = en(ji,jj,1)           * tmask(ji,jj,1) &
+                  &                  + MAX( zetop, rn_emin ) * (1._wp - tmask(ji,jj,1)) * ssmask(ji,jj)
             END_2D
          ENDIF
@@ -424,5 +434,5 @@
       REAL(wp) ::   zrn2, zraug, zcoef, zav   ! local scalars
       REAL(wp) ::   zdku,   zdkv, zsqen       !   -      -
-      REAL(wp) ::   zemxl, zemlm, zemlp       !   -      -
+      REAL(wp) ::   zemxl, zemlm, zemlp, zmaxice       !   -      -
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zmxlm, zmxld   ! 3D workspace
       !!--------------------------------------------------------------------
@@ -438,12 +448,51 @@
       zmxld(:,:,:)  = rmxl_min
       !
-      IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rho0*g)
+     IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rho0*g)
+         !
          zraug = vkarmn * 2.e5_wp / ( rho0 * grav )
+#if ! defined key_si3 && ! defined key_cice
          DO_2D_00_00
-            zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) )
+            zmxlm(ji,jj,1) =  zraug * taum(ji,jj) * tmask(ji,jj,1)
          END_2D
-      ELSE 
+#else
+         SELECT CASE( nn_mxlice )             ! Type of scaling under sea-ice
+         !
+         CASE( 0 )                      ! No scaling under sea-ice
+            DO_2D_00_00
+               zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1)
+            END_2D
+            !
+         CASE( 1 )                           ! scaling with constant sea-ice thickness
+            DO_2D_00_00
+               zmxlm(ji,jj,1) =  ( ( 1. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * rn_mxlice ) * tmask(ji,jj,1)
+            END_2D
+            !
+         CASE( 2 )                                 ! scaling with mean sea-ice thickness
+            DO_2D_00_00
+#if defined key_si3
+               zmxlm(ji,jj,1) = ( ( 1. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * hm_i(ji,jj) * 2. ) * tmask(ji,jj,1)
+#elif defined key_cice
+               zmaxice = MAXVAL( h_i(ji,jj,:) )
+               zmxlm(ji,jj,1) = ( ( 1. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1)
+#endif
+            END_2D
+            !
+         CASE( 3 )                                 ! scaling with max sea-ice thickness
+            DO_2D_00_00
+               zmaxice = MAXVAL( h_i(ji,jj,:) )
+               zmxlm(ji,jj,1) = ( ( 1. - fr_i(ji,jj) ) * zraug * taum(ji,jj) + fr_i(ji,jj) * zmaxice ) * tmask(ji,jj,1)
+            END_2D
+            !
+         END SELECT
+#endif
+         !
+         DO_2D_00_00
+            zmxlm(ji,jj,1) = MAX( rn_mxl0, zmxlm(ji,jj,1) )
+         END_2D
+         !
+      ELSE
          zmxlm(:,:,1) = rn_mxl0
       ENDIF
+
       !
       DO_3D_00_00( 2, jpkm1 )
@@ -554,8 +603,9 @@
       INTEGER             ::   ios
       !!
-      NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin  ,          &
-         &                 rn_emin0, rn_bshear, nn_mxl , ln_mxl0  ,          &
-         &                 rn_mxl0 , nn_pdl   , ln_drg , ln_lc    , rn_lc,   &
-         &                 nn_etau , nn_htau  , rn_efr , rn_eice  
+      NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb   , rn_emin  ,  &
+         &                 rn_emin0, rn_bshear, nn_mxl   , ln_mxl0  ,  &
+         &                 rn_mxl0 , nn_mxlice, rn_mxlice,             &
+         &                 nn_pdl  , ln_drg   , ln_lc    , rn_lc,      &
+         &                 nn_etau , nn_htau  , rn_efr   , rn_eice  
       !!----------------------------------------------------------------------
       !
@@ -583,4 +633,9 @@
          WRITE(numout,*) '      mixing length type                          nn_mxl    = ', nn_mxl
          WRITE(numout,*) '         surface mixing length = F(stress) or not    ln_mxl0   = ', ln_mxl0
+         IF( ln_mxl0 ) THEN
+            WRITE(numout,*) '      type of scaling under sea-ice               nn_mxlice = ', nn_mxlice
+            IF( nn_mxlice == 1 ) &
+            WRITE(numout,*) '      ice thickness when scaling under sea-ice    rn_mxlice = ', rn_mxlice
+         ENDIF         
          WRITE(numout,*) '         surface  mixing length minimum value        rn_mxl0   = ', rn_mxl0
          WRITE(numout,*) '      top/bottom friction forcing flag            ln_drg    = ', ln_drg
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/nemogcm.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/nemogcm.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/nemogcm.F90	(revision 13193)
@@ -194,11 +194,4 @@
       END DO
       !
-      IF( .NOT. Agrif_Root() ) THEN
-         CALL Agrif_ParentGrid_To_ChildGrid()
-         IF( ln_diaobs )   CALL dia_obs_wri
-         IF( ln_timing )   CALL timing_finalize
-         CALL Agrif_ChildGrid_To_ParentGrid()
-      ENDIF
-      !
 # else
       !
@@ -249,5 +242,13 @@
       IF( nstop /= 0 .AND. lwp ) THEN        ! error print
          WRITE(ctmp1,*) '   ==>>>   nemo_gcm: a total of ', nstop, ' errors have been found'
-         CALL ctl_stop( ctmp1 )
+         IF( ngrdstop > 0 ) THEN
+            WRITE(ctmp9,'(i2)') ngrdstop
+            WRITE(ctmp2,*) '           E R R O R detected in Agrif grid '//TRIM(ctmp9)
+            WRITE(ctmp3,*) '           Look for "E R R O R" messages in all existing '//TRIM(ctmp9)//'_ocean_output* files'
+            CALL ctl_stop( ' ', ctmp1, ' ', ctmp2, ' ', ctmp3 )
+         ELSE
+            WRITE(ctmp2,*) '           Look for "E R R O R" messages in all existing ocean_output* files'
+            CALL ctl_stop( ' ', ctmp1, ' ', ctmp2 )
+         ENDIF
       ENDIF
       !
@@ -261,5 +262,5 @@
 #else
       IF    ( lk_oasis ) THEN   ;   CALL cpl_finalize   ! end coupling and mpp communications with OASIS
-      ELSEIF( lk_mpp   ) THEN   ;   CALL mppstop      ! end mpp communications
+      ELSEIF( lk_mpp   ) THEN   ;   CALL mppstop        ! end mpp communications
       ENDIF
 #endif
@@ -347,20 +348,6 @@
       !
       ! finalize the definition of namctl variables
-      IF( sn_cfctl%l_allon ) THEN
-         ! Turn on all options.
-         CALL nemo_set_cfctl( sn_cfctl, .TRUE., .TRUE. )
-         ! Ensure all processors are active
-         sn_cfctl%procmin = 0 ; sn_cfctl%procmax = 1000000 ; sn_cfctl%procincr = 1
-      ELSEIF( sn_cfctl%l_config ) THEN
-         ! Activate finer control of report outputs
-         ! optionally switch off output from selected areas (note this only
-         ! applies to output which does not involve global communications)
-         IF( ( narea < sn_cfctl%procmin .OR. narea > sn_cfctl%procmax  ) .OR. &
-           & ( MOD( narea - sn_cfctl%procmin, sn_cfctl%procincr ) /= 0 ) )    &
-           &   CALL nemo_set_cfctl( sn_cfctl, .FALSE., .FALSE. )
-      ELSE
-         ! turn off all options.
-         CALL nemo_set_cfctl( sn_cfctl, .FALSE., .TRUE. )
-      ENDIF
+      IF( narea < sn_cfctl%procmin .OR. narea > sn_cfctl%procmax .OR. MOD( narea - sn_cfctl%procmin, sn_cfctl%procincr ) /= 0 )   &
+         &   CALL nemo_set_cfctl( sn_cfctl, .FALSE. )
       !
       lwp = (narea == 1) .OR. sn_cfctl%l_oceout    ! control of all listing output print
@@ -546,7 +533,4 @@
          WRITE(numout,*) '~~~~~~~~'
          WRITE(numout,*) '   Namelist namctl'
-         WRITE(numout,*) '                              sn_cfctl%l_glochk  = ', sn_cfctl%l_glochk
-         WRITE(numout,*) '                              sn_cfctl%l_allon   = ', sn_cfctl%l_allon
-         WRITE(numout,*) '       finer control over o/p sn_cfctl%l_config  = ', sn_cfctl%l_config
          WRITE(numout,*) '                              sn_cfctl%l_runstat = ', sn_cfctl%l_runstat
          WRITE(numout,*) '                              sn_cfctl%l_trcstat = ', sn_cfctl%l_trcstat
@@ -696,25 +680,18 @@
 
 
-   SUBROUTINE nemo_set_cfctl(sn_cfctl, setto, for_all )
+   SUBROUTINE nemo_set_cfctl(sn_cfctl, setto )
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE nemo_set_cfctl  ***
       !!
       !! ** Purpose :   Set elements of the output control structure to setto.
-      !!                for_all should be .false. unless all areas are to be
-      !!                treated identically.
       !!
       !! ** Method  :   Note this routine can be used to switch on/off some
-      !!                types of output for selected areas but any output types
-      !!                that involve global communications (e.g. mpp_max, glob_sum)
-      !!                should be protected from selective switching by the
-      !!                for_all argument
-      !!----------------------------------------------------------------------
-      LOGICAL :: setto, for_all
-      TYPE(sn_ctl) :: sn_cfctl
-      !!----------------------------------------------------------------------
-      IF( for_all ) THEN
-         sn_cfctl%l_runstat = setto
-         sn_cfctl%l_trcstat = setto
-      ENDIF
+      !!                types of output for selected areas.
+      !!----------------------------------------------------------------------
+      TYPE(sn_ctl), INTENT(inout) :: sn_cfctl
+      LOGICAL     , INTENT(in   ) :: setto
+      !!----------------------------------------------------------------------
+      sn_cfctl%l_runstat = setto
+      sn_cfctl%l_trcstat = setto
       sn_cfctl%l_oceout  = setto
       sn_cfctl%l_layout  = setto
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/step.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/step.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/step.F90	(revision 13193)
@@ -86,10 +86,9 @@
       !!----------------------------------------------------------------------
       INTEGER ::   ji, jj, jk   ! dummy loop indice
-      INTEGER ::   indic        ! error indicator if < 0
 !!gm kcall can be removed, I guess
       INTEGER ::   kcall        ! optional integer argument (dom_vvl_sf_nxt)
       !! ---------------------------------------------------------------------
 #if defined key_agrif
-      IF( nstop > 0 ) return   ! avoid to go further if an error was detected during previous time step 
+      IF( nstop > 0 ) RETURN   ! avoid to go further if an error was detected during previous time step (child grid)
       kstp = nit000 + Agrif_Nb_Step()
       Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs   ! agrif_oce module copies of time level indices
@@ -119,8 +118,6 @@
       ! update I/O and calendar 
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                             indic = 0                ! reset to no error condition
-                             
       IF( kstp == nit000 ) THEN                       ! initialize IOM context (must be done after nemo_init for AGRIF+XIOS+OASIS)
-                             CALL iom_init( cxios_context, ld_closedef=.FALSE. )   ! for model grid (including passible AGRIF zoom)
+                             CALL iom_init( cxios_context, ld_closedef=.FALSE. )   ! for model grid (including possible AGRIF zoom)
          IF( lk_diamlr   )   CALL dia_mlr_iom_init    ! with additional setup for multiple-linear-regression analysis
                              CALL iom_init_closedef
@@ -318,9 +315,11 @@
                          Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs      ! agrif_oce module copies of time level indices
                          CALL Agrif_Integrate_ChildGrids( stp )       ! allows to finish all the Child Grids before updating
+
 #endif
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Control
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                         CALL stp_ctl      ( kstp, Nbb, Nnn, indic )
+                         CALL stp_ctl      ( kstp, Nnn )
+
 #if defined key_agrif
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -330,6 +329,7 @@
                          CALL Agrif_update_all( )                  ! Update all components
       ENDIF
-#endif
-      IF( ln_diaobs  )   CALL dia_obs      ( kstp, Nnn )      ! obs-minus-model (assimilation) diagnostics (call after dynamics update)
+
+#endif
+      IF( ln_diaobs .AND. nstop == 0 )  CALL dia_obs( kstp, Nnn )  ! obs-minus-model (assimilation) diags (after dynamics update)
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -345,6 +345,5 @@
       ! Coupled mode
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-!!gm why lk_oasis and not lk_cpl ????
-      IF( lk_oasis   )   CALL sbc_cpl_snd( kstp, Nbb, Nnn )        ! coupled mode : field exchanges
+      IF( lk_oasis .AND. nstop == 0 )   CALL sbc_cpl_snd( kstp, Nbb, Nnn )     ! coupled mode : field exchanges
       !
 #if defined key_iomput
@@ -352,5 +351,5 @@
       ! Finalize contextes if end of simulation or error detected
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<                         
-      IF( kstp == nitend .OR. indic < 0 ) THEN 
+      IF( kstp == nitend .OR. nstop > 0 ) THEN 
                       CALL iom_context_finalize(      cxios_context          ) ! needed for XIOS+AGRIF
          IF( lrxios ) CALL iom_context_finalize(      crxios_context         )
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/stpctl.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/stpctl.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/stpctl.F90	(revision 13193)
@@ -19,12 +19,12 @@
    USE dom_oce         ! ocean space and time domain variables 
    USE c1d             ! 1D vertical configuration
+   USE zdf_oce ,  ONLY : ln_zad_Aimp       ! ocean vertical physics variables
+   USE wet_dry,   ONLY : ll_wd, ssh_ref    ! reference depth for negative bathy
+   !  
    USE diawri          ! Standard run outputs       (dia_wri_state routine)
-   !
    USE in_out_manager  ! I/O manager
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE lib_mpp         ! distributed memory computing
-   USE zdf_oce ,  ONLY : ln_zad_Aimp       ! ocean vertical physics variables
-   USE wet_dry,   ONLY : ll_wd, ssh_ref    ! reference depth for negative bathy
-
+   !
    USE netcdf          ! NetCDF library
    IMPLICIT NONE
@@ -33,6 +33,6 @@
    PUBLIC stp_ctl           ! routine called by step.F90
 
-   INTEGER  ::   idrun, idtime, idssh, idu, ids1, ids2, idt1, idt2, idc1, idw1, istatus
-   LOGICAL  ::   lsomeoce
+   INTEGER                ::   nrunid   ! netcdf file id
+   INTEGER, DIMENSION(8)  ::   nvarid   ! netcdf variable id
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -42,5 +42,5 @@
 CONTAINS
 
-   SUBROUTINE stp_ctl( kt, Kbb, Kmm, kindic )
+   SUBROUTINE stp_ctl( kt, Kmm )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE stp_ctl  ***
@@ -50,5 +50,5 @@
       !! ** Method  : - Save the time step in numstp
       !!              - Print it each 50 time steps
-      !!              - Stop the run IF problem encountered by setting indic=-3
+      !!              - Stop the run IF problem encountered by setting nstop > 0
       !!                Problems checked: |ssh| maximum larger than 10 m
       !!                                  |U|   maximum larger than 10 m/s 
@@ -57,143 +57,251 @@
       !! ** Actions :   "time.step" file = last ocean time-step
       !!                "run.stat"  file = run statistics
-      !!                nstop indicator sheared among all local domain (lk_mpp=T)
+      !!                 nstop indicator sheared among all local domain
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in   ) ::   kt       ! ocean time-step index
-      INTEGER, INTENT(in   ) ::   Kbb, Kmm      ! ocean time level index
-      INTEGER, INTENT(inout) ::   kindic   ! error indicator
-      !!
-      INTEGER                ::   ji, jj, jk          ! dummy loop indices
-      INTEGER, DIMENSION(2)  ::   ih                  ! min/max loc indices
-      INTEGER, DIMENSION(3)  ::   iu, is1, is2        ! min/max loc indices
-      REAL(wp)               ::   zzz                 ! local real 
-      REAL(wp), DIMENSION(9) ::   zmax
-      LOGICAL                ::   ll_wrtstp, ll_colruns, ll_wrtruns
-      CHARACTER(len=20) :: clname
-      !!----------------------------------------------------------------------
-      !
-      ll_wrtstp  = ( MOD( kt, sn_cfctl%ptimincr ) == 0 ) .OR. ( kt == nitend )
-      ll_colruns = ll_wrtstp .AND. ( sn_cfctl%l_runstat )
-      ll_wrtruns = ll_colruns .AND. lwm
-      IF( kt == nit000 .AND. lwp ) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'stp_ctl : time-stepping control'
-         WRITE(numout,*) '~~~~~~~'
-         !                                ! open time.step file
-         IF( lwm ) CALL ctl_opn( numstp, 'time.step', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
-         !                                ! open run.stat file(s) at start whatever
-         !                                ! the value of sn_cfctl%ptimincr
-         IF( lwm .AND. ( sn_cfctl%l_runstat ) ) THEN
+      INTEGER, INTENT(in   ) ::   Kmm      ! ocean time level index
+      !!
+      INTEGER                         ::   ji                                    ! dummy loop indices
+      INTEGER                         ::   idtime, istatus
+      INTEGER , DIMENSION(9)          ::   iareasum, iareamin, iareamax
+      INTEGER , DIMENSION(3,4)        ::   iloc                                  ! min/max loc indices
+      REAL(wp)                        ::   zzz                                   ! local real 
+      REAL(wp), DIMENSION(9)          ::   zmax, zmaxlocal
+      LOGICAL                         ::   ll_wrtstp, ll_colruns, ll_wrtruns
+      LOGICAL, DIMENSION(jpi,jpj,jpk) ::   llmsk
+      CHARACTER(len=20)               ::   clname
+      !!----------------------------------------------------------------------
+      IF( nstop > 0 .AND. ngrdstop > -1 )   RETURN   !   stpctl was already called by a child grid
+      !
+      ll_wrtstp  = ( MOD( kt-nit000, sn_cfctl%ptimincr ) == 0 ) .OR. ( kt == nitend )
+      ll_colruns = ll_wrtstp .AND. sn_cfctl%l_runstat .AND. jpnij > 1 
+      ll_wrtruns = ( ll_colruns .OR. jpnij == 1 ) .AND. lwm
+      !
+      IF( kt == nit000 ) THEN
+         !
+         IF( lwp ) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) 'stp_ctl : time-stepping control'
+            WRITE(numout,*) '~~~~~~~'
+         ENDIF
+         !                                ! open time.step    ascii file, done only by 1st subdomain
+         IF( lwm )   CALL ctl_opn( numstp, 'time.step', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
+         !
+         IF( ll_wrtruns ) THEN
+            !                             ! open run.stat     ascii file, done only by 1st subdomain
             CALL ctl_opn( numrun, 'run.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
+            !                             ! open run.stat.nc netcdf file, done only by 1st subdomain
             clname = 'run.stat.nc'
             IF( .NOT. Agrif_Root() )   clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname)
-            istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, idrun )
-            istatus = NF90_DEF_DIM( idrun, 'time', NF90_UNLIMITED, idtime )
-            istatus = NF90_DEF_VAR( idrun, 'abs_ssh_max', NF90_DOUBLE, (/ idtime /), idssh )
-            istatus = NF90_DEF_VAR( idrun,   'abs_u_max', NF90_DOUBLE, (/ idtime /), idu   )
-            istatus = NF90_DEF_VAR( idrun,       's_min', NF90_DOUBLE, (/ idtime /), ids1  )
-            istatus = NF90_DEF_VAR( idrun,       's_max', NF90_DOUBLE, (/ idtime /), ids2  )
-            istatus = NF90_DEF_VAR( idrun,       't_min', NF90_DOUBLE, (/ idtime /), idt1  )
-            istatus = NF90_DEF_VAR( idrun,       't_max', NF90_DOUBLE, (/ idtime /), idt2  )
+            istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, nrunid )
+            istatus = NF90_DEF_DIM( nrunid, 'time', NF90_UNLIMITED, idtime )
+            istatus = NF90_DEF_VAR( nrunid, 'abs_ssh_max', NF90_DOUBLE, (/ idtime /), nvarid(1) )
+            istatus = NF90_DEF_VAR( nrunid,   'abs_u_max', NF90_DOUBLE, (/ idtime /), nvarid(2) )
+            istatus = NF90_DEF_VAR( nrunid,       's_min', NF90_DOUBLE, (/ idtime /), nvarid(3) )
+            istatus = NF90_DEF_VAR( nrunid,       's_max', NF90_DOUBLE, (/ idtime /), nvarid(4) )
+            istatus = NF90_DEF_VAR( nrunid,       't_min', NF90_DOUBLE, (/ idtime /), nvarid(5) )
+            istatus = NF90_DEF_VAR( nrunid,       't_max', NF90_DOUBLE, (/ idtime /), nvarid(6) )
             IF( ln_zad_Aimp ) THEN
-               istatus = NF90_DEF_VAR( idrun,   'abs_wi_max', NF90_DOUBLE, (/ idtime /), idw1  )
-               istatus = NF90_DEF_VAR( idrun,       'Cf_max', NF90_DOUBLE, (/ idtime /), idc1  )
+               istatus = NF90_DEF_VAR( nrunid,   'Cf_max', NF90_DOUBLE, (/ idtime /), nvarid(7) )
+               istatus = NF90_DEF_VAR( nrunid,'abs_wi_max',NF90_DOUBLE, (/ idtime /), nvarid(8) )
             ENDIF
-            istatus = NF90_ENDDEF(idrun)
-            zmax(8:9) = 0._wp    ! initialise to zero in case ln_zad_Aimp option is not in use
-         ENDIF
-      ENDIF
-      IF( kt == nit000 )   lsomeoce = COUNT( ssmask(:,:) == 1._wp ) > 0
-      !
-      IF(lwm .AND. ll_wrtstp) THEN        !==  current time step  ==!   ("time.step" file)
+            istatus = NF90_ENDDEF(nrunid)
+         ENDIF
+         !    
+      ENDIF
+      !
+      !                                   !==              write current time step              ==!
+      !                                   !==  done only by 1st subdomain at writting timestep  ==!
+      IF( lwm .AND. ll_wrtstp ) THEN
          WRITE ( numstp, '(1x, i8)' )   kt
          REWIND( numstp )
       ENDIF
-      !
-      !                                   !==  test of extrema  ==!
-      IF( ll_wd ) THEN
-         zmax(1) = MAXVAL(  ABS( ssh(:,:,Kmm) + ssh_ref*tmask(:,:,1) )  )        ! ssh max 
-      ELSE
-         zmax(1) = MAXVAL(  ABS( ssh(:,:,Kmm) )  )                               ! ssh max
-      ENDIF
-      zmax(2) = MAXVAL(  ABS( uu(:,:,:,Kmm) )  )                                  ! velocity max (zonal only)
-      zmax(3) = MAXVAL( -ts(:,:,:,jp_sal,Kmm) , mask = tmask(:,:,:) == 1._wp )   ! minus salinity max
-      zmax(4) = MAXVAL(  ts(:,:,:,jp_sal,Kmm) , mask = tmask(:,:,:) == 1._wp )   !       salinity max
-      zmax(5) = MAXVAL( -ts(:,:,:,jp_tem,Kmm) , mask = tmask(:,:,:) == 1._wp )   ! minus temperature max
-      zmax(6) = MAXVAL(  ts(:,:,:,jp_tem,Kmm) , mask = tmask(:,:,:) == 1._wp )   !       temperature max
-      zmax(7) = REAL( nstop , wp )                                            ! stop indicator
-      IF( ln_zad_Aimp ) THEN
-         zmax(8) = MAXVAL(  ABS( wi(:,:,:) ) , mask = wmask(:,:,:) == 1._wp ) ! implicit vertical vel. max
-         zmax(9) = MAXVAL(   Cu_adv(:,:,:)   , mask = tmask(:,:,:) == 1._wp ) ! partitioning coeff. max
-      ENDIF
-      !
+      !                                   !==            test of local extrema           ==!
+      !                                   !==  done by all processes at every time step  ==!
+      !
+      ! define zmax default value. needed for land processors
+      IF( ll_colruns ) THEN    ! default value: must not be kept when calling mpp_max -> must be as small as possible
+         zmax(:) = -HUGE(1._wp)
+      ELSE                     ! default value: must not give true for any of the tests bellow (-> avoid manipulating HUGE...)
+         zmax(:) =  0._wp
+         zmax(3) = -1._wp      ! avoid salinity minimum at 0.
+      ENDIF
+      !
+      llmsk(:,:,1) = ssmask(:,:) == 1._wp
+      IF( COUNT( llmsk(:,:,1) ) > 0 ) THEN   ! avoid huge values sent back for land processors...
+         IF( ll_wd ) THEN
+            zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm) + ssh_ref ), mask = llmsk(:,:,1) )   ! ssh max
+         ELSE
+            zmax(1) = MAXVAL( ABS( ssh(:,:,Kmm)           ), mask = llmsk(:,:,1) )   ! ssh max
+         ENDIF
+      ENDIF
+      zmax(2) = MAXVAL( ABS( uu(:,:,:,Kmm) ) )                                       ! velocity max (zonal only)
+      llmsk(:,:,:) = tmask(:,:,:) == 1._wp
+      IF( COUNT( llmsk(:,:,:) ) > 0 ) THEN   ! avoid huge values sent back for land processors...
+         zmax(3) = MAXVAL( -ts(:,:,:,jp_sal,Kmm), mask = llmsk )                     ! minus salinity max
+         zmax(4) = MAXVAL(  ts(:,:,:,jp_sal,Kmm), mask = llmsk )                     !       salinity max
+         IF( ll_colruns .OR. jpnij == 1 ) THEN     ! following variables are used only in the netcdf file
+            zmax(5) = MAXVAL( -ts(:,:,:,jp_tem,Kmm), mask = llmsk )                  ! minus temperature max
+            zmax(6) = MAXVAL(  ts(:,:,:,jp_tem,Kmm), mask = llmsk )                  !       temperature max
+            IF( ln_zad_Aimp ) THEN
+               zmax(7) = MAXVAL(   Cu_adv(:,:,:)   , mask = llmsk )                  ! partitioning coeff. max
+               llmsk(:,:,:) = wmask(:,:,:) == 1._wp
+               IF( COUNT( llmsk(:,:,:) ) > 0 ) THEN   ! avoid huge values sent back for land processors...
+                  zmax(8) = MAXVAL(ABS( wi(:,:,:) ), mask = llmsk )                  ! implicit vertical vel. max
+               ENDIF
+            ENDIF
+         ENDIF
+      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
-         nstop = NINT( zmax(7) )                 ! nstop indicator sheared among all local domains
-      ENDIF
-      !                                   !==  run statistics  ==!   ("run.stat" files)
+         nstop = NINT( zmax(9) )                 ! update nstop indicator (now sheared among all local domains)
+      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( idrun, idssh, (/ zmax(1)/), (/kt/), (/1/) )
-         istatus = NF90_PUT_VAR( idrun,   idu, (/ zmax(2)/), (/kt/), (/1/) )
-         istatus = NF90_PUT_VAR( idrun,  ids1, (/-zmax(3)/), (/kt/), (/1/) )
-         istatus = NF90_PUT_VAR( idrun,  ids2, (/ zmax(4)/), (/kt/), (/1/) )
-         istatus = NF90_PUT_VAR( idrun,  idt1, (/-zmax(5)/), (/kt/), (/1/) )
-         istatus = NF90_PUT_VAR( idrun,  idt2, (/ zmax(6)/), (/kt/), (/1/) )
+         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( idrun,  idw1, (/ zmax(8)/), (/kt/), (/1/) )
-            istatus = NF90_PUT_VAR( idrun,  idc1, (/ zmax(9)/), (/kt/), (/1/) )
-         ENDIF
-         IF( MOD( kt , 100 ) == 0 ) istatus = NF90_SYNC(idrun)
-         IF( kt == nitend         ) istatus = NF90_CLOSE(idrun)
+            istatus = NF90_PUT_VAR( nrunid, nvarid(7), (/ zmax(7)/), (/kt/), (/1/) )
+            istatus = NF90_PUT_VAR( nrunid, nvarid(8), (/ zmax(8)/), (/kt/), (/1/) )
+         ENDIF
+         IF( kt == nitend )   istatus = NF90_CLOSE(nrunid)
       END IF
-      !                                   !==  error handling  ==!
-      IF( ( sn_cfctl%l_glochk .OR. lsomeoce ) .AND. (   &  ! domain contains some ocean points, check for sensible ranges
-         &  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) ) ) ) THEN   ! NaN encounter in the tests
-         IF( lk_mpp .AND. sn_cfctl%l_glochk ) THEN
-            ! have use mpp_max (because sn_cfctl%l_glochk=.T. and distributed)
-            CALL mpp_maxloc( 'stpctl', ABS(ssh(:,:,Kmm))        , ssmask(:,:)  , zzz, ih  )
-            CALL mpp_maxloc( 'stpctl', ABS(uu(:,:,:,Kmm))          , umask (:,:,:), zzz, iu  )
-            CALL mpp_minloc( 'stpctl', ts(:,:,:,jp_sal,Kmm), tmask (:,:,:), zzz, is1 )
-            CALL mpp_maxloc( 'stpctl', ts(:,:,:,jp_sal,Kmm), tmask (:,:,:), zzz, is2 )
+      !                                   !==               error handling               ==!
+      !                                   !==  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
+         !
+         iloc(:,:) = 0
+         IF( ll_colruns ) THEN   ! zmax is global, so it is the same on all subdomains -> no dead lock with mpp_maxloc
+            ! first: close the netcdf file, so we can read it
+            IF( lwm .AND. kt /= nitend )   istatus = NF90_CLOSE(nrunid)
+            ! get global loc on the min/max
+            CALL mpp_maxloc( 'stpctl', ABS(ssh(:,:,         Kmm)), ssmask(:,:  ), zzz, iloc(1:2,1) )   ! mpp_maxloc ok if mask = F 
+            CALL mpp_maxloc( 'stpctl', ABS( uu(:,:,:,       Kmm)),  umask(:,:,:), zzz, iloc(1:3,2) )
+            CALL mpp_minloc( 'stpctl',      ts(:,:,:,jp_sal,Kmm) ,  tmask(:,:,:), zzz, iloc(1:3,3) )
+            CALL mpp_maxloc( 'stpctl',      ts(:,:,:,jp_sal,Kmm) ,  tmask(:,:,:), zzz, iloc(1:3,4) )
+            ! find which subdomain has the max.
+            iareamin(:) = jpnij+1   ;   iareamax(:) = 0   ;   iareasum(:) = 0
+            DO ji = 1, 9
+               IF( zmaxlocal(ji) == zmax(ji) ) THEN
+                  iareamin(ji) = narea   ;   iareamax(ji) = narea   ;   iareasum(ji) = 1
+               ENDIF
+            END DO
+            CALL mpp_min( "stpctl", iareamin )         ! min over the global domain
+            CALL mpp_max( "stpctl", iareamax )         ! max over the global domain
+            CALL mpp_sum( "stpctl", iareasum )         ! sum over the global domain
+         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
+            iloc(1:2,1) = MAXLOC( ABS( ssh(:,:,         Kmm)), mask = ssmask(:,:  ) == 1._wp ) + (/ nimpp - 1, njmpp - 1    /)
+            iloc(1:3,2) = MAXLOC( ABS(  uu(:,:,:,       Kmm)), mask =  umask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /)
+            iloc(1:3,3) = MINLOC(       ts(:,:,:,jp_sal,Kmm) , mask =  tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /)
+            iloc(1:3,4) = MAXLOC(       ts(:,:,:,jp_sal,Kmm) , mask =  tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /)
+            iareamin(:) = narea   ;   iareamax(:) = narea   ;   iareasum(:) = 1         ! this is local information
+         ENDIF
+         !
+         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) )
+         IF( Agrif_Root() ) THEN
+            WRITE(ctmp6,*) '      ===> output of last computed fields in output.abort* files'
          ELSE
-            ! find local min and max locations
-            ih(:)  = MAXLOC( ABS( ssh(:,:,Kmm)   )                              ) + (/ nimpp - 1, njmpp - 1    /)
-            iu(:)  = MAXLOC( ABS( uu  (:,:,:,Kmm) )                              ) + (/ nimpp - 1, njmpp - 1, 0 /)
-            is1(:) = MINLOC( ts(:,:,:,jp_sal,Kmm), mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /)
-            is2(:) = MAXLOC( ts(:,:,:,jp_sal,Kmm), mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /)
-         ENDIF
-         
-         WRITE(ctmp1,*) ' stp_ctl: |ssh| > 20 m  or  |U| > 10 m/s  or  S <= 0  or  S >= 100  or  NaN encounter in the tests'
-         WRITE(ctmp2,9100) kt,   zmax(1), ih(1) , ih(2)
-         WRITE(ctmp3,9200) kt,   zmax(2), iu(1) , iu(2) , iu(3)
-         WRITE(ctmp4,9300) kt, - zmax(3), is1(1), is1(2), is1(3)
-         WRITE(ctmp5,9400) kt,   zmax(4), is2(1), is2(2), is2(3)
-         WRITE(ctmp6,*) '      ===> output of last computed fields in output.abort.nc file'
-         
+            WRITE(ctmp6,*) '      ===> output of last computed fields in '//TRIM(Agrif_CFixed())//'_output.abort* files'
+         ENDIF
+         !
          CALL dia_wri_state( Kmm, 'output.abort' )     ! create an output.abort file
-         
-         IF( .NOT. sn_cfctl%l_glochk ) THEN
-            WRITE(ctmp8,*) 'E R R O R message from sub-domain: ', narea
-            CALL ctl_stop( 'STOP', ctmp1, ' ', ctmp8, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ctmp6 )
-         ELSE
-            CALL ctl_stop( ctmp1, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ' ', ctmp6, ' ' )
-         ENDIF
-
-         kindic = -3
-         !
-      ENDIF
-      !
-9100  FORMAT (' kt=',i8,'   |ssh| max: ',1pg11.4,', at  i j  : ',2i5)
-9200  FORMAT (' kt=',i8,'   |U|   max: ',1pg11.4,', at  i j k: ',3i5)
-9300  FORMAT (' kt=',i8,'   S     min: ',1pg11.4,', at  i j k: ',3i5)
-9400  FORMAT (' kt=',i8,'   S     max: ',1pg11.4,', at  i j k: ',3i5)
+         !
+         IF( ll_colruns .or. jpnij == 1 ) THEN   ! all processes synchronized -> use lwp to print in opened ocean.output files
+            IF(lwp) THEN   ;   CALL ctl_stop( ctmp1, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ' ', ctmp6 )
+            ELSE           ;   nstop = MAX(1, nstop)   ! make sure nstop > 0 (automatically done when calling ctl_stop)
+            ENDIF
+         ELSE                                    ! only mpi subdomains with errors are here -> STOP now
+            CALL ctl_stop( 'STOP', ctmp1, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ' ', ctmp6 )
+         ENDIF
+         !
+      ENDIF
+      !
+      IF( nstop > 0 ) THEN                                                  ! an error was detected and we did not abort yet...
+         ngrdstop = Agrif_Fixed()                                           ! store which grid got this error
+         IF( .NOT. ll_colruns .AND. jpnij > 1 )   CALL ctl_stop( 'STOP' )   ! we must abort here to avoid MPI deadlock
+      ENDIF
+      !
 9500  FORMAT(' it :', i8, '    |ssh|_max: ', D23.16, ' |U|_max: ', D23.16,' S_min: ', D23.16,' S_max: ', D23.16)
       !
    END SUBROUTINE stp_ctl
+
+
+   SUBROUTINE wrt_line( cdline, kt, cdprefix, pval, kloc, ksum, kmin, kmax )
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE wrt_line  ***
+      !!
+      !! ** Purpose :   write information line
+      !!
+      !!----------------------------------------------------------------------
+      CHARACTER(len=*),      INTENT(  out) ::   cdline
+      CHARACTER(len=*),      INTENT(in   ) ::   cdprefix
+      REAL(wp),              INTENT(in   ) ::   pval
+      INTEGER, DIMENSION(3), INTENT(in   ) ::   kloc
+      INTEGER,               INTENT(in   ) ::   kt, ksum, kmin, kmax
+      !
+      CHARACTER(len=80) ::   clsuff
+      CHARACTER(len=9 ) ::   clkt, clsum, clmin, clmax
+      CHARACTER(len=9 ) ::   cli, clj, clk
+      CHARACTER(len=1 ) ::   clfmt
+      CHARACTER(len=4 ) ::   cl4   ! needed to be able to compile with Agrif, I don't know why
+      INTEGER           ::   ifmtk
+      !!----------------------------------------------------------------------
+      WRITE(clkt , '(i9)') kt
+      
+      WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpnij  ,wp))) + 1     ! how many digits to we need to write ? (we decide max = 9)
+      !!! WRITE(clsum, '(i'//clfmt//')') ksum                   ! this is creating a compilation error with AGRIF
+      cl4 = '(i'//clfmt//')'   ;   WRITE(clsum, cl4) ksum
+      WRITE(clfmt, '(i1)') INT(LOG10(REAL(MAX(1,jpnij-1),wp))) + 1    ! how many digits to we need to write ? (we decide max = 9)
+      cl4 = '(i'//clfmt//')'   ;   WRITE(clmin, cl4) kmin-1
+                                   WRITE(clmax, cl4) kmax-1
+      !
+      WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpiglo,wp))) + 1      ! how many digits to we need to write jpiglo? (we decide max = 9)
+      cl4 = '(i'//clfmt//')'   ;   WRITE(cli, cl4) kloc(1)      ! this is ok with AGRIF
+      WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpjglo,wp))) + 1      ! how many digits to we need to write jpjglo? (we decide max = 9)
+      cl4 = '(i'//clfmt//')'   ;   WRITE(clj, cl4) kloc(2)      ! this is ok with AGRIF
+      !
+      IF( ksum == 1 ) THEN   ;   WRITE(clsuff,9100) TRIM(clmin)
+      ELSE                   ;   WRITE(clsuff,9200) TRIM(clsum), TRIM(clmin), TRIM(clmax)
+      ENDIF
+      IF(kloc(3) == 0) THEN
+         ifmtk = INT(LOG10(REAL(jpk,wp))) + 1                   ! how many digits to we need to write jpk? (we decide max = 9)
+         clk = REPEAT(' ', ifmtk)                               ! create the equivalent in blank string
+         WRITE(cdline,9300) TRIM(ADJUSTL(clkt)), TRIM(ADJUSTL(cdprefix)), pval, TRIM(cli), TRIM(clj), clk(1:ifmtk), TRIM(clsuff)
+      ELSE
+         WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpk,wp))) + 1      ! how many digits to we need to write jpk? (we decide max = 9)
+         !!! WRITE(clk, '(i'//clfmt//')') kloc(3)               ! this is creating a compilation error with AGRIF
+         cl4 = '(i'//clfmt//')'   ;   WRITE(clk, cl4) kloc(3)   ! this is ok with AGRIF
+         WRITE(cdline,9400) TRIM(ADJUSTL(clkt)), TRIM(ADJUSTL(cdprefix)), pval, TRIM(cli), TRIM(clj),    TRIM(clk), TRIM(clsuff)
+      ENDIF
+      !
+9100  FORMAT('MPI rank ', a)
+9200  FORMAT('found in ', a, ' MPI tasks, spread out among ranks ', a, ' to ', a)
+9300  FORMAT('kt ', a, ' ', a, ' ', 1pg11.4, ' at i j   ', a, ' ', a, ' ', a, ' ', a)
+9400  FORMAT('kt ', a, ' ', a, ' ', 1pg11.4, ' at i j k ', a, ' ', a, ' ', a, ' ', a)
+      !
+   END SUBROUTINE wrt_line
+
 
    !!======================================================================
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OFF/nemogcm.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OFF/nemogcm.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OFF/nemogcm.F90	(revision 13193)
@@ -31,4 +31,6 @@
    USE domqco         ! tools for scale factor         (dom_qco_r3c  routine)
 #endif
+   USE bdy_oce,  ONLY : ln_bdy
+   USE bdyini         ! open boundary cond. setting       (bdy_init routine)
    !              ! ocean physics
    USE ldftra         ! lateral diffusivity setting    (ldf_tra_init routine)
@@ -93,5 +95,5 @@
       !!              Madec, 2008, internal report, IPSL.
       !!----------------------------------------------------------------------
-      INTEGER :: istp, indic       ! time step index
+      INTEGER :: istp       ! time step index
       !!----------------------------------------------------------------------
 
@@ -145,5 +147,5 @@
 # endif
 #endif         
-                                CALL stp_ctl    ( istp, indic )  ! Time loop: control and print
+                                CALL stp_ctl    ( istp )             ! Time loop: control and print
          istp = istp + 1
       END DO
@@ -160,5 +162,6 @@
       IF( nstop /= 0 .AND. lwp ) THEN                 ! error print
          WRITE(ctmp1,*) '   ==>>>   nemo_gcm: a total of ', nstop, ' errors have been found'
-         CALL ctl_stop( ctmp1 )
+         WRITE(ctmp2,*) '           Look for "E R R O R" messages in all existing ocean_output* files'
+         CALL ctl_stop( ' ', ctmp1, ' ', ctmp2 )
       ENDIF
       !
@@ -242,20 +245,6 @@
       !
       ! finalize the definition of namctl variables
-      IF( sn_cfctl%l_allon ) THEN
-         ! Turn on all options.
-         CALL nemo_set_cfctl( sn_cfctl, .TRUE., .TRUE. )
-         ! Ensure all processors are active
-         sn_cfctl%procmin = 0 ; sn_cfctl%procmax = 1000000 ; sn_cfctl%procincr = 1
-      ELSEIF( sn_cfctl%l_config ) THEN
-         ! Activate finer control of report outputs
-         ! optionally switch off output from selected areas (note this only
-         ! applies to output which does not involve global communications)
-         IF( ( narea < sn_cfctl%procmin .OR. narea > sn_cfctl%procmax  ) .OR. &
-           & ( MOD( narea - sn_cfctl%procmin, sn_cfctl%procincr ) /= 0 ) )    &
-           &   CALL nemo_set_cfctl( sn_cfctl, .FALSE., .FALSE. )
-      ELSE
-         ! turn off all options.
-         CALL nemo_set_cfctl( sn_cfctl, .FALSE., .TRUE. )
-      ENDIF
+      IF( narea < sn_cfctl%procmin .OR. narea > sn_cfctl%procmax .OR. MOD( narea - sn_cfctl%procmin, sn_cfctl%procincr ) /= 0 )   &
+         &   CALL nemo_set_cfctl( sn_cfctl, .FALSE. )
       !
       lwp = (narea == 1) .OR. sn_cfctl%l_oceout    ! control of all listing output print
@@ -322,5 +311,4 @@
       ! Initialise time level indices
       Nbb = 1; Nnn = 2; Naa = 3; Nrhs = Naa
-	
 
       !                             !-------------------------------!
@@ -344,4 +332,5 @@
 
                            CALL     sbc_init( Nbb, Nnn, Naa )    ! Forcings : surface module
+                           CALL     bdy_init    ! Open boundaries initialisation    
 
       !                                      ! Tracer physics
@@ -386,7 +375,4 @@
          WRITE(numout,*) '~~~~~~~~'
          WRITE(numout,*) '   Namelist namctl'
-         WRITE(numout,*) '                              sn_cfctl%l_glochk  = ', sn_cfctl%l_glochk
-         WRITE(numout,*) '                              sn_cfctl%l_allon   = ', sn_cfctl%l_allon
-         WRITE(numout,*) '       finer control over o/p sn_cfctl%l_config  = ', sn_cfctl%l_config
          WRITE(numout,*) '                              sn_cfctl%l_runstat = ', sn_cfctl%l_runstat
          WRITE(numout,*) '                              sn_cfctl%l_trcstat = ', sn_cfctl%l_trcstat
@@ -507,4 +493,5 @@
       USE zdf_oce,   ONLY : zdf_oce_alloc
       USE trc_oce,   ONLY : trc_oce_alloc
+      USE bdy_oce,   ONLY : bdy_oce_alloc
       !
       INTEGER :: ierr
@@ -516,4 +503,5 @@
       ierr = ierr + zdf_oce_alloc()          ! ocean vertical physics
       ierr = ierr + trc_oce_alloc()          ! shared TRC / TRA arrays
+      ierr = ierr + bdy_oce_alloc()    ! bdy masks (incl. initialization)      
       !
       CALL mpp_sum( 'nemogcm', ierr )
@@ -522,25 +510,18 @@
    END SUBROUTINE nemo_alloc
 
-   SUBROUTINE nemo_set_cfctl(sn_cfctl, setto, for_all )
+   SUBROUTINE nemo_set_cfctl(sn_cfctl, setto )
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE nemo_set_cfctl  ***
       !!
       !! ** Purpose :   Set elements of the output control structure to setto.
-      !!                for_all should be .false. unless all areas are to be
-      !!                treated identically.
-      !!
+     !!
       !! ** Method  :   Note this routine can be used to switch on/off some
-      !!                types of output for selected areas but any output types
-      !!                that involve global communications (e.g. mpp_max, glob_sum)
-      !!                should be protected from selective switching by the
-      !!                for_all argument
-      !!----------------------------------------------------------------------
-      LOGICAL :: setto, for_all
-      TYPE(sn_ctl) :: sn_cfctl
-      !!----------------------------------------------------------------------
-      IF( for_all ) THEN
-         sn_cfctl%l_runstat = setto
-         sn_cfctl%l_trcstat = setto
-      ENDIF
+      !!                types of output for selected areas.
+      !!----------------------------------------------------------------------
+      TYPE(sn_ctl), INTENT(inout) :: sn_cfctl
+      LOGICAL     , INTENT(in   ) :: setto
+      !!----------------------------------------------------------------------
+      sn_cfctl%l_runstat = setto
+      sn_cfctl%l_trcstat = setto
       sn_cfctl%l_oceout  = setto
       sn_cfctl%l_layout  = setto
@@ -572,5 +553,5 @@
 
 
-   SUBROUTINE stp_ctl( kt, kindic )
+   SUBROUTINE stp_ctl( kt )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE stp_ctl  ***
@@ -583,5 +564,4 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in   ) ::   kt      ! ocean time-step index
-      INTEGER, INTENT(inout) ::   kindic  ! indicator of solver convergence
       !!----------------------------------------------------------------------
       !
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/SAO/nemogcm.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/SAO/nemogcm.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/SAO/nemogcm.F90	(revision 13193)
@@ -158,20 +158,6 @@
       !
       ! finalize the definition of namctl variables
-      IF( sn_cfctl%l_allon ) THEN
-         ! Turn on all options.
-         CALL nemo_set_cfctl( sn_cfctl, .TRUE., .TRUE. )
-         ! Ensure all processors are active
-         sn_cfctl%procmin = 0 ; sn_cfctl%procmax = 1000000 ; sn_cfctl%procincr = 1
-      ELSEIF( sn_cfctl%l_config ) THEN
-         ! Activate finer control of report outputs
-         ! optionally switch off output from selected areas (note this only
-         ! applies to output which does not involve global communications)
-         IF( ( narea < sn_cfctl%procmin .OR. narea > sn_cfctl%procmax  ) .OR. &
-           & ( MOD( narea - sn_cfctl%procmin, sn_cfctl%procincr ) /= 0 ) )    &
-           &   CALL nemo_set_cfctl( sn_cfctl, .FALSE., .FALSE. )
-      ELSE
-         ! turn off all options.
-         CALL nemo_set_cfctl( sn_cfctl, .FALSE., .TRUE. )
-      ENDIF
+      IF( narea < sn_cfctl%procmin .OR. narea > sn_cfctl%procmax .OR. MOD( narea - sn_cfctl%procmin, sn_cfctl%procincr ) /= 0 )   &
+         &   CALL nemo_set_cfctl( sn_cfctl, .FALSE. )
       !
       lwp = (narea == 1) .OR. sn_cfctl%l_oceout    ! control of all listing output print
@@ -270,7 +256,4 @@
          WRITE(numout,*) '~~~~~~~~'
          WRITE(numout,*) '   Namelist namctl'
-         WRITE(numout,*) '                              sn_cfctl%l_glochk  = ', sn_cfctl%l_glochk
-         WRITE(numout,*) '                              sn_cfctl%l_allon   = ', sn_cfctl%l_allon
-         WRITE(numout,*) '       finer control over o/p sn_cfctl%l_config  = ', sn_cfctl%l_config
          WRITE(numout,*) '                              sn_cfctl%l_runstat = ', sn_cfctl%l_runstat
          WRITE(numout,*) '                              sn_cfctl%l_trcstat = ', sn_cfctl%l_trcstat
@@ -410,25 +393,18 @@
    END SUBROUTINE nemo_alloc
 
-   SUBROUTINE nemo_set_cfctl(sn_cfctl, setto, for_all )
+   SUBROUTINE nemo_set_cfctl(sn_cfctl, setto )
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE nemo_set_cfctl  ***
       !!
       !! ** Purpose :   Set elements of the output control structure to setto.
-      !!                for_all should be .false. unless all areas are to be
-      !!                treated identically.
       !!
       !! ** Method  :   Note this routine can be used to switch on/off some
-      !!                types of output for selected areas but any output types
-      !!                that involve global communications (e.g. mpp_max, glob_sum)
-      !!                should be protected from selective switching by the
-      !!                for_all argument
-      !!----------------------------------------------------------------------
-      LOGICAL :: setto, for_all
-      TYPE(sn_ctl) :: sn_cfctl
-      !!----------------------------------------------------------------------
-      IF( for_all ) THEN
-         sn_cfctl%l_runstat = setto
-         sn_cfctl%l_trcstat = setto
-      ENDIF
+      !!                types of output for selected areas.
+      !!----------------------------------------------------------------------
+      TYPE(sn_ctl), INTENT(inout) :: sn_cfctl
+      LOGICAL     , INTENT(in   ) :: setto
+      !!----------------------------------------------------------------------
+      sn_cfctl%l_runstat = setto
+      sn_cfctl%l_trcstat = setto
       sn_cfctl%l_oceout  = setto
       sn_cfctl%l_layout  = setto
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/SAS/diawri.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/SAS/diawri.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/SAS/diawri.F90	(revision 13193)
@@ -138,7 +138,6 @@
       !!      Each nn_write time step, output the instantaneous or mean fields
       !!----------------------------------------------------------------------
-      !!
       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
-      INTEGER, INTENT( in ) ::   Kmm  ! ocean time level index
+      INTEGER, INTENT( in ) ::   Kmm     ! ocean time level index
       !!
       LOGICAL ::   ll_print = .FALSE.                        ! =T print and flush numout
@@ -462,6 +461,6 @@
          CALL iom_close( inum )
       ENDIF
-#endif
-
+      !
+#endif
    END SUBROUTINE dia_wri_state
 
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/SAS/nemogcm.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/SAS/nemogcm.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/SAS/nemogcm.F90	(revision 13193)
@@ -126,10 +126,4 @@
       END DO
       !
-      IF( .NOT. Agrif_Root() ) THEN
-         CALL Agrif_ParentGrid_To_ChildGrid()
-         IF( ln_timing )   CALL timing_finalize
-         CALL Agrif_ChildGrid_To_ParentGrid()
-      ENDIF
-      !
 #else
       !
@@ -166,5 +160,13 @@
       IF( nstop /= 0 .AND. lwp ) THEN        ! error print
          WRITE(ctmp1,*) '   ==>>>   nemo_gcm: a total of ', nstop, ' errors have been found'
-         CALL ctl_stop( ctmp1 )
+         IF( ngrdstop > 0 ) THEN
+            WRITE(ctmp9,'(i2)') ngrdstop
+            WRITE(ctmp2,*) '           E R R O R detected in Agrif grid '//TRIM(ctmp9)
+            WRITE(ctmp3,*) '           Look for "E R R O R" messages in all existing '//TRIM(ctmp9)//'_ocean_output* files'
+            CALL ctl_stop( ' ', ctmp1, ' ', ctmp2, ' ', ctmp3 )
+         ELSE
+            WRITE(ctmp2,*) '           Look for "E R R O R" messages in all existing ocean_output* files'
+            CALL ctl_stop( ' ', ctmp1, ' ', ctmp2 )
+         ENDIF
       ENDIF
       !
@@ -275,20 +277,6 @@
       !
       ! finalize the definition of namctl variables
-      IF( sn_cfctl%l_allon ) THEN
-         ! Turn on all options.
-         CALL nemo_set_cfctl( sn_cfctl, .TRUE., .TRUE. )
-         ! Ensure all processors are active
-         sn_cfctl%procmin = 0 ; sn_cfctl%procmax = 1000000 ; sn_cfctl%procincr = 1
-      ELSEIF( sn_cfctl%l_config ) THEN
-         ! Activate finer control of report outputs
-         ! optionally switch off output from selected areas (note this only
-         ! applies to output which does not involve global communications)
-         IF( ( narea < sn_cfctl%procmin .OR. narea > sn_cfctl%procmax  ) .OR. &
-           & ( MOD( narea - sn_cfctl%procmin, sn_cfctl%procincr ) /= 0 ) )    &
-           &   CALL nemo_set_cfctl( sn_cfctl, .FALSE., .FALSE. )
-      ELSE
-         ! turn off all options.
-         CALL nemo_set_cfctl( sn_cfctl, .FALSE., .TRUE. )
-      ENDIF
+      IF( narea < sn_cfctl%procmin .OR. narea > sn_cfctl%procmax .OR. MOD( narea - sn_cfctl%procmin, sn_cfctl%procincr ) /= 0 )   &
+         &   CALL nemo_set_cfctl( sn_cfctl, .FALSE. )
       !
       lwp = (narea == 1) .OR. sn_cfctl%l_oceout    ! control of all listing output print
@@ -408,7 +396,4 @@
          WRITE(numout,*) '~~~~~~~~'
          WRITE(numout,*) '   Namelist namctl'
-         WRITE(numout,*) '                              sn_cfctl%l_glochk  = ', sn_cfctl%l_glochk
-         WRITE(numout,*) '                              sn_cfctl%l_allon   = ', sn_cfctl%l_allon
-         WRITE(numout,*) '       finer control over o/p sn_cfctl%l_config  = ', sn_cfctl%l_config
          WRITE(numout,*) '                              sn_cfctl%l_runstat = ', sn_cfctl%l_runstat
          WRITE(numout,*) '                              sn_cfctl%l_trcstat = ', sn_cfctl%l_trcstat
@@ -552,25 +537,18 @@
    END SUBROUTINE nemo_alloc
 
-   SUBROUTINE nemo_set_cfctl(sn_cfctl, setto, for_all )
+   SUBROUTINE nemo_set_cfctl(sn_cfctl, setto )
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE nemo_set_cfctl  ***
       !!
       !! ** Purpose :   Set elements of the output control structure to setto.
-      !!                for_all should be .false. unless all areas are to be
-      !!                treated identically.
       !!
       !! ** Method  :   Note this routine can be used to switch on/off some
-      !!                types of output for selected areas but any output types
-      !!                that involve global communications (e.g. mpp_max, glob_sum)
-      !!                should be protected from selective switching by the
-      !!                for_all argument
-      !!----------------------------------------------------------------------
-      LOGICAL :: setto, for_all
-      TYPE(sn_ctl) :: sn_cfctl
-      !!----------------------------------------------------------------------
-      IF( for_all ) THEN
-         sn_cfctl%l_runstat = setto
-         sn_cfctl%l_trcstat = setto
-      ENDIF
+      !!                types of output for selected areas.
+      !!----------------------------------------------------------------------
+      TYPE(sn_ctl), INTENT(inout) :: sn_cfctl
+      LOGICAL     , INTENT(in   ) :: setto
+      !!----------------------------------------------------------------------
+      sn_cfctl%l_runstat = setto
+      sn_cfctl%l_trcstat = setto
       sn_cfctl%l_oceout  = setto
       sn_cfctl%l_layout  = setto
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/SAS/step.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/SAS/step.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/SAS/step.F90	(revision 13193)
@@ -74,24 +74,18 @@
       !!              -2- Outputs and diagnostics
       !!----------------------------------------------------------------------
-      INTEGER ::   indic    ! error indicator if < 0
-      !! ---------------------------------------------------------------------
 
 #if defined key_agrif
-      IF( nstop > 0 ) return   ! avoid to go further if an error was detected during previous time step 
+      IF( nstop > 0 ) RETURN   ! avoid to go further if an error was detected during previous time step (child grid)
       kstp = nit000 + Agrif_Nb_Step()
       Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs   ! agrif_oce module copies of time level indices
-      IF ( lk_agrif_debug ) THEN
-         IF ( Agrif_Root() .and. lwp) Write(*,*) '---'
-         IF (lwp) Write(*,*) 'Grid Number',Agrif_Fixed(),' time step ',kstp, 'int tstep',Agrif_NbStepint()
+      IF( lk_agrif_debug ) THEN
+         IF( Agrif_Root() .and. lwp)   WRITE(*,*) '---'
+         IF(lwp)   WRITE(*,*) 'Grid Number', Agrif_Fixed(),' time step ', kstp, 'int tstep', Agrif_NbStepint()
       ENDIF
-
-      IF ( kstp == (nit000 + 1) ) lk_agrif_fstep = .FALSE.
-
+      IF( kstp == nit000 + 1 )   lk_agrif_fstep = .FALSE.
 # if defined key_iomput
       IF( Agrif_Nbstepint() == 0 )   CALL iom_swap( cxios_context )
 # endif   
 #endif   
-                             indic = 0                    ! although indic is not changed in stp_ctl
-                                                          ! need to keep the same interface 
       IF( kstp == nit000 )   CALL iom_init( cxios_context ) ! iom_put initialization (must be done after nemo_init for AGRIF+XIOS+OASIS)
       IF( kstp /= nit000 )   CALL day( kstp )             ! Calendar (day was already called at nit000 in day_init)
@@ -112,15 +106,12 @@
       ! AGRIF recursive integration
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<      
-                             CALL Agrif_Integrate_ChildGrids( stp )  
-#endif
+                             CALL Agrif_Integrate_ChildGrids( stp )
                              
+#endif                             
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Control
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                             CALL stp_ctl( kstp, indic )
-      IF( indic < 0  )  THEN
-                             CALL ctl_stop( 'step: indic < 0' )
-                             CALL dia_wri_state( Nnn, 'output.abort' )
-      ENDIF
+                             CALL stp_ctl( kstp, Nnn )
+
 #if defined key_agrif
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -132,4 +123,5 @@
 #endif
       ENDIF
+
 #endif
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -141,5 +133,5 @@
       ! Coupled mode
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( lk_oasis    )  CALL sbc_cpl_snd( kstp, Nbb, Nnn )     ! coupled mode : field exchanges if OASIS-coupled ice
+      IF( lk_oasis .AND. nstop == 0 ) CALL sbc_cpl_snd( kstp, Nbb, Nnn )       ! coupled mode : field exchanges if OASIS-coupled ice
 
 #if defined key_iomput
@@ -152,6 +144,6 @@
          lrst_oce = .FALSE.
       ENDIF
-      IF( kstp == nitend .OR. indic < 0 ) THEN
-                             CALL iom_context_finalize( cxios_context ) ! needed for XIOS+AGRIF
+      IF( kstp == nitend .OR. nstop > 0 ) THEN
+         CALL iom_context_finalize( cxios_context ) ! needed for XIOS+AGRIF
       ENDIF
 #endif
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/SAS/stpctl.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/SAS/stpctl.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/SAS/stpctl.F90	(revision 13193)
@@ -21,8 +21,9 @@
    USE ice      , ONLY : vt_i, u_ice, tm_i
    !
+   USE diawri          ! Standard run outputs       (dia_wri_state routine)
    USE in_out_manager  ! I/O manager
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE lib_mpp         ! distributed memory computing
-
+   !
    USE netcdf          ! NetCDF library
    IMPLICIT NONE
@@ -31,6 +32,6 @@
    PUBLIC stp_ctl           ! routine called by step.F90
 
-   INTEGER  ::   idrun, idtime, idssh, idu, ids, istatus
-   LOGICAL  ::   lsomeoce
+   INTEGER                ::   nrunid   ! netcdf file id
+   INTEGER, DIMENSION(3)  ::   nvarid   ! netcdf variable id
    !!----------------------------------------------------------------------
    !! NEMO/SAS 4.0 , NEMO Consortium (2018)
@@ -38,8 +39,7 @@
    !! Software governed by the CeCILL license (see ./LICENSE)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
-   SUBROUTINE stp_ctl( kt, kindic )
+   SUBROUTINE stp_ctl( kt, Kmm )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE stp_ctl  ***
@@ -49,66 +49,217 @@
       !! ** Method  : - Save the time step in numstp
       !!              - Print it each 50 time steps
+      !!              - Stop the run IF problem encountered by setting nstop > 0
+      !!                Problems checked: ice thickness maximum > 100 m
+      !!                                  ice velocity  maximum > 10 m/s 
+      !!                                  min ice temperature   < -100 degC
       !!
       !! ** Actions :   "time.step" file = last ocean time-step
       !!                "run.stat"  file = run statistics
-      !!                
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT( in    ) ::   kt       ! ocean time-step index
-      INTEGER, INTENT( inout ) ::   kindic   ! indicator of solver convergence
-      !!
-      REAL(wp), DIMENSION(3) ::   zmax
-      LOGICAL                ::   ll_wrtstp, ll_colruns, ll_wrtruns
-      CHARACTER(len=20) :: clname
-      !!----------------------------------------------------------------------
-      !
-      ll_wrtstp  = ( MOD( kt, sn_cfctl%ptimincr ) == 0 ) .OR. ( kt == nitend )
-      ll_colruns = ll_wrtstp .AND. ( sn_cfctl%l_runstat )
-      ll_wrtruns = ll_colruns .AND. lwm
-      IF( kt == nit000 .AND. lwp ) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'stp_ctl : time-stepping control'
-         WRITE(numout,*) '~~~~~~~'
-         !                                ! open time.step file
-         IF( lwm ) CALL ctl_opn( numstp, 'time.step', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
-         !                                ! open run.stat file(s) at start whatever
-         !                                ! the value of sn_cfctl%ptimincr
-         IF( lwm .AND. ( sn_cfctl%l_runstat ) ) THEN
+      !!                 nstop indicator sheared among all local domain
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in   ) ::   kt       ! ocean time-step index
+      INTEGER, INTENT(in   ) ::   Kmm      ! ocean time level index
+      !!
+      INTEGER                         ::   ji                                    ! dummy loop indices
+      INTEGER                         ::   idtime, istatus
+      INTEGER , DIMENSION(4)          ::   iareasum, iareamin, iareamax
+      INTEGER , DIMENSION(3,3)        ::   iloc                                  ! min/max loc indices
+      REAL(wp)                        ::   zzz                                   ! local real 
+      REAL(wp), DIMENSION(4)          ::   zmax, zmaxlocal
+      LOGICAL                         ::   ll_wrtstp, ll_colruns, ll_wrtruns
+      LOGICAL, DIMENSION(jpi,jpj)     ::   llmsk
+      CHARACTER(len=20)               ::   clname
+      !!----------------------------------------------------------------------
+      IF( nstop > 0 .AND. ngrdstop > -1 )   RETURN   !   stpctl was already called by a child grid
+      !
+      ll_wrtstp  = ( MOD( kt-nit000, sn_cfctl%ptimincr ) == 0 ) .OR. ( kt == nitend )
+      ll_colruns = ll_wrtstp .AND. sn_cfctl%l_runstat .AND. jpnij > 1 
+      ll_wrtruns = ( ll_colruns .OR. jpnij == 1 ) .AND. lwm
+      !
+      IF( kt == nit000 ) THEN
+         !
+         IF( lwp ) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) 'stp_ctl : time-stepping control'
+            WRITE(numout,*) '~~~~~~~'
+         ENDIF
+         !                                ! open time.step    ascii file, done only by 1st subdomain
+         IF( lwm )   CALL ctl_opn( numstp, 'time.step', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
+         !
+         IF( ll_wrtruns ) THEN
+            !                             ! open run.stat     ascii file, done only by 1st subdomain
             CALL ctl_opn( numrun, 'run.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
+            !                             ! open run.stat.nc netcdf file, done only by 1st subdomain
             clname = 'run.stat.nc'
             IF( .NOT. Agrif_Root() )   clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname)
-            istatus = NF90_CREATE( 'run.stat.nc', NF90_CLOBBER, idrun )
-            istatus = NF90_DEF_DIM( idrun, 'time'     , NF90_UNLIMITED, idtime )
-            istatus = NF90_DEF_VAR( idrun, 'vt_i_max' , NF90_DOUBLE, (/ idtime /), idssh )
-            istatus = NF90_DEF_VAR( idrun, 'abs_u_max', NF90_DOUBLE, (/ idtime /), idu )
-            istatus = NF90_DEF_VAR( idrun, 'tm_i_min' , NF90_DOUBLE, (/ idtime /), ids )
-            istatus = NF90_ENDDEF(idrun)
-         ENDIF
-      ENDIF
-      IF( kt == nit000 )   lsomeoce = COUNT( ssmask(:,:) == 1._wp ) > 0
-      !
-      IF(lwm .AND. ll_wrtstp) THEN        !==  current time step  ==!   ("time.step" file)
+            istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, nrunid )
+            istatus = NF90_DEF_DIM( nrunid, 'time'     , NF90_UNLIMITED, idtime )
+            istatus = NF90_DEF_VAR( nrunid, 'vt_i_max' , NF90_DOUBLE, (/ idtime /), nvarid(1) )
+            istatus = NF90_DEF_VAR( nrunid, 'abs_u_max', NF90_DOUBLE, (/ idtime /), nvarid(2) )
+            istatus = NF90_DEF_VAR( nrunid, 'tm_i_min' , NF90_DOUBLE, (/ idtime /), nvarid(3) )
+            istatus = NF90_ENDDEF(nrunid)
+         ENDIF
+         !    
+      ENDIF
+      !
+      !                                   !==              write current time step              ==!
+      !                                   !==  done only by 1st subdomain at writting timestep  ==!
+      IF( lwm .AND. ll_wrtstp ) THEN
          WRITE ( numstp, '(1x, i8)' )   kt
          REWIND( numstp )
       ENDIF
-      !                                   !==  test of extrema  ==!
+      !                                   !==            test of local extrema           ==!
+      !                                   !==  done by all processes at every time step  ==!
+      llmsk(:,:) = tmask(:,:,1) == 1._wp
+      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 (:,:) + 273.15_wp, mask = llmsk )   ! min ice temperature
+      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
+      !                                   !==               get global extrema             ==!
+      !                                   !==  done by all processes if writting run.stat  ==!
       IF( ll_colruns ) THEN
-         zmax(1) = MAXVAL(      vt_i (:,:) )                                           ! max ice thickness
-         zmax(2) = MAXVAL( ABS( u_ice(:,:) ) )                                         ! max ice velocity (zonal only)
-         zmax(3) = MAXVAL(     -tm_i (:,:)+273.15_wp , mask = ssmask(:,:) == 1._wp )   ! min ice temperature
-         CALL mpp_max( "stpctl", zmax )                                   ! max over the global domain
+         zmaxlocal(:) = zmax(:)
+         CALL mpp_max( "stpctl", zmax )          ! max over the global domain
+         nstop = NINT( zmax(4) )                 ! update nstop indicator (now sheared among all local domains)
+      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)
+         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/) )
+         IF( kt == nitend )   istatus = NF90_CLOSE(nrunid)
       END IF
-      !                                            !==  run statistics  ==!   ("run.stat" file)
-      IF( ll_wrtruns ) THEN
-         WRITE(numrun,9500) kt, zmax(1), zmax(2), - zmax(3)
-         istatus = NF90_PUT_VAR( idrun, idssh, (/ zmax(1)/), (/kt/), (/1/) )
-         istatus = NF90_PUT_VAR( idrun,   idu, (/ zmax(2)/), (/kt/), (/1/) )
-         istatus = NF90_PUT_VAR( idrun,   ids, (/-zmax(3)/), (/kt/), (/1/) )
-         IF( MOD( kt , 100 ) == 0 ) istatus = NF90_SYNC(idrun)
-         IF( kt == nitend         ) istatus = NF90_CLOSE(idrun)
-      END IF
+      !                                   !==               error handling               ==!
+      !                                   !==  done by all processes at every time step  ==!
+      !
+      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)
+         &  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
+         IF( ll_colruns ) THEN   ! zmax is global, so it is the same on all subdomains -> no dead lock with mpp_maxloc
+            ! first: close the netcdf file, so we can read it
+            IF( lwm .AND. kt /= nitend )   istatus = NF90_CLOSE(nrunid)
+            ! get global loc on the min/max
+            CALL mpp_maxloc( 'stpctl',      vt_i(:,:)            , tmask(:,:,1), zzz, iloc(1:2,1) )   ! mpp_maxloc ok if mask = F 
+            CALL mpp_maxloc( 'stpctl',ABS( u_ice(:,:) )          , tmask(:,:,1), zzz, iloc(1:2,2) )
+            CALL mpp_minloc( 'stpctl',      tm_i(:,:) - 273.15_wp, tmask(:,:,1), zzz, iloc(1:2,3) )
+            ! find which subdomain has the max.
+            iareamin(:) = jpnij+1   ;   iareamax(:) = 0   ;   iareasum(:) = 0
+            DO ji = 1, 4
+               IF( zmaxlocal(ji) == zmax(ji) ) THEN
+                  iareamin(ji) = narea   ;   iareamax(ji) = narea   ;   iareasum(ji) = 1
+               ENDIF
+            END DO
+            CALL mpp_min( "stpctl", iareamin )         ! min over the global domain
+            CALL mpp_max( "stpctl", iareamax )         ! max over the global domain
+            CALL mpp_sum( "stpctl", iareasum )         ! sum over the global domain
+         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
+            iloc(1:2,1) = MAXLOC(       vt_i(:,:)            , mask = llmsk ) + (/ nimpp - 1, njmpp - 1/)
+            iloc(1:2,2) = MAXLOC( ABS( u_ice(:,:) )          , mask = llmsk ) + (/ nimpp - 1, njmpp - 1/)
+            iloc(1:2,3) = MINLOC(       tm_i(:,:) - 273.15_wp, mask = llmsk ) + (/ nimpp - 1, njmpp - 1/)
+            iareamin(:) = narea   ;   iareamax(:) = narea   ;   iareasum(:) = 1         ! this is local information
+         ENDIF
+         !
+         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) )
+         IF( Agrif_Root() ) THEN
+            WRITE(ctmp6,*) '      ===> output of last computed fields in output.abort* files'
+         ELSE
+            WRITE(ctmp6,*) '      ===> output of last computed fields in '//TRIM(Agrif_CFixed())//'_output.abort* files'
+         ENDIF
+         !
+         CALL dia_wri_state( Kmm, 'output.abort' )     ! create an output.abort file
+         !
+         IF( ll_colruns .or. jpnij == 1 ) THEN   ! all processes synchronized -> use lwp to print in opened ocean.output files
+            IF(lwp) THEN   ;   CALL ctl_stop( ctmp1, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ' ', ctmp6 )
+            ELSE           ;   nstop = MAX(1, nstop)   ! make sure nstop > 0 (automatically done when calling ctl_stop)
+            ENDIF
+         ELSE                                    ! only mpi subdomains with errors are here -> STOP now
+            CALL ctl_stop( 'STOP', ctmp1, ' ', ctmp2, ctmp3, ctmp4, ctmp5, ' ', ctmp6 )
+         ENDIF
+         !
+      ENDIF
+      !
+      IF( nstop > 0 ) THEN                                                  ! an error was detected and we did not abort yet...
+         ngrdstop = Agrif_Fixed()                                           ! store which grid got this error
+         IF( .NOT. ll_colruns .AND. jpnij > 1 )   CALL ctl_stop( 'STOP' )   ! we must abort here to avoid MPI deadlock
+      ENDIF
       !
 9500  FORMAT(' it :', i8, '    vt_i_max: ', D23.16, ' |u|_max: ', D23.16,' tm_i_min: ', D23.16)
       !
    END SUBROUTINE stp_ctl
+
+
+   SUBROUTINE wrt_line( cdline, kt, cdprefix, pval, kloc, ksum, kmin, kmax )
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE wrt_line  ***
+      !!
+      !! ** Purpose :   write information line
+      !!
+      !!----------------------------------------------------------------------
+      CHARACTER(len=*),      INTENT(  out) ::   cdline
+      CHARACTER(len=*),      INTENT(in   ) ::   cdprefix
+      REAL(wp),              INTENT(in   ) ::   pval
+      INTEGER, DIMENSION(3), INTENT(in   ) ::   kloc
+      INTEGER,               INTENT(in   ) ::   kt, ksum, kmin, kmax
+      !
+      CHARACTER(len=80) ::   clsuff
+      CHARACTER(len=9 ) ::   clkt, clsum, clmin, clmax
+      CHARACTER(len=9 ) ::   cli, clj, clk
+      CHARACTER(len=1 ) ::   clfmt
+      CHARACTER(len=4 ) ::   cl4   ! needed to be able to compile with Agrif, I don't know why
+      INTEGER           ::   ifmtk
+      !!----------------------------------------------------------------------
+      WRITE(clkt , '(i9)') kt
+      
+      WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpnij  ,wp))) + 1     ! how many digits to we need to write ? (we decide max = 9)
+      !!! WRITE(clsum, '(i'//clfmt//')') ksum                   ! this is creating a compilation error with AGRIF
+      cl4 = '(i'//clfmt//')'   ;   WRITE(clsum, cl4) ksum
+      WRITE(clfmt, '(i1)') INT(LOG10(REAL(MAX(1,jpnij-1),wp))) + 1    ! how many digits to we need to write ? (we decide max = 9)
+      cl4 = '(i'//clfmt//')'   ;   WRITE(clmin, cl4) kmin-1
+                                   WRITE(clmax, cl4) kmax-1
+      !
+      WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpiglo,wp))) + 1      ! how many digits to we need to write jpiglo? (we decide max = 9)
+      cl4 = '(i'//clfmt//')'   ;   WRITE(cli, cl4) kloc(1)      ! this is ok with AGRIF
+      WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpjglo,wp))) + 1      ! how many digits to we need to write jpjglo? (we decide max = 9)
+      cl4 = '(i'//clfmt//')'   ;   WRITE(clj, cl4) kloc(2)      ! this is ok with AGRIF
+      !
+      IF( ksum == 1 ) THEN   ;   WRITE(clsuff,9100) TRIM(clmin)
+      ELSE                   ;   WRITE(clsuff,9200) TRIM(clsum), TRIM(clmin), TRIM(clmax)
+      ENDIF
+      IF(kloc(3) == 0) THEN
+         ifmtk = INT(LOG10(REAL(jpk,wp))) + 1                   ! how many digits to we need to write jpk? (we decide max = 9)
+         clk = REPEAT(' ', ifmtk)                               ! create the equivalent in blank string
+         WRITE(cdline,9300) TRIM(ADJUSTL(clkt)), TRIM(ADJUSTL(cdprefix)), pval, TRIM(cli), TRIM(clj), clk(1:ifmtk), TRIM(clsuff)
+      ELSE
+         WRITE(clfmt, '(i1)') INT(LOG10(REAL(jpk,wp))) + 1      ! how many digits to we need to write jpk? (we decide max = 9)
+         !!! WRITE(clk, '(i'//clfmt//')') kloc(3)               ! this is creating a compilation error with AGRIF
+         cl4 = '(i'//clfmt//')'   ;   WRITE(clk, cl4) kloc(3)   ! this is ok with AGRIF
+         WRITE(cdline,9400) TRIM(ADJUSTL(clkt)), TRIM(ADJUSTL(cdprefix)), pval, TRIM(cli), TRIM(clj),    TRIM(clk), TRIM(clsuff)
+      ENDIF
+      !
+9100  FORMAT('MPI rank ', a)
+9200  FORMAT('found in ', a, ' MPI tasks, spread out among ranks ', a, ' to ', a)
+9300  FORMAT('kt ', a, ' ', a, ' ', 1pg11.4, ' at i j   ', a, ' ', a, ' ', a, ' ', a)
+9400  FORMAT('kt ', a, ' ', a, ' ', 1pg11.4, ' at i j k ', a, ' ', a, ' ', a, ' ', a)
+      !
+   END SUBROUTINE wrt_line
+
 
    !!======================================================================
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P4Z/p4zmeso.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P4Z/p4zmeso.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P4Z/p4zmeso.F90	(revision 13193)
@@ -69,5 +69,6 @@
       REAL(wp) :: zfact   , zfood, zfoodlim, zproport, zbeta
       REAL(wp) :: zmortzgoc, zfrac, zfracfe, zratio, zratio2, zfracal, zgrazcal
-      REAL(wp) :: zepsherf, zepshert, zepsherv, zgrarsig, zgraztotc, zgraztotn, zgraztotf
+      REAL(wp) :: zepsherf, zepshert, zepsherv, zepsherq 
+      REAL(wp) :: zgrarsig, zgraztotc, zgraztotn, zgraztotf
       REAL(wp) :: zgrarem2, zgrafer2, zgrapoc2, zprcaca, zmortz, zgrasrat, zgrasratn
       REAL(wp) :: zrespz, ztortz, zgrazd, zgrazz, zgrazpof
@@ -156,6 +157,12 @@
          zgrazing2(ji,jj,jk) = zgraztotc
 
-         !    Mesozooplankton efficiency
-         !    --------------------------
+         ! Mesozooplankton efficiency. 
+         ! We adopt a formulation proposed by Mitra et al. (2007)
+         ! The gross growth efficiency is controled by the most limiting nutrient.
+         ! Growth is also further decreased when the food quality is poor. This is currently
+         ! hard coded : it can be decreased by up to 50% (zepsherq)
+         ! GGE can also be decreased when food quantity is high, zepsherf (Montagnes and 
+         ! Fulton, 2012)
+         ! -----------------------------------------------------------------------------------
          zgrasrat  =  ( zgraztotf + rtrn )/ ( zgraztotc + rtrn )
          zgrasratn =  ( zgraztotn + rtrn )/ ( zgraztotc + rtrn )
@@ -163,5 +170,6 @@
          zbeta     = MAX(0., (epsher2 - epsher2min) )
          zepsherf  = epsher2min + zbeta / ( 1.0 + 0.04E6 * 12. * zfood * zbeta ) 
-         zepsherv  = zepsherf * zepshert 
+         zepsherq  = 0.5 + (1.0 - 0.5) * zepshert * ( 1.0 + 1.0 ) / ( zepshert + 1.0 )
+         zepsherv  = zepsherf * zepshert * zepsherq 
 
          zgrarem2  = zgraztotc * ( 1. - zepsherv - unass2 ) &
@@ -170,4 +178,5 @@
          &         + ferat3 * ( ( 1. - epsher2 - unass2 ) /( 1. - epsher2 ) * ztortz )
          zgrapoc2  = zgraztotc * unass2
+
 
          !   Update the arrays TRA which contain the biological sources and sinks
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P4Z/p4zmicro.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P4Z/p4zmicro.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P4Z/p4zmicro.F90	(revision 13193)
@@ -67,5 +67,6 @@
       REAL(wp) :: zgraze  , zdenom, zdenom2
       REAL(wp) :: zfact   , zfood, zfoodlim, zbeta
-      REAL(wp) :: zepsherf, zepshert, zepsherv, zgrarsig, zgraztotc, zgraztotn, zgraztotf
+      REAL(wp) :: zepsherf, zepshert, zepsherv, zepsherq
+      REAL(wp) :: zgrarsig, zgraztotc, zgraztotn, zgraztotf
       REAL(wp) :: zgrarem, zgrafer, zgrapoc, zprcaca, zmortz
       REAL(wp) :: zrespz, ztortz, zgrasrat, zgrasratn
@@ -119,6 +120,13 @@
          zgrazing(ji,jj,jk) = zgraztotc
 
-         !    Various remineralization and excretion terms
-         !    --------------------------------------------
+
+         ! Microzooplankton efficiency. 
+         ! We adopt a formulation proposed by Mitra et al. (2007)
+         ! The gross growth efficiency is controled by the most limiting nutrient.
+         ! Growth is also further decreased when the food quality is poor. This is currently
+         ! hard coded : it can be decreased by up to 50% (zepsherq)
+         ! GGE can also be decreased when food quantity is high, zepsherf (Montagnes and 
+         ! Fulton, 2012)
+         ! -----------------------------------------------------------------------------
          zgrasrat  = ( zgraztotf + rtrn ) / ( zgraztotc + rtrn )
          zgrasratn = ( zgraztotn + rtrn ) / ( zgraztotc + rtrn )
@@ -126,5 +134,6 @@
          zbeta     = MAX(0., (epsher - epshermin) )
          zepsherf  = epshermin + zbeta / ( 1.0 + 0.04E6 * 12. * zfood * zbeta )
-         zepsherv  = zepsherf * zepshert 
+         zepsherq  = 0.5 + (1.0 - 0.5) * zepshert * ( 1.0 + 1.0 ) / ( zepshert + 1.0 )
+         zepsherv  = zepsherf * zepshert * zepsherq 
 
          zgrafer   = zgraztotc * MAX( 0. , ( 1. - unass ) * zgrasrat - ferat3 * zepsherv ) 
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P4Z/p4zsms.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P4Z/p4zsms.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/P4Z/p4zsms.F90	(revision 13193)
@@ -207,5 +207,5 @@
       IF( l_trdtrc ) THEN
          DO jn = jp_pcs0, jp_pcs1
-           ztrdt(:,:,:,jn) = ( tr(:,:,:,jn,Kbb) - ztrdt(:,:,:,jn) ) * rfact2r 
+           ztrdt(:,:,:,jn) = ( tr(:,:,:,jn,Kbb) - ztrdt(:,:,:,jn) ) * rfactr 
            CALL trd_trc( tr(:,:,:,jn,Krhs), jn, jptra_sms, kt, Kmm )   ! save trends
          END DO
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/SED/sedchem.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/SED/sedchem.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/SED/sedchem.F90	(revision 13193)
@@ -577,5 +577,5 @@
          saltprac(:) = salt(:) * 35.0 / 35.16504
       ELSE
-         saltprac(:) = temp(:)
+         saltprac(:) = salt(:)
       ENDIF
 
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/SED/sedinorg.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/SED/sedinorg.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/PISCES/SED/sedinorg.F90	(revision 13193)
@@ -89,4 +89,5 @@
             zsolcpcl = zsolcpcl + solcp(ji,jk,jsclay) * dz(jk)
          END DO
+         zsolcpsi = MAX( zsolcpsi, rtrn )
          zsieq(ji) = sieqs(ji) * MAX(0.25, 1.0 - (0.045 * zsolcpcl / zsolcpsi )**0.58 )
          zsieq(ji) = MAX( rtrn, sieqs(ji) )
Index: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/trcbc.F90
===================================================================
--- NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/trcbc.F90	(revision 13167)
+++ NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/TOP/trcbc.F90	(revision 13193)
@@ -152,5 +152,5 @@
                IF(trcdta_bdy(jn,ib)%cn_obc == 'frs' .AND. nn_trcdmp_bdy(ib) /= 0 )  &
                    & CALL ctl_stop( 'trc_bc_ini: Use FRS OR relaxation' )
-               IF(  .NOT.( 0 < nn_trcdmp_bdy(ib)  .AND.  nn_trcdmp_bdy(ib) <= 2 )  )   &
+               IF(  .NOT.( 0 <= nn_trcdmp_bdy(ib)  .AND.  nn_trcdmp_bdy(ib) <= 2 )  )   &
                    & CALL ctl_stop( 'trc_bc_ini: Not a valid option for nn_trcdmp_bdy. Allowed: 0,1,2.' )
             END DO
