Index: NEMO/branches/2020/r12377_ticket2386/tests/ISOMIP+/MY_SRC/dtatsd.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/tests/ISOMIP+/MY_SRC/dtatsd.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/tests/ISOMIP+/MY_SRC/dtatsd.F90	(revision 13540)
@@ -36,4 +36,6 @@
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_tsddmp ! structure of input SST (file informations, fields read)
 
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -67,8 +69,6 @@
       ierr0 = 0  ;  ierr1 = 0  ;  ierr2 = 0  ;  ierr3 = 0
       !
-      REWIND( numnam_ref )              ! Namelist namtsd in reference namelist : 
       READ  ( numnam_ref, namtsd, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtsd in reference namelist' )
-      REWIND( numnam_cfg )              ! Namelist namtsd in configuration namelist : Parameters of the run
       READ  ( numnam_cfg, namtsd, IOSTAT = ios, ERR = 902 )
 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtsd in configuration namelist' )
@@ -191,32 +191,30 @@
          ENDIF
          !
-         DO jj = 1, jpj                         ! vertical interpolation of T & S
-            DO ji = 1, jpi
-               DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points
-                  zl = gdept_0(ji,jj,jk)
-                  IF(     zl < gdept_1d(1  ) ) THEN          ! above the first level of data
-                     ztp(jk) =  ptsd(ji,jj,1    ,jp_tem)
-                     zsp(jk) =  ptsd(ji,jj,1    ,jp_sal)
-                  ELSEIF( zl > gdept_1d(jpk) ) THEN          ! below the last level of data
-                     ztp(jk) =  ptsd(ji,jj,jpkm1,jp_tem)
-                     zsp(jk) =  ptsd(ji,jj,jpkm1,jp_sal)
-                  ELSE                                      ! inbetween : vertical interpolation between jkk & jkk+1
-                     DO jkk = 1, jpkm1                                  ! when  gdept(jkk) < zl < gdept(jkk+1)
-                        IF( (zl-gdept_1d(jkk)) * (zl-gdept_1d(jkk+1)) <= 0._wp ) THEN
-                           zi = ( zl - gdept_1d(jkk) ) / (gdept_1d(jkk+1)-gdept_1d(jkk))
-                           ztp(jk) = ptsd(ji,jj,jkk,jp_tem) + ( ptsd(ji,jj,jkk+1,jp_tem) - ptsd(ji,jj,jkk,jp_tem) ) * zi 
-                           zsp(jk) = ptsd(ji,jj,jkk,jp_sal) + ( ptsd(ji,jj,jkk+1,jp_sal) - ptsd(ji,jj,jkk,jp_sal) ) * zi
-                        ENDIF
-                     END DO
-                  ENDIF
-               END DO
-               DO jk = 1, jpkm1
-                  ptsd(ji,jj,jk,jp_tem) = ztp(jk) * tmask(ji,jj,jk)     ! mask required for mixed zps-s-coord
-                  ptsd(ji,jj,jk,jp_sal) = zsp(jk) * tmask(ji,jj,jk)
-               END DO
-               ptsd(ji,jj,jpk,jp_tem) = 0._wp
-               ptsd(ji,jj,jpk,jp_sal) = 0._wp
+         DO_2D( 1, 1, 1, 1 )
+            DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points
+               zl = gdept_0(ji,jj,jk)
+               IF(     zl < gdept_1d(1  ) ) THEN          ! above the first level of data
+                  ztp(jk) =  ptsd(ji,jj,1    ,jp_tem)
+                  zsp(jk) =  ptsd(ji,jj,1    ,jp_sal)
+               ELSEIF( zl > gdept_1d(jpk) ) THEN          ! below the last level of data
+                  ztp(jk) =  ptsd(ji,jj,jpkm1,jp_tem)
+                  zsp(jk) =  ptsd(ji,jj,jpkm1,jp_sal)
+               ELSE                                      ! inbetween : vertical interpolation between jkk & jkk+1
+                  DO jkk = 1, jpkm1                                  ! when  gdept(jkk) < zl < gdept(jkk+1)
+                     IF( (zl-gdept_1d(jkk)) * (zl-gdept_1d(jkk+1)) <= 0._wp ) THEN
+                        zi = ( zl - gdept_1d(jkk) ) / (gdept_1d(jkk+1)-gdept_1d(jkk))
+                        ztp(jk) = ptsd(ji,jj,jkk,jp_tem) + ( ptsd(ji,jj,jkk+1,jp_tem) - ptsd(ji,jj,jkk,jp_tem) ) * zi 
+                        zsp(jk) = ptsd(ji,jj,jkk,jp_sal) + ( ptsd(ji,jj,jkk+1,jp_sal) - ptsd(ji,jj,jkk,jp_sal) ) * zi
+                     ENDIF
+                  END DO
+               ENDIF
             END DO
-         END DO
+            DO jk = 1, jpkm1
+               ptsd(ji,jj,jk,jp_tem) = ztp(jk) * tmask(ji,jj,jk)     ! mask required for mixed zps-s-coord
+               ptsd(ji,jj,jk,jp_sal) = zsp(jk) * tmask(ji,jj,jk)
+            END DO
+            ptsd(ji,jj,jpk,jp_tem) = 0._wp
+            ptsd(ji,jj,jpk,jp_sal) = 0._wp
+         END_2D
          ! 
       ELSE                                !==   z- or zps- coordinate   ==!
@@ -226,20 +224,18 @@
          !
          IF( ln_zps ) THEN                      ! zps-coordinate (partial steps) interpolation at the last ocean level
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ik = mbkt(ji,jj) 
-                  IF( ik > 1 ) THEN
-                     zl = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
-                     ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik-1,jp_tem)
-                     ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik-1,jp_sal)
-                  ENDIF
-                  ik = mikt(ji,jj)
-                  IF( ik > 1 ) THEN
-                     zl = ( gdept_0(ji,jj,ik) - gdept_1d(ik) ) / ( gdept_1d(ik+1) - gdept_1d(ik) ) 
-                     ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik+1,jp_tem)
-                     ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik+1,jp_sal)
-                  END IF
-               END DO
-            END DO
+            DO_2D( 1, 1, 1, 1 )
+               ik = mbkt(ji,jj) 
+               IF( ik > 1 ) THEN
+                  zl = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
+                  ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik-1,jp_tem)
+                  ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik-1,jp_sal)
+               ENDIF
+               ik = mikt(ji,jj)
+               IF( ik > 1 ) THEN
+                  zl = ( gdept_0(ji,jj,ik) - gdept_1d(ik) ) / ( gdept_1d(ik+1) - gdept_1d(ik) ) 
+                  ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik+1,jp_tem)
+                  ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik+1,jp_sal)
+               END IF
+            END_2D
          ENDIF
          !
Index: NEMO/branches/2020/r12377_ticket2386/tests/ISOMIP+/MY_SRC/eosbn2.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/tests/ISOMIP+/MY_SRC/eosbn2.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/tests/ISOMIP+/MY_SRC/eosbn2.F90	(revision 13540)
@@ -180,4 +180,6 @@
    REAL(wp) ::   BPE002
 
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -241,12 +243,129 @@
       CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
          !
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  !
-                  zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth
-                  zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature
-                  zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
-                  ztm = tmask(ji,jj,jk)                                         ! tmask
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+            !
+            zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth
+            zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature
+            zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
+            ztm = tmask(ji,jj,jk)                                         ! tmask
+            !
+            zn3 = EOS013*zt   &
+               &   + EOS103*zs+EOS003
+               !
+            zn2 = (EOS022*zt   &
+               &   + EOS112*zs+EOS012)*zt   &
+               &   + (EOS202*zs+EOS102)*zs+EOS002
+               !
+            zn1 = (((EOS041*zt   &
+               &   + EOS131*zs+EOS031)*zt   &
+               &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   &
+               &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   &
+               &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
+               !
+            zn0 = (((((EOS060*zt   &
+               &   + EOS150*zs+EOS050)*zt   &
+               &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   &
+               &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   &
+               &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   &
+               &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   &
+               &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
+               !
+            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+            !
+            prd(ji,jj,jk) = (  zn * r1_rho0 - 1._wp  ) * ztm  ! density anomaly (masked)
+            !
+         END_3D
+         !
+      CASE( np_seos )                !==  simplified EOS  ==!
+         !
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+            zt  = pts  (ji,jj,jk,jp_tem) - 10._wp
+            zs  = pts  (ji,jj,jk,jp_sal) - 35._wp
+            zh  = pdep (ji,jj,jk)
+            ztm = tmask(ji,jj,jk)
+            !
+            zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt   &
+               &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   &
+               &  - rn_nu * zt * zs
+               !                                 
+            prd(ji,jj,jk) = zn * r1_rho0 * ztm                ! density anomaly (masked)
+         END_3D
+         !
+      CASE( np_leos )                !==  linear ISOMIP EOS  ==!
+         !
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+            zt  = pts  (ji,jj,jk,jp_tem) - (-1._wp)
+            zs  = pts  (ji,jj,jk,jp_sal) - 34.2_wp
+            zh  = pdep (ji,jj,jk)
+            ztm = tmask(ji,jj,jk)
+            !
+            zn =  rho0 * ( - rn_a0 * zt + rn_b0 * zs )
+            !                                 
+            prd(ji,jj,jk) = zn * r1_rho0 * ztm                ! density anomaly (masked)
+         END_3D
+         !
+      END SELECT
+      !
+      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-insitu  : ', kdim=jpk )
+      !
+      IF( ln_timing )   CALL timing_stop('eos-insitu')
+      !
+   END SUBROUTINE eos_insitu
+
+
+   SUBROUTINE eos_insitu_pot( pts, prd, prhop, pdep )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE eos_insitu_pot  ***
+      !!
+      !! ** Purpose :   Compute the in situ density (ratio rho/rho0) and the
+      !!      potential volumic mass (Kg/m3) from potential temperature and
+      !!      salinity fields using an equation of state selected in the
+      !!     namelist.
+      !!
+      !! ** Action  : - prd  , the in situ density (no units)
+      !!              - prhop, the potential volumic mass (Kg/m3)
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celsius]
+      !                                                                ! 2 : salinity               [psu]
+      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prd    ! in situ density            [-]
+      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prhop  ! potential density (surface referenced)
+      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pdep   ! depth                      [m]
+      !
+      INTEGER  ::   ji, jj, jk, jsmp             ! dummy loop indices
+      INTEGER  ::   jdof
+      REAL(wp) ::   zt , zh , zstemp, zs , ztm   ! local scalars
+      REAL(wp) ::   zn , zn0, zn1, zn2, zn3      !   -      -
+      REAL(wp), DIMENSION(:), ALLOCATABLE :: zn0_sto, zn_sto, zsign    ! local vectors
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('eos-pot')
+      !
+      SELECT CASE ( neos )
+      !
+      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
+         !
+         ! Stochastic equation of state
+         IF ( ln_sto_eos ) THEN
+            ALLOCATE(zn0_sto(1:2*nn_sto_eos))
+            ALLOCATE(zn_sto(1:2*nn_sto_eos))
+            ALLOCATE(zsign(1:2*nn_sto_eos))
+            DO jsmp = 1, 2*nn_sto_eos, 2
+              zsign(jsmp)   = 1._wp
+              zsign(jsmp+1) = -1._wp
+            END DO
+            !
+            DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+               !
+               ! compute density (2*nn_sto_eos) times:
+               ! (1) for t+dt, s+ds (with the random TS fluctutation computed in sto_pts)
+               ! (2) for t-dt, s-ds (with the opposite fluctuation)
+               DO jsmp = 1, nn_sto_eos*2
+                  jdof   = (jsmp + 1) / 2
+                  zh     = pdep(ji,jj,jk) * r1_Z0                                  ! depth
+                  zt     = (pts (ji,jj,jk,jp_tem) + pts_ran(ji,jj,jk,jp_tem,jdof) * zsign(jsmp)) * r1_T0    ! temperature
+                  zstemp = pts  (ji,jj,jk,jp_sal) + pts_ran(ji,jj,jk,jp_sal,jdof) * zsign(jsmp)
+                  zs     = SQRT( ABS( zstemp + rdeltaS ) * r1_S0 )   ! square root salinity
+                  ztm    = tmask(ji,jj,jk)                                         ! tmask
                   !
                   zn3 = EOS013*zt   &
@@ -263,5 +382,5 @@
                      &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
                      !
-                  zn0 = (((((EOS060*zt   &
+                  zn0_sto(jsmp) = (((((EOS060*zt   &
                      &   + EOS150*zs+EOS050)*zt   &
                      &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   &
@@ -271,278 +390,26 @@
                      &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
                      !
-                  zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+                  zn_sto(jsmp)  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0_sto(jsmp)
+               END DO
+               !
+               ! compute stochastic density as the mean of the (2*nn_sto_eos) densities
+               prhop(ji,jj,jk) = 0._wp ; prd(ji,jj,jk) = 0._wp
+               DO jsmp = 1, nn_sto_eos*2
+                  prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp)                      ! potential density referenced at the surface
                   !
-                  prd(ji,jj,jk) = (  zn * r1_rho0 - 1._wp  ) * ztm  ! density anomaly (masked)
-                  !
+                  prd(ji,jj,jk) = prd(ji,jj,jk) + (  zn_sto(jsmp) * r1_rho0 - 1._wp  )   ! density anomaly (masked)
                END DO
-            END DO
-         END DO
-         !
-      CASE( np_seos )                !==  simplified EOS  ==!
-         !
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zt  = pts  (ji,jj,jk,jp_tem) - 10._wp
-                  zs  = pts  (ji,jj,jk,jp_sal) - 35._wp
-                  zh  = pdep (ji,jj,jk)
-                  ztm = tmask(ji,jj,jk)
-                  !
-                  zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt   &
-                     &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   &
-                     &  - rn_nu * zt * zs
-                     !                                 
-                  prd(ji,jj,jk) = zn * r1_rho0 * ztm                ! density anomaly (masked)
-               END DO
-            END DO
-         END DO
-         !
-      CASE( np_leos )                !==  linear ISOMIP EOS  ==!
-         !
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zt  = pts  (ji,jj,jk,jp_tem) - (-1._wp)
-                  zs  = pts  (ji,jj,jk,jp_sal) - 34.2_wp
-                  zh  = pdep (ji,jj,jk)
-                  ztm = tmask(ji,jj,jk)
-                  !
-                  zn =  rho0 * ( - rn_a0 * zt + rn_b0 * zs )
-                  !                                 
-                  prd(ji,jj,jk) = zn * r1_rho0 * ztm                ! density anomaly (masked)
-               END DO
-            END DO
-         END DO
-         !
-      END SELECT
-      !
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-insitu  : ', kdim=jpk )
-      !
-      IF( ln_timing )   CALL timing_stop('eos-insitu')
-      !
-   END SUBROUTINE eos_insitu
-
-
-   SUBROUTINE eos_insitu_pot( pts, prd, prhop, pdep )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE eos_insitu_pot  ***
-      !!
-      !! ** Purpose :   Compute the in situ density (ratio rho/rho0) and the
-      !!      potential volumic mass (Kg/m3) from potential temperature and
-      !!      salinity fields using an equation of state selected in the
-      !!     namelist.
-      !!
-      !! ** Action  : - prd  , the in situ density (no units)
-      !!              - prhop, the potential volumic mass (Kg/m3)
-      !!
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celsius]
-      !                                                                ! 2 : salinity               [psu]
-      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prd    ! in situ density            [-]
-      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(  out) ::   prhop  ! potential density (surface referenced)
-      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pdep   ! depth                      [m]
-      !
-      INTEGER  ::   ji, jj, jk, jsmp             ! dummy loop indices
-      INTEGER  ::   jdof
-      REAL(wp) ::   zt , zh , zstemp, zs , ztm   ! local scalars
-      REAL(wp) ::   zn , zn0, zn1, zn2, zn3      !   -      -
-      REAL(wp), DIMENSION(:), ALLOCATABLE :: zn0_sto, zn_sto, zsign    ! local vectors
-      !!----------------------------------------------------------------------
-      !
-      IF( ln_timing )   CALL timing_start('eos-pot')
-      !
-      SELECT CASE ( neos )
-      !
-      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
-         !
-         ! Stochastic equation of state
-         IF ( ln_sto_eos ) THEN
-            ALLOCATE(zn0_sto(1:2*nn_sto_eos))
-            ALLOCATE(zn_sto(1:2*nn_sto_eos))
-            ALLOCATE(zsign(1:2*nn_sto_eos))
-            DO jsmp = 1, 2*nn_sto_eos, 2
-              zsign(jsmp)   = 1._wp
-              zsign(jsmp+1) = -1._wp
-            END DO
-            !
-            DO jk = 1, jpkm1
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     !
-                     ! compute density (2*nn_sto_eos) times:
-                     ! (1) for t+dt, s+ds (with the random TS fluctutation computed in sto_pts)
-                     ! (2) for t-dt, s-ds (with the opposite fluctuation)
-                     DO jsmp = 1, nn_sto_eos*2
-                        jdof   = (jsmp + 1) / 2
-                        zh     = pdep(ji,jj,jk) * r1_Z0                                  ! depth
-                        zt     = (pts (ji,jj,jk,jp_tem) + pts_ran(ji,jj,jk,jp_tem,jdof) * zsign(jsmp)) * r1_T0    ! temperature
-                        zstemp = pts  (ji,jj,jk,jp_sal) + pts_ran(ji,jj,jk,jp_sal,jdof) * zsign(jsmp)
-                        zs     = SQRT( ABS( zstemp + rdeltaS ) * r1_S0 )   ! square root salinity
-                        ztm    = tmask(ji,jj,jk)                                         ! tmask
-                        !
-                        zn3 = EOS013*zt   &
-                           &   + EOS103*zs+EOS003
-                           !
-                        zn2 = (EOS022*zt   &
-                           &   + EOS112*zs+EOS012)*zt   &
-                           &   + (EOS202*zs+EOS102)*zs+EOS002
-                           !
-                        zn1 = (((EOS041*zt   &
-                           &   + EOS131*zs+EOS031)*zt   &
-                           &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   &
-                           &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   &
-                           &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
-                           !
-                        zn0_sto(jsmp) = (((((EOS060*zt   &
-                           &   + EOS150*zs+EOS050)*zt   &
-                           &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   &
-                           &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   &
-                           &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   &
-                           &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   &
-                           &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
-                           !
-                        zn_sto(jsmp)  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0_sto(jsmp)
-                     END DO
-                     !
-                     ! compute stochastic density as the mean of the (2*nn_sto_eos) densities
-                     prhop(ji,jj,jk) = 0._wp ; prd(ji,jj,jk) = 0._wp
-                     DO jsmp = 1, nn_sto_eos*2
-                        prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp)                      ! potential density referenced at the surface
-                        !
-                        prd(ji,jj,jk) = prd(ji,jj,jk) + (  zn_sto(jsmp) * r1_rho0 - 1._wp  )   ! density anomaly (masked)
-                     END DO
-                     prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos
-                     prd  (ji,jj,jk) = 0.5_wp * prd  (ji,jj,jk) * ztm / nn_sto_eos
-                  END DO
-               END DO
-            END DO
+               prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos
+               prd  (ji,jj,jk) = 0.5_wp * prd  (ji,jj,jk) * ztm / nn_sto_eos
+            END_3D
             DEALLOCATE(zn0_sto,zn_sto,zsign)
          ! Non-stochastic equation of state
          ELSE
-            DO jk = 1, jpkm1
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     !
-                     zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth
-                     zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature
-                     zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
-                     ztm = tmask(ji,jj,jk)                                         ! tmask
-                     !
-                     zn3 = EOS013*zt   &
-                        &   + EOS103*zs+EOS003
-                        !
-                     zn2 = (EOS022*zt   &
-                        &   + EOS112*zs+EOS012)*zt   &
-                        &   + (EOS202*zs+EOS102)*zs+EOS002
-                        !
-                     zn1 = (((EOS041*zt   &
-                        &   + EOS131*zs+EOS031)*zt   &
-                        &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   &
-                        &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   &
-                        &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
-                        !
-                     zn0 = (((((EOS060*zt   &
-                        &   + EOS150*zs+EOS050)*zt   &
-                        &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   &
-                        &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   &
-                        &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   &
-                        &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   &
-                        &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
-                        !
-                     zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
-                     !
-                     prhop(ji,jj,jk) = zn0 * ztm                           ! potential density referenced at the surface
-                     !
-                     prd(ji,jj,jk) = (  zn * r1_rho0 - 1._wp  ) * ztm      ! density anomaly (masked)
-                  END DO
-               END DO
-            END DO
-         ENDIF
-         
-      CASE( np_seos )                !==  simplified EOS  ==!
-         !
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zt  = pts  (ji,jj,jk,jp_tem) - 10._wp
-                  zs  = pts  (ji,jj,jk,jp_sal) - 35._wp
-                  zh  = pdep (ji,jj,jk)
-                  ztm = tmask(ji,jj,jk)
-                  !                                                     ! potential density referenced at the surface
-                  zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt ) * zt   &
-                     &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs ) * zs   &
-                     &  - rn_nu * zt * zs
-                  prhop(ji,jj,jk) = ( rho0 + zn ) * ztm
-                  !                                                     ! density anomaly (masked)
-                  zn = zn - ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zh
-                  prd(ji,jj,jk) = zn * r1_rho0 * ztm
-                  !
-               END DO
-            END DO
-         END DO
-         !
-      CASE( np_leos )                !==  linear ISOMIP EOS  ==!
-         !
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zt  = pts  (ji,jj,jk,jp_tem) - (-1._wp)
-                  zs  = pts  (ji,jj,jk,jp_sal) - 34.2_wp
-                  zh  = pdep (ji,jj,jk)
-                  ztm = tmask(ji,jj,jk)
-                  !                                                     ! potential density referenced at the surface
-                  zn =  rho0 * ( - rn_a0 * zt + rn_b0 * zs )
-                  prhop(ji,jj,jk) = ( rho0 + zn ) * ztm
-                  !                                                     ! density anomaly (masked)
-                  prd(ji,jj,jk) = zn * r1_rho0 * ztm
-                  !
-               END DO
-            END DO
-         END DO
-         !
-      END SELECT
-      !
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-pot: ', tab3d_2=prhop, clinfo2=' pot : ', kdim=jpk )
-      !
-      IF( ln_timing )   CALL timing_stop('eos-pot')
-      !
-   END SUBROUTINE eos_insitu_pot
-
-
-   SUBROUTINE eos_insitu_2d( pts, pdep, prd )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE eos_insitu_2d  ***
-      !!
-      !! ** Purpose :   Compute the in situ density (ratio rho/rho0) from
-      !!      potential temperature and salinity using an equation of state
-      !!      selected in the nameos namelist. * 2D field case
-      !!
-      !! ** Action  : - prd , the in situ density (no units) (unmasked)
-      !!
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celsius]
-      !                                                           ! 2 : salinity               [psu]
-      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(in   ) ::   pdep  ! depth                      [m]
-      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(  out) ::   prd   ! in situ density
-      !
-      INTEGER  ::   ji, jj, jk                ! dummy loop indices
-      REAL(wp) ::   zt , zh , zs              ! local scalars
-      REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      -
-      !!----------------------------------------------------------------------
-      !
-      IF( ln_timing )   CALL timing_start('eos2d')
-      !
-      prd(:,:) = 0._wp
-      !
-      SELECT CASE( neos )
-      !
-      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
-         !
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               !
-               zh  = pdep(ji,jj) * r1_Z0                                  ! depth
-               zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature
-               zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
+            DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+               !
+               zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth
+               zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature
+               zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
+               ztm = tmask(ji,jj,jk)                                         ! tmask
                !
                zn3 = EOS013*zt   &
@@ -569,52 +436,148 @@
                zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
                !
-               prd(ji,jj) = zn * r1_rho0 - 1._wp               ! unmasked in situ density anomaly
-               !
-            END DO
-         END DO
-         !
-         CALL lbc_lnk( 'eosbn2', prd, 'T', 1. )                    ! Lateral boundary conditions
-         !
+               prhop(ji,jj,jk) = zn0 * ztm                           ! potential density referenced at the surface
+               !
+               prd(ji,jj,jk) = (  zn * r1_rho0 - 1._wp  ) * ztm      ! density anomaly (masked)
+            END_3D
+         ENDIF
+         
       CASE( np_seos )                !==  simplified EOS  ==!
          !
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               !
-               zt    = pts  (ji,jj,jp_tem)  - 10._wp
-               zs    = pts  (ji,jj,jp_sal)  - 35._wp
-               zh    = pdep (ji,jj)                         ! depth at the partial step level
-               !
-               zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt   &
-                  &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   &
-                  &  - rn_nu * zt * zs
-                  !
-               prd(ji,jj) = zn * r1_rho0               ! unmasked in situ density anomaly
-               !
-            END DO
-         END DO
-         !
-         CALL lbc_lnk( 'eosbn2', prd, 'T', 1. )                    ! Lateral boundary conditions
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+            zt  = pts  (ji,jj,jk,jp_tem) - 10._wp
+            zs  = pts  (ji,jj,jk,jp_sal) - 35._wp
+            zh  = pdep (ji,jj,jk)
+            ztm = tmask(ji,jj,jk)
+            !                                                     ! potential density referenced at the surface
+            zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt ) * zt   &
+               &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs ) * zs   &
+               &  - rn_nu * zt * zs
+            prhop(ji,jj,jk) = ( rho0 + zn ) * ztm
+            !                                                     ! density anomaly (masked)
+            zn = zn - ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zh
+            prd(ji,jj,jk) = zn * r1_rho0 * ztm
+            !
+         END_3D
+         !
+      CASE( np_leos )                !==  linear ISOMIP EOS  ==!
+         !
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+            zt  = pts  (ji,jj,jk,jp_tem) - (-1._wp)
+            zs  = pts  (ji,jj,jk,jp_sal) - 34.2_wp
+            zh  = pdep (ji,jj,jk)
+            ztm = tmask(ji,jj,jk)
+            !                                                     ! potential density referenced at the surface
+            zn =  rho0 * ( - rn_a0 * zt + rn_b0 * zs )
+            prhop(ji,jj,jk) = ( rho0 + zn ) * ztm
+            !                                                     ! density anomaly (masked)
+            prd(ji,jj,jk) = zn * r1_rho0 * ztm
+            !
+         END_3D
+         !
+      END SELECT
+      !
+      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-pot: ', tab3d_2=prhop, clinfo2=' pot : ', kdim=jpk )
+      !
+      IF( ln_timing )   CALL timing_stop('eos-pot')
+      !
+   END SUBROUTINE eos_insitu_pot
+
+
+   SUBROUTINE eos_insitu_2d( pts, pdep, prd )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE eos_insitu_2d  ***
+      !!
+      !! ** Purpose :   Compute the in situ density (ratio rho/rho0) from
+      !!      potential temperature and salinity using an equation of state
+      !!      selected in the nameos namelist. * 2D field case
+      !!
+      !! ** Action  : - prd , the in situ density (no units) (unmasked)
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) ::   pts   ! 1 : potential temperature  [Celsius]
+      !                                                           ! 2 : salinity               [psu]
+      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(in   ) ::   pdep  ! depth                      [m]
+      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(  out) ::   prd   ! in situ density
+      !
+      INTEGER  ::   ji, jj, jk                ! dummy loop indices
+      REAL(wp) ::   zt , zh , zs              ! local scalars
+      REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      -
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('eos2d')
+      !
+      prd(:,:) = 0._wp
+      !
+      SELECT CASE( neos )
+      !
+      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
+         !
+         DO_2D( 1, 1, 1, 1 )
+            !
+            zh  = pdep(ji,jj) * r1_Z0                                  ! depth
+            zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature
+            zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
+            !
+            zn3 = EOS013*zt   &
+               &   + EOS103*zs+EOS003
+               !
+            zn2 = (EOS022*zt   &
+               &   + EOS112*zs+EOS012)*zt   &
+               &   + (EOS202*zs+EOS102)*zs+EOS002
+               !
+            zn1 = (((EOS041*zt   &
+               &   + EOS131*zs+EOS031)*zt   &
+               &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   &
+               &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   &
+               &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001
+               !
+            zn0 = (((((EOS060*zt   &
+               &   + EOS150*zs+EOS050)*zt   &
+               &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   &
+               &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   &
+               &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   &
+               &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   &
+               &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
+               !
+            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+            !
+            prd(ji,jj) = zn * r1_rho0 - 1._wp               ! unmasked in situ density anomaly
+            !
+         END_2D
+         !
+      CASE( np_seos )                !==  simplified EOS  ==!
+         !
+         DO_2D( 1, 1, 1, 1 )
+            !
+            zt    = pts  (ji,jj,jp_tem)  - 10._wp
+            zs    = pts  (ji,jj,jp_sal)  - 35._wp
+            zh    = pdep (ji,jj)                         ! depth at the partial step level
+            !
+            zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt   &
+               &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs   &
+               &  - rn_nu * zt * zs
+               !
+            prd(ji,jj) = zn * r1_rho0               ! unmasked in situ density anomaly
+            !
+         END_2D
          !
       CASE( np_leos )                !==  ISOMIP EOS  ==!
          !
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               !
-               zt    = pts  (ji,jj,jp_tem)  - (-1._wp)
-               zs    = pts  (ji,jj,jp_sal)  - 34.2_wp
-               zh    = pdep (ji,jj)                         ! depth at the partial step level
-               !
-               zn =  rho0 * ( - rn_a0 * zt + rn_b0 * zs )
-                  !
-               prd(ji,jj) = zn * r1_rho0               ! unmasked in situ density anomaly
-               !
-            END DO
-         END DO
-         !
-         CALL lbc_lnk( 'eosbn2', prd, 'T', 1. )                    ! Lateral boundary conditions
+         DO_2D( 1, 1, 1, 1 )
+            !
+            zt    = pts  (ji,jj,jp_tem)  - (-1._wp)
+            zs    = pts  (ji,jj,jp_sal)  - 34.2_wp
+            zh    = pdep (ji,jj)                         ! depth at the partial step level
+            !
+            zn =  rho0 * ( - rn_a0 * zt + rn_b0 * zs )
+            !
+            prd(ji,jj) = zn * r1_rho0               ! unmasked in situ density anomaly
+            !
+         END_2D
+         !
          !
       END SELECT
       !
-      IF(ln_ctl)   CALL prt_ctl( tab2d_1=prd, clinfo1=' eos2d: ' )
+      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=prd, clinfo1=' eos2d: ' )
       !
       IF( ln_timing )   CALL timing_stop('eos2d')
@@ -648,98 +611,86 @@
       CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
          !
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  !
-                  zh  = gdept(ji,jj,jk,Kmm) * r1_Z0                                ! depth
-                  zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature
-                  zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
-                  ztm = tmask(ji,jj,jk)                                         ! tmask
-                  !
-                  ! alpha
-                  zn3 = ALP003
-                  !
-                  zn2 = ALP012*zt + ALP102*zs+ALP002
-                  !
-                  zn1 = ((ALP031*zt   &
-                     &   + ALP121*zs+ALP021)*zt   &
-                     &   + (ALP211*zs+ALP111)*zs+ALP011)*zt   &
-                     &   + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001
-                     !
-                  zn0 = ((((ALP050*zt   &
-                     &   + ALP140*zs+ALP040)*zt   &
-                     &   + (ALP230*zs+ALP130)*zs+ALP030)*zt   &
-                     &   + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt   &
-                     &   + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt   &
-                     &   + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000
-                     !
-                  zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
-                  !
-                  pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm
-                  !
-                  ! beta
-                  zn3 = BET003
-                  !
-                  zn2 = BET012*zt + BET102*zs+BET002
-                  !
-                  zn1 = ((BET031*zt   &
-                     &   + BET121*zs+BET021)*zt   &
-                     &   + (BET211*zs+BET111)*zs+BET011)*zt   &
-                     &   + ((BET301*zs+BET201)*zs+BET101)*zs+BET001
-                     !
-                  zn0 = ((((BET050*zt   &
-                     &   + BET140*zs+BET040)*zt   &
-                     &   + (BET230*zs+BET130)*zs+BET030)*zt   &
-                     &   + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt   &
-                     &   + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt   &
-                     &   + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000
-                     !
-                  zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
-                  !
-                  pab(ji,jj,jk,jp_sal) = zn / zs * r1_rho0 * ztm
-                  !
-               END DO
-            END DO
-         END DO
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+            !
+            zh  = gdept(ji,jj,jk,Kmm) * r1_Z0                                ! depth
+            zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature
+            zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
+            ztm = tmask(ji,jj,jk)                                         ! tmask
+            !
+            ! alpha
+            zn3 = ALP003
+            !
+            zn2 = ALP012*zt + ALP102*zs+ALP002
+            !
+            zn1 = ((ALP031*zt   &
+               &   + ALP121*zs+ALP021)*zt   &
+               &   + (ALP211*zs+ALP111)*zs+ALP011)*zt   &
+               &   + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001
+               !
+            zn0 = ((((ALP050*zt   &
+               &   + ALP140*zs+ALP040)*zt   &
+               &   + (ALP230*zs+ALP130)*zs+ALP030)*zt   &
+               &   + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt   &
+               &   + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt   &
+               &   + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000
+               !
+            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+            !
+            pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm
+            !
+            ! beta
+            zn3 = BET003
+            !
+            zn2 = BET012*zt + BET102*zs+BET002
+            !
+            zn1 = ((BET031*zt   &
+               &   + BET121*zs+BET021)*zt   &
+               &   + (BET211*zs+BET111)*zs+BET011)*zt   &
+               &   + ((BET301*zs+BET201)*zs+BET101)*zs+BET001
+               !
+            zn0 = ((((BET050*zt   &
+               &   + BET140*zs+BET040)*zt   &
+               &   + (BET230*zs+BET130)*zs+BET030)*zt   &
+               &   + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt   &
+               &   + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt   &
+               &   + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000
+               !
+            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+            !
+            pab(ji,jj,jk,jp_sal) = zn / zs * r1_rho0 * ztm
+            !
+         END_3D
          !
       CASE( np_seos )                  !==  simplified EOS  ==!
          !
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zt  = pts (ji,jj,jk,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0)
-                  zs  = pts (ji,jj,jk,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0)
-                  zh  = gdept(ji,jj,jk,Kmm)                ! depth in meters at t-point
-                  ztm = tmask(ji,jj,jk)                  ! land/sea bottom mask = surf. mask
-                  !
-                  zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs
-                  pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm   ! alpha
-                  !
-                  zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt
-                  pab(ji,jj,jk,jp_sal) = zn * r1_rho0 * ztm   ! beta
-                  !
-               END DO
-            END DO
-         END DO
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+            zt  = pts (ji,jj,jk,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0)
+            zs  = pts (ji,jj,jk,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0)
+            zh  = gdept(ji,jj,jk,Kmm)                ! depth in meters at t-point
+            ztm = tmask(ji,jj,jk)                  ! land/sea bottom mask = surf. mask
+            !
+            zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs
+            pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm   ! alpha
+            !
+            zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt
+            pab(ji,jj,jk,jp_sal) = zn * r1_rho0 * ztm   ! beta
+            !
+         END_3D
          !
       CASE( np_leos )                  !==  linear ISOMIP EOS  ==!
          !
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zt  = pts (ji,jj,jk,jp_tem) - (-1._wp)
-                  zs  = pts (ji,jj,jk,jp_sal) - 34.2_wp   ! abs. salinity anomaly (s-S0)
-                  zh  = gdept(ji,jj,jk,Kmm)                 ! depth in meters at t-point
-                  ztm = tmask(ji,jj,jk)                   ! land/sea bottom mask = surf. mask
-                  !
-                  zn  = rn_a0 * rho0
-                  pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm   ! alpha
-                  !
-                  zn  = rn_b0 * rho0
-                  pab(ji,jj,jk,jp_sal) = zn * r1_rho0 * ztm   ! beta
-                  !
-               END DO
-            END DO
-         END DO
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+            zt  = pts (ji,jj,jk,jp_tem) - (-1._wp)
+            zs  = pts (ji,jj,jk,jp_sal) - 34.2_wp   ! abs. salinity anomaly (s-S0)
+            zh  = gdept(ji,jj,jk,Kmm)                 ! depth in meters at t-point
+            ztm = tmask(ji,jj,jk)                   ! land/sea bottom mask = surf. mask
+            !
+            zn  = rn_a0 * rho0
+            pab(ji,jj,jk,jp_tem) = zn * r1_rho0 * ztm   ! alpha
+            !
+            zn  = rn_b0 * rho0
+            pab(ji,jj,jk,jp_sal) = zn * r1_rho0 * ztm   ! beta
+            !
+         END_3D
          !
       CASE DEFAULT
@@ -749,6 +700,6 @@
       END SELECT
       !
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pab(:,:,:,jp_tem), clinfo1=' rab_3d_t: ', &
-         &                       tab3d_2=pab(:,:,:,jp_sal), clinfo2=' rab_3d_s : ', kdim=jpk )
+      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pab(:,:,:,jp_tem), clinfo1=' rab_3d_t: ', &
+         &                                  tab3d_2=pab(:,:,:,jp_sal), clinfo2=' rab_3d_s : ', kdim=jpk )
       !
       IF( ln_timing )   CALL timing_stop('rab_3d')
@@ -783,98 +734,86 @@
       CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
          !
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               !
-               zh  = pdep(ji,jj) * r1_Z0                                  ! depth
-               zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature
-               zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
-               !
-               ! alpha
-               zn3 = ALP003
-               !
-               zn2 = ALP012*zt + ALP102*zs+ALP002
-               !
-               zn1 = ((ALP031*zt   &
-                  &   + ALP121*zs+ALP021)*zt   &
-                  &   + (ALP211*zs+ALP111)*zs+ALP011)*zt   &
-                  &   + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001
-                  !
-               zn0 = ((((ALP050*zt   &
-                  &   + ALP140*zs+ALP040)*zt   &
-                  &   + (ALP230*zs+ALP130)*zs+ALP030)*zt   &
-                  &   + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt   &
-                  &   + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt   &
-                  &   + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000
-                  !
-               zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
-               !
-               pab(ji,jj,jp_tem) = zn * r1_rho0
-               !
-               ! beta
-               zn3 = BET003
-               !
-               zn2 = BET012*zt + BET102*zs+BET002
-               !
-               zn1 = ((BET031*zt   &
-                  &   + BET121*zs+BET021)*zt   &
-                  &   + (BET211*zs+BET111)*zs+BET011)*zt   &
-                  &   + ((BET301*zs+BET201)*zs+BET101)*zs+BET001
-                  !
-               zn0 = ((((BET050*zt   &
-                  &   + BET140*zs+BET040)*zt   &
-                  &   + (BET230*zs+BET130)*zs+BET030)*zt   &
-                  &   + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt   &
-                  &   + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt   &
-                  &   + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000
-                  !
-               zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
-               !
-               pab(ji,jj,jp_sal) = zn / zs * r1_rho0
-               !
-               !
-            END DO
-         END DO
-         !                            ! Lateral boundary conditions
-         CALL lbc_lnk_multi( 'eosbn2', pab(:,:,jp_tem), 'T', 1. , pab(:,:,jp_sal), 'T', 1. )                    
+         DO_2D( 1, 1, 1, 1 )
+            !
+            zh  = pdep(ji,jj) * r1_Z0                                  ! depth
+            zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature
+            zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
+            !
+            ! alpha
+            zn3 = ALP003
+            !
+            zn2 = ALP012*zt + ALP102*zs+ALP002
+            !
+            zn1 = ((ALP031*zt   &
+               &   + ALP121*zs+ALP021)*zt   &
+               &   + (ALP211*zs+ALP111)*zs+ALP011)*zt   &
+               &   + ((ALP301*zs+ALP201)*zs+ALP101)*zs+ALP001
+               !
+            zn0 = ((((ALP050*zt   &
+               &   + ALP140*zs+ALP040)*zt   &
+               &   + (ALP230*zs+ALP130)*zs+ALP030)*zt   &
+               &   + ((ALP320*zs+ALP220)*zs+ALP120)*zs+ALP020)*zt   &
+               &   + (((ALP410*zs+ALP310)*zs+ALP210)*zs+ALP110)*zs+ALP010)*zt   &
+               &   + ((((ALP500*zs+ALP400)*zs+ALP300)*zs+ALP200)*zs+ALP100)*zs+ALP000
+               !
+            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+            !
+            pab(ji,jj,jp_tem) = zn * r1_rho0
+            !
+            ! beta
+            zn3 = BET003
+            !
+            zn2 = BET012*zt + BET102*zs+BET002
+            !
+            zn1 = ((BET031*zt   &
+               &   + BET121*zs+BET021)*zt   &
+               &   + (BET211*zs+BET111)*zs+BET011)*zt   &
+               &   + ((BET301*zs+BET201)*zs+BET101)*zs+BET001
+               !
+            zn0 = ((((BET050*zt   &
+               &   + BET140*zs+BET040)*zt   &
+               &   + (BET230*zs+BET130)*zs+BET030)*zt   &
+               &   + ((BET320*zs+BET220)*zs+BET120)*zs+BET020)*zt   &
+               &   + (((BET410*zs+BET310)*zs+BET210)*zs+BET110)*zs+BET010)*zt   &
+               &   + ((((BET500*zs+BET400)*zs+BET300)*zs+BET200)*zs+BET100)*zs+BET000
+               !
+            zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
+            !
+            pab(ji,jj,jp_sal) = zn / zs * r1_rho0
+            !
+            !
+         END_2D
          !
       CASE( np_seos )                  !==  simplified EOS  ==!
          !
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               !
-               zt    = pts  (ji,jj,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0)
-               zs    = pts  (ji,jj,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0)
-               zh    = pdep (ji,jj)                   ! depth at the partial step level
-               !
-               zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs
-               pab(ji,jj,jp_tem) = zn * r1_rho0   ! alpha
-               !
-               zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt
-               pab(ji,jj,jp_sal) = zn * r1_rho0   ! beta
-               !
-            END DO
-         END DO
-         !                            ! Lateral boundary conditions
-         CALL lbc_lnk_multi( 'eosbn2', pab(:,:,jp_tem), 'T', 1. , pab(:,:,jp_sal), 'T', 1. )                    
+         DO_2D( 1, 1, 1, 1 )
+            !
+            zt    = pts  (ji,jj,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0)
+            zs    = pts  (ji,jj,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0)
+            zh    = pdep (ji,jj)                   ! depth at the partial step level
+            !
+            zn  = rn_a0 * ( 1._wp + rn_lambda1*zt + rn_mu1*zh ) + rn_nu*zs
+            pab(ji,jj,jp_tem) = zn * r1_rho0   ! alpha
+            !
+            zn  = rn_b0 * ( 1._wp - rn_lambda2*zs - rn_mu2*zh ) - rn_nu*zt
+            pab(ji,jj,jp_sal) = zn * r1_rho0   ! beta
+            !
+         END_2D
          !
       CASE( np_leos )                  !==  linear ISOMIP EOS  ==!
          !
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               !
-               zt    = pts  (ji,jj,jp_tem) - (-1._wp)   ! pot. temperature anomaly (t-T0)
-               zs    = pts  (ji,jj,jp_sal) - 34.2_wp   ! abs. salinity anomaly (s-S0)
-               zh    = pdep (ji,jj)                   ! depth at the partial step level
-               !
-               zn  = rn_a0 * rho0
-               pab(ji,jj,jp_tem) = zn * r1_rho0   ! alpha
-               !
-               zn  = rn_b0 * rho0
-               pab(ji,jj,jp_sal) = zn * r1_rho0   ! beta
-               !
-            END DO
-         END DO
-         !
-         CALL lbc_lnk_multi( 'eosbn2', pab(:,:,jp_tem), 'T', 1. , pab(:,:,jp_sal), 'T', 1. )                    ! Lateral boundary conditions
+         DO_2D( 1, 1, 1, 1 )
+            !
+            zt    = pts  (ji,jj,jp_tem) - (-1._wp)   ! pot. temperature anomaly (t-T0)
+            zs    = pts  (ji,jj,jp_sal) - 34.2_wp   ! abs. salinity anomaly (s-S0)
+            zh    = pdep (ji,jj)                   ! depth at the partial step level
+            !
+            zn  = rn_a0 * rho0
+            pab(ji,jj,jp_tem) = zn * r1_rho0   ! alpha
+            !
+            zn  = rn_b0 * rho0
+            pab(ji,jj,jp_sal) = zn * r1_rho0   ! beta
+            !
+         END_2D
          !
       CASE DEFAULT
@@ -884,6 +823,6 @@
       END SELECT
       !
-      IF(ln_ctl)   CALL prt_ctl( tab2d_1=pab(:,:,jp_tem), clinfo1=' rab_2d_t: ', &
-         &                       tab2d_2=pab(:,:,jp_sal), clinfo2=' rab_2d_s : ' )
+      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=pab(:,:,jp_tem), clinfo1=' rab_2d_t: ', &
+         &                                  tab2d_2=pab(:,:,jp_sal), clinfo2=' rab_2d_s : ' )
       !
       IF( ln_timing )   CALL timing_stop('rab_2d')
@@ -1026,21 +965,17 @@
       IF( ln_timing )   CALL timing_start('bn2')
       !
-      DO jk = 2, jpkm1           ! interior points only (2=< jk =< jpkm1 )
-         DO jj = 1, jpj          ! surface and bottom value set to zero one for all in istate.F90
-            DO ji = 1, jpi
-               zrw =   ( gdepw(ji,jj,jk  ,Kmm) - gdept(ji,jj,jk,Kmm) )   &
-                  &  / ( gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm) ) 
-                  !
-               zaw = pab(ji,jj,jk,jp_tem) * (1. - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw 
-               zbw = pab(ji,jj,jk,jp_sal) * (1. - zrw) + pab(ji,jj,jk-1,jp_sal) * zrw
-               !
-               pn2(ji,jj,jk) = grav * (  zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )     &
-                  &                    - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )  )  &
-                  &            / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-      !
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2  : ', kdim=jpk )
+      DO_3D( 1, 1, 1, 1, 2, jpkm1 )
+         zrw =   ( gdepw(ji,jj,jk  ,Kmm) - gdept(ji,jj,jk,Kmm) )   &
+            &  / ( gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm) ) 
+            !
+         zaw = pab(ji,jj,jk,jp_tem) * (1. - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw 
+         zbw = pab(ji,jj,jk,jp_sal) * (1. - zrw) + pab(ji,jj,jk-1,jp_sal) * zrw
+         !
+         pn2(ji,jj,jk) = grav * (  zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )     &
+            &                    - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )  )  &
+            &            / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)
+      END_3D
+      !
+      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2  : ', kdim=jpk )
       !
       IF( ln_timing )   CALL timing_stop('bn2')
@@ -1078,28 +1013,26 @@
       z1_T0   = 1._wp/40._wp
       !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            !
-            zt  = ctmp   (ji,jj) * z1_T0
-            zs  = SQRT( ABS( psal(ji,jj) + zdeltaS ) * r1_S0 )
-            ztm = tmask(ji,jj,1)
-            !
-            zn = ((((-2.1385727895e-01_wp*zt   &
-               &   - 2.7674419971e-01_wp*zs+1.0728094330_wp)*zt   &
-               &   + (2.6366564313_wp*zs+3.3546960647_wp)*zs-7.8012209473_wp)*zt   &
-               &   + ((1.8835586562_wp*zs+7.3949191679_wp)*zs-3.3937395875_wp)*zs-5.6414948432_wp)*zt   &
-               &   + (((3.5737370589_wp*zs-1.5512427389e+01_wp)*zs+2.4625741105e+01_wp)*zs   &
-               &      +1.9912291000e+01_wp)*zs-3.2191146312e+01_wp)*zt   &
-               &   + ((((5.7153204649e-01_wp*zs-3.0943149543_wp)*zs+9.3052495181_wp)*zs   &
-               &      -9.4528934807_wp)*zs+3.1066408996_wp)*zs-4.3504021262e-01_wp
-               !
-            zd = (2.0035003456_wp*zt   &
-               &   -3.4570358592e-01_wp*zs+5.6471810638_wp)*zt   &
-               &   + (1.5393993508_wp*zs-6.9394762624_wp)*zs+1.2750522650e+01_wp
-               !
-            ptmp(ji,jj) = ( zt / z1_T0 + zn / zd ) * ztm
-               !
-         END DO
-      END DO
+      DO_2D( 1, 1, 1, 1 )
+         !
+         zt  = ctmp   (ji,jj) * z1_T0
+         zs  = SQRT( ABS( psal(ji,jj) + zdeltaS ) * r1_S0 )
+         ztm = tmask(ji,jj,1)
+         !
+         zn = ((((-2.1385727895e-01_wp*zt   &
+            &   - 2.7674419971e-01_wp*zs+1.0728094330_wp)*zt   &
+            &   + (2.6366564313_wp*zs+3.3546960647_wp)*zs-7.8012209473_wp)*zt   &
+            &   + ((1.8835586562_wp*zs+7.3949191679_wp)*zs-3.3937395875_wp)*zs-5.6414948432_wp)*zt   &
+            &   + (((3.5737370589_wp*zs-1.5512427389e+01_wp)*zs+2.4625741105e+01_wp)*zs   &
+            &      +1.9912291000e+01_wp)*zs-3.2191146312e+01_wp)*zt   &
+            &   + ((((5.7153204649e-01_wp*zs-3.0943149543_wp)*zs+9.3052495181_wp)*zs   &
+            &      -9.4528934807_wp)*zs+3.1066408996_wp)*zs-4.3504021262e-01_wp
+            !
+         zd = (2.0035003456_wp*zt   &
+            &   -3.4570358592e-01_wp*zs+5.6471810638_wp)*zt   &
+            &   + (1.5393993508_wp*zs-6.9394762624_wp)*zs+1.2750522650e+01_wp
+            !
+         ptmp(ji,jj) = ( zt / z1_T0 + zn / zd ) * ztm
+            !
+      END_2D
       !
       IF( ln_timing )   CALL timing_stop('eos_pt_from_ct')
@@ -1133,16 +1066,14 @@
          !
          z1_S0 = 1._wp / 35.16504_wp
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zs= SQRT( ABS( psal(ji,jj) ) * z1_S0 )           ! square root salinity
-               ptf(ji,jj) = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs &
-                  &          - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp
-            END DO
-         END DO
+         DO_2D( 1, 1, 1, 1 )
+            zs= SQRT( ABS( psal(ji,jj) ) * z1_S0 )           ! square root salinity
+            ptf(ji,jj) = ((((1.46873e-03_wp*zs-9.64972e-03_wp)*zs+2.28348e-02_wp)*zs &
+               &          - 3.12775e-02_wp)*zs+2.07679e-02_wp)*zs-5.87701e-02_wp
+         END_2D
          ptf(:,:) = ptf(:,:) * psal(:,:)
          !
          IF( PRESENT( pdep ) )   ptf(:,:) = ptf(:,:) - 7.53e-4 * pdep(:,:)
          !
-      CASE ( np_eos80, np_leos )                !==  PT,SP (UNESCO formulation)  ==!
+      CASE ( np_eos80 )                !==  PT,SP (UNESCO formulation)  ==!
          !
          ptf(:,:) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(:,:) )   &
@@ -1190,5 +1121,5 @@
          IF( PRESENT( pdep ) )   ptf = ptf - 7.53e-4 * pdep
          !
-      CASE ( np_eos80, np_leos )                !==  PT,SP (UNESCO formulation)  ==!
+      CASE ( np_eos80 )                !==  PT,SP (UNESCO formulation)  ==!
          !
          ptf = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal )   &
@@ -1242,104 +1173,92 @@
       CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
          !
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  !
-                  zh  = gdept(ji,jj,jk,Kmm) * r1_Z0                                ! depth
-                  zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature
-                  zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
-                  ztm = tmask(ji,jj,jk)                                         ! tmask
-                  !
-                  ! potential energy non-linear anomaly
-                  zn2 = (PEN012)*zt   &
-                     &   + PEN102*zs+PEN002
-                     !
-                  zn1 = ((PEN021)*zt   &
-                     &   + PEN111*zs+PEN011)*zt   &
-                     &   + (PEN201*zs+PEN101)*zs+PEN001
-                     !
-                  zn0 = ((((PEN040)*zt   &
-                     &   + PEN130*zs+PEN030)*zt   &
-                     &   + (PEN220*zs+PEN120)*zs+PEN020)*zt   &
-                     &   + ((PEN310*zs+PEN210)*zs+PEN110)*zs+PEN010)*zt   &
-                     &   + (((PEN400*zs+PEN300)*zs+PEN200)*zs+PEN100)*zs+PEN000
-                     !
-                  zn  = ( zn2 * zh + zn1 ) * zh + zn0
-                  !
-                  ppen(ji,jj,jk)  = zn * zh * r1_rho0 * ztm
-                  !
-                  ! alphaPE non-linear anomaly
-                  zn2 = APE002
-                  !
-                  zn1 = (APE011)*zt   &
-                     &   + APE101*zs+APE001
-                     !
-                  zn0 = (((APE030)*zt   &
-                     &   + APE120*zs+APE020)*zt   &
-                     &   + (APE210*zs+APE110)*zs+APE010)*zt   &
-                     &   + ((APE300*zs+APE200)*zs+APE100)*zs+APE000
-                     !
-                  zn  = ( zn2 * zh + zn1 ) * zh + zn0
-                  !                              
-                  pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_rho0 * ztm
-                  !
-                  ! betaPE non-linear anomaly
-                  zn2 = BPE002
-                  !
-                  zn1 = (BPE011)*zt   &
-                     &   + BPE101*zs+BPE001
-                     !
-                  zn0 = (((BPE030)*zt   &
-                     &   + BPE120*zs+BPE020)*zt   &
-                     &   + (BPE210*zs+BPE110)*zs+BPE010)*zt   &
-                     &   + ((BPE300*zs+BPE200)*zs+BPE100)*zs+BPE000
-                     !
-                  zn  = ( zn2 * zh + zn1 ) * zh + zn0
-                  !                              
-                  pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_rho0 * ztm
-                  !
-               END DO
-            END DO
-         END DO
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+            !
+            zh  = gdept(ji,jj,jk,Kmm) * r1_Z0                                ! depth
+            zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature
+            zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
+            ztm = tmask(ji,jj,jk)                                         ! tmask
+            !
+            ! potential energy non-linear anomaly
+            zn2 = (PEN012)*zt   &
+               &   + PEN102*zs+PEN002
+               !
+            zn1 = ((PEN021)*zt   &
+               &   + PEN111*zs+PEN011)*zt   &
+               &   + (PEN201*zs+PEN101)*zs+PEN001
+               !
+            zn0 = ((((PEN040)*zt   &
+               &   + PEN130*zs+PEN030)*zt   &
+               &   + (PEN220*zs+PEN120)*zs+PEN020)*zt   &
+               &   + ((PEN310*zs+PEN210)*zs+PEN110)*zs+PEN010)*zt   &
+               &   + (((PEN400*zs+PEN300)*zs+PEN200)*zs+PEN100)*zs+PEN000
+               !
+            zn  = ( zn2 * zh + zn1 ) * zh + zn0
+            !
+            ppen(ji,jj,jk)  = zn * zh * r1_rho0 * ztm
+            !
+            ! alphaPE non-linear anomaly
+            zn2 = APE002
+            !
+            zn1 = (APE011)*zt   &
+               &   + APE101*zs+APE001
+               !
+            zn0 = (((APE030)*zt   &
+               &   + APE120*zs+APE020)*zt   &
+               &   + (APE210*zs+APE110)*zs+APE010)*zt   &
+               &   + ((APE300*zs+APE200)*zs+APE100)*zs+APE000
+               !
+            zn  = ( zn2 * zh + zn1 ) * zh + zn0
+            !                              
+            pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_rho0 * ztm
+            !
+            ! betaPE non-linear anomaly
+            zn2 = BPE002
+            !
+            zn1 = (BPE011)*zt   &
+               &   + BPE101*zs+BPE001
+               !
+            zn0 = (((BPE030)*zt   &
+               &   + BPE120*zs+BPE020)*zt   &
+               &   + (BPE210*zs+BPE110)*zs+BPE010)*zt   &
+               &   + ((BPE300*zs+BPE200)*zs+BPE100)*zs+BPE000
+               !
+            zn  = ( zn2 * zh + zn1 ) * zh + zn0
+            !                              
+            pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_rho0 * ztm
+            !
+         END_3D
          !
       CASE( np_seos )                !==  Vallis (2006) simplified EOS  ==!
          !
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zt  = pts(ji,jj,jk,jp_tem) - 10._wp  ! temperature anomaly (t-T0)
-                  zs = pts (ji,jj,jk,jp_sal) - 35._wp  ! abs. salinity anomaly (s-S0)
-                  zh  = gdept(ji,jj,jk,Kmm)              ! depth in meters  at t-point
-                  ztm = tmask(ji,jj,jk)                ! tmask
-                  zn  = 0.5_wp * zh * r1_rho0 * ztm
-                  !                                    ! Potential Energy
-                  ppen(ji,jj,jk) = ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zn
-                  !                                    ! alphaPE
-                  pab_pe(ji,jj,jk,jp_tem) = - rn_a0 * rn_mu1 * zn
-                  pab_pe(ji,jj,jk,jp_sal) =   rn_b0 * rn_mu2 * zn
-                  !
-               END DO
-            END DO
-         END DO
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+            zt  = pts(ji,jj,jk,jp_tem) - 10._wp  ! temperature anomaly (t-T0)
+            zs = pts (ji,jj,jk,jp_sal) - 35._wp  ! abs. salinity anomaly (s-S0)
+            zh  = gdept(ji,jj,jk,Kmm)              ! depth in meters  at t-point
+            ztm = tmask(ji,jj,jk)                ! tmask
+            zn  = 0.5_wp * zh * r1_rho0 * ztm
+            !                                    ! Potential Energy
+            ppen(ji,jj,jk) = ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zn
+            !                                    ! alphaPE
+            pab_pe(ji,jj,jk,jp_tem) = - rn_a0 * rn_mu1 * zn
+            pab_pe(ji,jj,jk,jp_sal) =   rn_b0 * rn_mu2 * zn
+            !
+         END_3D
          !
       CASE( np_leos )                !==  linear ISOMIP EOS  ==!
          !
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zt  = pts(ji,jj,jk,jp_tem) - (-1._wp)  ! temperature anomaly (t-T0)
-                  zs = pts (ji,jj,jk,jp_sal) - 34.2_wp   ! abs. salinity anomaly (s-S0)
-                  zh  = gdept(ji,jj,jk,Kmm)                ! depth in meters  at t-point
-                  ztm = tmask(ji,jj,jk)                  ! tmask
-                  zn  = 0.5_wp * zh * r1_rho0 * ztm
-                  !                                    ! Potential Energy
-                  ppen(ji,jj,jk) = 0.
-                  !                                    ! alphaPE
-                  pab_pe(ji,jj,jk,jp_tem) = 0.
-                  pab_pe(ji,jj,jk,jp_sal) = 0.
-                  !
-               END DO
-            END DO
-         END DO
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+            zt  = pts(ji,jj,jk,jp_tem) - (-1._wp)  ! temperature anomaly (t-T0)
+            zs = pts (ji,jj,jk,jp_sal) - 34.2_wp   ! abs. salinity anomaly (s-S0)
+            zh  = gdept(ji,jj,jk,Kmm)                ! depth in meters  at t-point
+            ztm = tmask(ji,jj,jk)                  ! tmask
+            zn  = 0.5_wp * zh * r1_rho0 * ztm
+            !                                    ! Potential Energy
+            ppen(ji,jj,jk) = 0.
+            !                                    ! alphaPE
+            pab_pe(ji,jj,jk,jp_tem) = 0.
+            pab_pe(ji,jj,jk,jp_sal) = 0.
+            !
+         END_3D
          !
       CASE DEFAULT
@@ -1365,14 +1284,11 @@
       INTEGER  ::   ioptio   ! local integer
       !!
-      NAMELIST/nameos/ ln_TEOS10, ln_EOS80, ln_SEOS   , ln_LEOS, &
-         &             rn_a0    , rn_b0   , rn_lambda1, rn_mu1 , &
-         &                                  rn_lambda2, rn_mu2 , rn_nu
-      !!----------------------------------------------------------------------
-      !
-      REWIND( numnam_ref )              ! Namelist nameos in reference namelist : equation of state
+      NAMELIST/nameos/ ln_TEOS10, ln_EOS80, ln_SEOS, ln_LEOS, rn_a0, rn_b0, &
+         &             rn_lambda1, rn_mu1, rn_lambda2, rn_mu2, rn_nu
+      !!----------------------------------------------------------------------
+      !
       READ  ( numnam_ref, nameos, IOSTAT = ios, ERR = 901 )
 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nameos in reference namelist' )
       !
-      REWIND( numnam_cfg )              ! Namelist nameos in configuration namelist : equation of state
       READ  ( numnam_cfg, nameos, IOSTAT = ios, ERR = 902 )
 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nameos in configuration namelist' )
Index: NEMO/branches/2020/r12377_ticket2386/tests/ISOMIP+/MY_SRC/isfcavgam.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/tests/ISOMIP+/MY_SRC/isfcavgam.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/tests/ISOMIP+/MY_SRC/isfcavgam.F90	(revision 13540)
@@ -91,7 +91,7 @@
          pgs(:,:) = rn_gammas0
       CASE ( 'vel' ) ! gamma is proportional to u*
-         CALL gammats_vel      (                   zutbl, zvtbl, rCd0_top, r_ke0_top,               pgt, pgs )
+         CALL gammats_vel      (                   zutbl, zvtbl, rCd0_top, rn_vtide**2,               pgt, pgs )
       CASE ( 'vel_stab' ) ! gamma depends of stability of boundary layer and u*
-         CALL gammats_vel_stab (Kmm, pttbl, pstbl, zutbl, zvtbl, rCd0_top, r_ke0_top, pqoce, pqfwf, pgt, pgs )
+         CALL gammats_vel_stab (Kmm, pttbl, pstbl, zutbl, zvtbl, rCd0_top, rn_vtide**2, pqoce, pqfwf, pgt, pgs )
       CASE DEFAULT
          CALL ctl_stop('STOP','method to compute gamma (cn_gammablk) is unknown (should not see this)')
Index: NEMO/branches/2020/r12377_ticket2386/tests/ISOMIP+/MY_SRC/isfstp.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/tests/ISOMIP+/MY_SRC/isfstp.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/tests/ISOMIP+/MY_SRC/isfstp.F90	(revision 13540)
@@ -250,5 +250,5 @@
       IF ( l_isfoasis .AND. ln_isf ) THEN
          !
-         CALL ctl_stop( ' ln_ctl and ice shelf not tested' )
+         CALL ctl_stop( 'namelist combination ln_cpl and ln_isf not tested' )
          !
          ! NEMO coupled to ATMO model with isf cavity need oasis method for melt computation 
@@ -291,9 +291,7 @@
       !!----------------------------------------------------------------------
       !
-      REWIND( numnam_ref )              ! Namelist namsbc_rnf in reference namelist : Runoffs 
       READ  ( numnam_ref, namisf, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namisf in reference namelist' )
       !
-      REWIND( numnam_cfg )              ! Namelist namsbc_rnf in configuration namelist : Runoffs
       READ  ( numnam_cfg, namisf, IOSTAT = ios, ERR = 902 )
 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namisf in configuration namelist' )
Index: NEMO/branches/2020/r12377_ticket2386/tests/ISOMIP+/MY_SRC/istate.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/tests/ISOMIP+/MY_SRC/istate.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/tests/ISOMIP+/MY_SRC/istate.F90	(revision 13540)
@@ -41,4 +41,6 @@
    PUBLIC   istate_init   ! routine called by step.F90
 
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -75,5 +77,5 @@
       rhd  (:,:,:  ) = 0._wp   ;   rhop (:,:,:  ) = 0._wp      ! set one for all to 0 at level jpk
       rn2b (:,:,:  ) = 0._wp   ;   rn2  (:,:,:  ) = 0._wp      ! set one for all to 0 at levels 1 and jpk
-      ts   (:,:,:,:,Kaa) = 0._wp                               ! set one for all to 0 at level jpk
+      ts  (:,:,:,:,Kaa) = 0._wp                                   ! set one for all to 0 at level jpk
       rab_b(:,:,:,:) = 0._wp   ;   rab_n(:,:,:,:) = 0._wp      ! set one for all to 0 at level jpk
 #if defined key_agrif
@@ -90,5 +92,5 @@
          !                                    ! ---------------
          numror = 0                           ! define numror = 0 -> no restart file to read
-         neuler = 0                           ! Set time-step indicator at nit000 (euler forward)
+         l_1st_euler = .true.                 ! Set time-step indicator at nit000 (euler forward)
          CALL day_init                        ! model calendar (using both namelist and restart infos)
          !                                    ! Initialization of ocean to zero
@@ -103,11 +105,9 @@
                ! Apply minimum wetdepth criterion
                !
-               DO jj = 1,jpj
-                  DO ji = 1,jpi
-                     IF( ht_0(ji,jj) + ssh(ji,jj,Kbb)  < rn_wdmin1 ) THEN
-                        ssh(ji,jj,Kbb) = tmask(ji,jj,1)*( rn_wdmin1 - (ht_0(ji,jj)) )
-                     ENDIF
-                  END DO
-               END DO 
+               DO_2D( 1, 1, 1, 1 )
+                  IF( ht_0(ji,jj) + ssh(ji,jj,Kbb)  < rn_wdmin1 ) THEN
+                     ssh(ji,jj,Kbb) = tmask(ji,jj,1)*( rn_wdmin1 - (ht_0(ji,jj)) )
+                  ENDIF
+               END_2D
             ENDIF 
             uu  (:,:,:,Kbb) = 0._wp
@@ -159,15 +159,11 @@
       !
 !!gm  the use of umsak & vmask is not necessary below as uu(:,:,:,Kmm), vv(:,:,:,Kmm), uu(:,:,:,Kbb), vv(:,:,:,Kbb) are always masked
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               uu_b(ji,jj,Kmm) = uu_b(ji,jj,Kmm) + e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) * umask(ji,jj,jk)
-               vv_b(ji,jj,Kmm) = vv_b(ji,jj,Kmm) + e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm) * vmask(ji,jj,jk)
-               !
-               uu_b(ji,jj,Kbb) = uu_b(ji,jj,Kbb) + e3u(ji,jj,jk,Kbb) * uu(ji,jj,jk,Kbb) * umask(ji,jj,jk)
-               vv_b(ji,jj,Kbb) = vv_b(ji,jj,Kbb) + e3v(ji,jj,jk,Kbb) * vv(ji,jj,jk,Kbb) * vmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+         uu_b(ji,jj,Kmm) = uu_b(ji,jj,Kmm) + e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) * umask(ji,jj,jk)
+         vv_b(ji,jj,Kmm) = vv_b(ji,jj,Kmm) + e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm) * vmask(ji,jj,jk)
+         !
+         uu_b(ji,jj,Kbb) = uu_b(ji,jj,Kbb) + e3u(ji,jj,jk,Kbb) * uu(ji,jj,jk,Kbb) * umask(ji,jj,jk)
+         vv_b(ji,jj,Kbb) = vv_b(ji,jj,Kbb) + e3v(ji,jj,jk,Kbb) * vv(ji,jj,jk,Kbb) * vmask(ji,jj,jk)
+      END_3D
       !
       uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * r1_hu(:,:,Kmm)
Index: NEMO/branches/2020/r12377_ticket2386/tests/ISOMIP+/MY_SRC/sbcfwb.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/tests/ISOMIP+/MY_SRC/sbcfwb.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/tests/ISOMIP+/MY_SRC/sbcfwb.F90	(revision 13540)
@@ -95,4 +95,5 @@
          snwice_mass_b(:,:) = 0.e0               ! no sea-ice model is being used : no snow+ice mass
          snwice_mass  (:,:) = 0.e0
+         snwice_fmass (:,:) = 0.e0
 #endif
          !
@@ -151,5 +152,5 @@
          ENDIF   
          !                                         ! Update fwfold if new year start
-         ikty = 365 * 86400 / rn_Dt               !!bug  use of 365 days leap year or 360d year !!!!!!!
+         ikty = 365 * 86400 / rn_Dt                  !!bug  use of 365 days leap year or 360d year !!!!!!!
          IF( MOD( kt, ikty ) == 0 ) THEN
             a_fwb_b = a_fwb                           ! mean sea level taking into account the ice+snow
@@ -211,5 +212,5 @@
             erp(:,:) = erp(:,:) + zerp_cor(:,:)
             !
-            IF( nprint == 1 .AND. lwp ) THEN                   ! control print
+            IF( lwp ) THEN                   ! control print
                IF( z_fwf < 0._wp ) THEN
                   WRITE(numout,*)'   z_fwf < 0'
Index: NEMO/branches/2020/r12377_ticket2386/tests/ISOMIP+/MY_SRC/tradmp.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/tests/ISOMIP+/MY_SRC/tradmp.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/tests/ISOMIP+/MY_SRC/tradmp.F90	(revision 13540)
@@ -51,4 +51,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   resto    !: restoring coeff. on T and S (s-1)
 
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -110,41 +112,29 @@
       CASE( 0 )                        !*  newtonian damping throughout the water column  *!
          DO jn = 1, jpts
-            DO jk = 1, jpkm1
-               DO jj = 2, jpjm1
-                  DO ji = fs_2, fs_jpim1   ! vector opt.
-                     pts(ji,jj,jk,jn,Krhs) = pts(ji,jj,jk,jn,Krhs)           &
-                        &                  + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jn) - pts(ji,jj,jk,jn,Kbb) )
-                  END DO
-               END DO
-            END DO
+            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+               pts(ji,jj,jk,jn,Krhs) = pts(ji,jj,jk,jn,Krhs)           &
+                  &                  + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jn) - pts(ji,jj,jk,jn,Kbb) )
+            END_3D
          END DO
          !
       CASE ( 1 )                       !*  no damping in the turbocline (avt > 5 cm2/s)  *!
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  IF( avt(ji,jj,jk) <= avt_c ) THEN
-                     pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   &
-                        &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - pts(ji,jj,jk,jp_tem,Kbb) )
-                     pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs)   &
-                        &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - pts(ji,jj,jk,jp_sal,Kbb) )
-                  ENDIF
-               END DO
-            END DO
-         END DO
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+            IF( avt(ji,jj,jk) <= avt_c ) THEN
+               pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   &
+                  &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - pts(ji,jj,jk,jp_tem,Kbb) )
+               pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs)   &
+                  &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - pts(ji,jj,jk,jp_sal,Kbb) )
+            ENDIF
+         END_3D
          !
       CASE ( 2 )                       !*  no damping in the mixed layer   *!
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  IF( gdept(ji,jj,jk,Kmm) >= hmlp (ji,jj) ) THEN
-                     pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   &
-                        &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - pts(ji,jj,jk,jp_tem,Kbb) )
-                     pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs)   &
-                        &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - pts(ji,jj,jk,jp_sal,Kbb) )
-                  ENDIF
-               END DO
-            END DO
-         END DO
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+            IF( gdept(ji,jj,jk,Kmm) >= hmlp (ji,jj) ) THEN
+               pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   &
+                  &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - pts(ji,jj,jk,jp_tem,Kbb) )
+               pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs)   &
+                  &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - pts(ji,jj,jk,jp_sal,Kbb) )
+            ENDIF
+         END_3D
          !
       END SELECT
@@ -157,6 +147,6 @@
       ENDIF
       !                           ! Control print
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' dmp  - Ta: ', mask1=tmask,   &
-         &                       tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
+      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' dmp  - Ta: ', mask1=tmask,   &
+         &                                  tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
       !
       IF( ln_timing )   CALL timing_stop('tra_dmp')
@@ -178,9 +168,7 @@
       !!----------------------------------------------------------------------
       !
-      REWIND( numnam_ref )   ! Namelist namtra_dmp in reference namelist : T & S relaxation
       READ  ( numnam_ref, namtra_dmp, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_dmp in reference namelist' )
       !
-      REWIND( numnam_cfg )   ! Namelist namtra_dmp in configuration namelist : T & S relaxation
       READ  ( numnam_cfg, namtra_dmp, IOSTAT = ios, ERR = 902 )
 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtra_dmp in configuration namelist' )
@@ -220,5 +208,5 @@
          !                          ! Read in mask from file
          CALL iom_open ( cn_resto, imask)
-         CALL iom_get  ( imask, jpdom_autoglo, 'resto', resto )
+         CALL iom_get  ( imask, jpdom_auto, 'resto', resto )
          CALL iom_close( imask )
       ENDIF
