Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/C14/trcatm_c14.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/C14/trcatm_c14.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/C14/trcatm_c14.F90	(revision 13710)
@@ -120,5 +120,5 @@
             IF( ierr3 /= 0 )   CALL ctl_stop( 'STOP', 'trc_atm_c14_ini: unable to allocate fareaz' )
       !
-            DO_2D_11_11
+            DO_2D( 1, 1, 1, 1 )                 ! from C14b package
               IF( gphit(ji,jj) >= yn40 ) THEN
                  fareaz(ji,jj,1) = 0.
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/C14/trcini_c14.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/C14/trcini_c14.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/C14/trcini_c14.F90	(revision 13710)
@@ -69,8 +69,8 @@
         ! 
         CALL iom_get( numrtr, 'co2sbc', co2sbc ) 
-        CALL iom_get( numrtr, jpdom_autoglo, 'c14sbc', c14sbc ) 
-        CALL iom_get( numrtr, jpdom_autoglo, 'exch_co2', exch_co2 ) 
-        CALL iom_get( numrtr, jpdom_autoglo, 'exch_c14', exch_c14 ) 
-        CALL iom_get( numrtr, jpdom_autoglo, 'qtr_c14', qtr_c14 )
+        CALL iom_get( numrtr, jpdom_auto, 'c14sbc', c14sbc ) 
+        CALL iom_get( numrtr, jpdom_auto, 'exch_co2', exch_co2 ) 
+        CALL iom_get( numrtr, jpdom_auto, 'exch_c14', exch_c14 ) 
+        CALL iom_get( numrtr, jpdom_auto, 'qtr_c14', qtr_c14 )
         !
       END IF
@@ -85,5 +85,5 @@
       ELSE
         !
-        CALL iom_get( numrtr, jpdom_autoglo, 'qint_c14', qint_c14 ) 
+        CALL iom_get( numrtr, jpdom_auto, 'qint_c14', qint_c14 ) 
         !
       ENDIF
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/C14/trcsms_c14.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/C14/trcsms_c14.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/C14/trcsms_c14.F90	(revision 13710)
@@ -28,4 +28,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -80,5 +81,5 @@
       ! -------------------------------------------------------------------
 
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          IF( tmask(ji,jj,1) >  0. ) THEN
             !
@@ -127,10 +128,10 @@
       !
       ! Add the surface flux to the trend of jp_c14
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          tr(ji,jj,1,jp_c14,Krhs) = tr(ji,jj,1,jp_c14,Krhs) + qtr_c14(ji,jj) / e3t(ji,jj,1,Kmm) 
       END_2D
       !
       ! Computation of decay effects on jp_c14
-      DO_3D_11_11( 1, jpk )
+      DO_3D( 1, 1, 1, 1, 1, jpk )
          !
          tr(ji,jj,jk,jp_c14,Krhs) = tr(ji,jj,jk,jp_c14,Krhs) - rlam14 * tr(ji,jj,jk,jp_c14,Kbb) * tmask(ji,jj,jk) 
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/C14/trcwri_c14.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/C14/trcwri_c14.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/C14/trcwri_c14.F90	(revision 13710)
@@ -60,5 +60,5 @@
          zz3d(:,:,:) = 0._wp
          !
-         DO_3D_11_11( 1, jpkm1 )
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
             IF( tmask(ji,jj,jk) > 0._wp) THEN
                z3d (ji,jj,jk) = tr(ji,jj,jk,jp_c14,Kmm)
@@ -71,5 +71,5 @@
          z2d(:,:) =0._wp
          jk = 1
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             ztemp = zres(ji,jj) / c14sbc(ji,jj)
             IF( ztemp > 0._wp .AND. tmask(ji,jj,jk) > 0._wp ) z2d(ji,jj) = LOG( ztemp )
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/CFC/trcini_cfc.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/CFC/trcini_cfc.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/CFC/trcini_cfc.F90	(revision 13710)
@@ -132,5 +132,5 @@
       !---------------------------------------------------------------------------------------
       zyd = ylatn - ylats      
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          IF(     gphit(ji,jj) >= ylatn ) THEN   ;   xphem(ji,jj) = 1.e0
          ELSEIF( gphit(ji,jj) <= ylats ) THEN   ;   xphem(ji,jj) = 0.e0
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/CFC/trcsms_cfc.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/CFC/trcsms_cfc.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/CFC/trcsms_cfc.F90	(revision 13710)
@@ -49,4 +49,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -125,6 +126,6 @@
          
          !                                                         !------------!
-         DO_2D_11_11
- 
+         DO_2D( 1, 1, 1, 1 )                                       !  i-j loop  !
+            !                                                      !------------!
             ! space interpolation
             zpp_cfc  =       xphem(ji,jj)   * zpatm(1,jl)   &
@@ -297,5 +298,5 @@
          DO jn = jp_cfc0, jp_cfc1
             jl = jl + 1
-            CALL iom_get( numrtr, jpdom_autoglo, 'qint_'//ctrcnm(jn), qint_cfc(:,:,jl) ) 
+            CALL iom_get( numrtr, jpdom_auto, 'qint_'//ctrcnm(jn), qint_cfc(:,:,jl) ) 
          END DO
       ENDIF
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P2Z/p2zbio.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P2Z/p2zbio.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P2Z/p2zbio.F90	(revision 13710)
@@ -19,5 +19,5 @@
    !
    USE lbclnk          ! 
-   USE prtctl_trc      ! Print control for debbuging
+   USE prtctl          ! Print control for debbuging
    USE iom             !
    
@@ -58,4 +58,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -121,5 +122,5 @@
       DO jk = 1, jpkbm1                      !  Upper ocean (bio-layers)  !
          !                                   ! -------------------------- !
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             ! trophic variables( det, zoo, phy, no3, nh4, dom)
             ! ------------------------------------------------
@@ -241,5 +242,5 @@
       DO jk = jpkb, jpkm1                    !  Upper ocean (bio-layers)  !
          !                                   ! -------------------------- !
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             ! remineralisation of all quantities towards nitrate 
 
@@ -338,6 +339,6 @@
       !
       IF( lk_iomput ) THEN
-         CALL lbc_lnk( 'p2zbio', zw2d(:,:,:),'T', 1. )
-         CALL lbc_lnk_multi( 'p2zbio', zw3d(:,:,:,1),'T', 1., zw3d(:,:,:,2),'T', 1., zw3d(:,:,:,3),'T', 1. )
+         CALL lbc_lnk( 'p2zbio', zw2d(:,:,:),'T', 1.0_wp )
+         CALL lbc_lnk_multi( 'p2zbio', zw3d(:,:,:,1),'T', 1.0_wp, zw3d(:,:,:,2),'T', 1.0_wp, zw3d(:,:,:,3),'T', 1.0_wp )
          ! Save diagnostics
          CALL iom_put( "TNO3PHY", zw2d(:,:,1) )
@@ -366,6 +367,6 @@
       IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
          WRITE(charout, FMT="('bio')")
-         CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+         CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
       ENDIF
       !
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P2Z/p2zexp.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P2Z/p2zexp.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P2Z/p2zexp.F90	(revision 13710)
@@ -17,5 +17,5 @@
    USE p2zsed
    USE lbclnk
-   USE prtctl_trc      ! Print control for debbuging
+   USE prtctl          ! Print control for debbuging
    USE trd_oce
    USE trdtrc
@@ -39,4 +39,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -81,5 +82,5 @@
       ! LAYERS IS DETERMINED BY DMIN3 DEFINED IN sms_p2z.F90
       ! ----------------------------------------------------------------------
-      DO_3D_00_00( 1, jpkm1 )
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
          ze3t = 1. / e3t(ji,jj,jk,Kmm)
          tr(ji,jj,jk,jpno3,Krhs) = tr(ji,jj,jk,jpno3,Krhs) + ze3t * dmin3(ji,jj,jk) * xksi(ji,jj)
@@ -92,5 +93,5 @@
       zgeolpoc = 0.e0         !     Initialization
       ! Release of nutrients from the "simple" sediment
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          ikt = mbkt(ji,jj) 
          tr(ji,jj,ikt,jpno3,Krhs) = tr(ji,jj,ikt,jpno3,Krhs) + sedlam * sedpocn(ji,jj) / e3t(ji,jj,ikt,Kmm) 
@@ -102,9 +103,9 @@
       END_2D
 
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          tr(ji,jj,1,jpno3,Krhs) = tr(ji,jj,1,jpno3,Krhs) + zgeolpoc * cmask(ji,jj) / areacot / e3t(ji,jj,1,Kmm)
       END_2D
 
-      CALL lbc_lnk( 'p2zexp', sedpocn, 'T', 1. )
+      CALL lbc_lnk( 'p2zexp', sedpocn, 'T', 1.0_wp )
  
       ! Oa & Ek: diagnostics depending on jpdia2d !          left as example
@@ -120,5 +121,5 @@
       ELSE
         !
-        DO_2D_11_11
+        DO_2D( 1, 1, 1, 1 )
            zsedpocd = zsedpoca(ji,jj) - 2. * sedpocn(ji,jj) + sedpocb(ji,jj)      ! time laplacian on tracers
            sedpocb(ji,jj) = sedpocn(ji,jj) + rn_atfp * zsedpocd                     ! sedpocb <-- filtered sedpocn
@@ -139,6 +140,6 @@
       IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
          WRITE(charout, FMT="('exp')")
-         CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+         CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
       ENDIF
       !
@@ -173,5 +174,5 @@
       zdm0 = 0._wp
       zrro = 1._wp
-      DO_3D_11_11( jpkb, jpkm1 )
+      DO_3D( 1, 1, 1, 1, jpkb, jpkm1 )
          zfluo = ( gdepw(ji,jj,jk  ,Kmm) / gdepw(ji,jj,jpkb,Kmm) )**xhr
          zfluu = ( gdepw(ji,jj,jk+1,Kmm) / gdepw(ji,jj,jpkb,Kmm) )**xhr
@@ -190,5 +191,5 @@
       dminl(:,:)   = 0._wp
       dmin3(:,:,:) = zdm0
-      DO_3D_11_11( 1, jpk )
+      DO_3D( 1, 1, 1, 1, 1, jpk )
          IF( tmask(ji,jj,jk) == 0._wp ) THEN
             dminl(ji,jj) = dminl(ji,jj) + dmin3(ji,jj,jk)
@@ -197,5 +198,5 @@
       END_3D
 
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          IF( tmask(ji,jj,1) == 0 )   dmin3(ji,jj,1) = 0._wp
       END_2D
@@ -203,5 +204,5 @@
       ! Coastal mask 
       cmask(:,:) = 0._wp
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          IF( tmask(ji,jj,1) /= 0. ) THEN
             zmaskt = tmask(ji+1,jj,1) * tmask(ji-1,jj,1) * tmask(ji,jj+1,1) * tmask(ji,jj-1,1) 
@@ -209,10 +210,10 @@
          END IF
       END_2D
-      CALL lbc_lnk( 'p2zexp', cmask , 'T', 1. )      ! lateral boundary conditions on cmask   (sign unchanged)
+      CALL lbc_lnk( 'p2zexp', cmask , 'T', 1.0_wp )      ! lateral boundary conditions on cmask   (sign unchanged)
       areacot = glob_sum( 'p2zexp', e1e2t(:,:) * cmask(:,:) )
       !
       IF( ln_rsttr ) THEN
-         CALL iom_get( numrtr, jpdom_autoglo, 'SEDB'//ctrcnm(jpdet), sedpocb(:,:) )
-         CALL iom_get( numrtr, jpdom_autoglo, 'SEDN'//ctrcnm(jpdet), sedpocn(:,:) )
+         CALL iom_get( numrtr, jpdom_auto, 'SEDB'//ctrcnm(jpdet), sedpocb(:,:) )
+         CALL iom_get( numrtr, jpdom_auto, 'SEDN'//ctrcnm(jpdet), sedpocn(:,:) )
       ELSE
          sedpocb(:,:) = 0._wp
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P2Z/p2zopt.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P2Z/p2zopt.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P2Z/p2zopt.F90	(revision 13710)
@@ -18,5 +18,5 @@
    USE trc
    USE sms_pisces
-   USE prtctl_trc      ! Print control for debbuging
+   USE prtctl          ! Print control for debbuging
 
    IMPLICIT NONE
@@ -40,4 +40,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -94,5 +95,5 @@
       !                                          ! Photosynthetically Available Radiation (PAR)
       zcoef = 12 * redf / rcchl / rpig           ! --------------------------------------
-      DO_3D_11_11( 2, jpk )
+      DO_3D( 1, 1, 1, 1, 2, jpk )                     ! local par at w-levels
          zpig = LOG(  MAX( TINY(0.), tr(ji,jj,jk-1,jpphy,Kmm) ) * zcoef  )
          zkr  = xkr0 + xkrp * EXP( xlr * zpig )
@@ -101,5 +102,5 @@
          zparg(ji,jj,jk) = zparg(ji,jj,jk-1) * EXP( -zkg * e3t(ji,jj,jk-1,Kmm) )
       END_3D
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )                   ! mean par at t-levels
          zpig = LOG(  MAX( TINY(0.), tr(ji,jj,jk,jpphy,Kmm) ) * zcoef  )
          zkr  = xkr0 + xkrp * EXP( xlr * zpig )
@@ -113,9 +114,9 @@
       !                                          ! --------------
       neln(:,:) = 1                                   ! euphotic layer level
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )                   ! (i.e. 1rst T-level strictly below EL bottom)
         IF( etot(ji,jj,jk) >= zpar100(ji,jj) )   neln(ji,jj) = jk + 1 
       END_3D
       !                                               ! Euphotic layer depth
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          heup(ji,jj) = gdepw(ji,jj,neln(ji,jj),Kmm)
       END_2D
@@ -124,6 +125,6 @@
       IF(sn_cfctl%l_prttrc) THEN      ! print mean trends (used for debugging)
          WRITE(charout, FMT="('opt')")
-         CALL prt_ctl_trc_info( charout )
-         CALL prt_ctl_trc( tab4d=tr(:,:,:,:,Kmm), mask=tmask, clinfo=ctrcnm )
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+         CALL prt_ctl( tab4d_1=tr(:,:,:,:,Kmm), mask1=tmask, clinfo=ctrcnm )
       ENDIF
       !
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P2Z/p2zsed.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P2Z/p2zsed.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P2Z/p2zsed.F90	(revision 13710)
@@ -18,5 +18,5 @@
    USE lbclnk          !
    USE iom             !
-   USE prtctl_trc      ! Print control for debbuging
+   USE prtctl          ! Print control for debbuging
 
    IMPLICIT NONE
@@ -33,4 +33,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -88,5 +89,5 @@
 
       ! tracer flux divergence at t-point added to the general trend
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          ztra(ji,jj,jk)  = - ( zwork(ji,jj,jk) - zwork(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm)
          tr(ji,jj,jk,jpdet,Krhs) = tr(ji,jj,jk,jpdet,Krhs) + ztra(ji,jj,jk) 
@@ -108,6 +109,6 @@
       IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
          WRITE(charout, FMT="('sed')")
-         CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+         CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
       ENDIF
       !
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zagg.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zagg.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zagg.F90	(revision 13710)
@@ -17,5 +17,5 @@
    USE trc             !  passive tracers common variables 
    USE sms_pisces      !  PISCES Source Minus Sink variables
-   USE prtctl_trc      !  print control for debugging
+   USE prtctl          !  print control for debugging
 
    IMPLICIT NONE
@@ -60,5 +60,5 @@
       IF( ln_p4z ) THEN
          !
-         DO_3D_11_11( 1, jpkm1 )
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
             !
             zfact = xstep * xdiss(ji,jj,jk)
@@ -102,5 +102,5 @@
       ELSE    ! ln_p5z
         !
-         DO_3D_11_11( 1, jpkm1 )
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
             !
             zfact = xstep * xdiss(ji,jj,jk)
@@ -170,6 +170,6 @@
       IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
          WRITE(charout, FMT="('agg')")
-         CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+         CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
       ENDIF
       !
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zbc.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zbc.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zbc.F90	(revision 13710)
@@ -48,4 +48,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -111,5 +112,5 @@
       IF( ll_river ) THEN
           jl = n_trc_indcbc(jpno3)
-          DO_2D_11_11
+          DO_2D( 1, 1, 1, 1 )
              DO jk = 1, nk_rnf(ji,jj)
                 zcoef = rn_rfact / ( e1e2t(ji,jj) * h_rnf(ji,jj) * rn_cbc_time ) * tmask(ji,jj,1)
@@ -144,5 +145,5 @@
          ALLOCATE( zironice(jpi,jpj) )
          !
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             zdep    = rfact / e3t(ji,jj,1,Kmm)
             zwflux  = fmmflx(ji,jj) / 1000._wp
@@ -287,5 +288,5 @@
          CALL iom_open ( TRIM( sn_ironsed%clname ), numiron )
          ALLOCATE( zcmask(jpi,jpj,jpk) )
-         CALL iom_get  ( numiron, jpdom_data, TRIM( sn_ironsed%clvar ), zcmask(:,:,:), 1 )
+         CALL iom_get  ( numiron, jpdom_global, TRIM( sn_ironsed%clvar ), zcmask(:,:,:), 1 )
          CALL iom_close( numiron )
          !
@@ -296,5 +297,5 @@
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) ' Level corresponding to 50m depth ',  ik50,' ', gdept_1d(ik50+1)
-         DO_3D_00_00( 1, ik50 )
+         DO_3D( 0, 0, 0, 0, 1, ik50 )
             ze3t   = e3t_0(ji,jj,jk)
             zsurfc =  e1u(ji,jj) * ( 1. - umask(ji  ,jj  ,jk) )   &
@@ -310,7 +311,7 @@
          END_3D
          !
-         CALL lbc_lnk( 'p4zbc', zcmask , 'T', 1. )      ! lateral boundary conditions on cmask   (sign unchanged)
-         !
-         DO_3D_11_11( 1, jpk )
+         CALL lbc_lnk( 'p4zbc', zcmask , 'T', 1.0_wp )      ! lateral boundary conditions on cmask   (sign unchanged)
+         !
+         DO_3D( 1, 1, 1, 1, 1, jpk )
             zexpide   = MIN( 8.,( gdept(ji,jj,jk,Kmm) / 500. )**(-1.5) )
             zdenitide = -0.9543 + 0.7662 * LOG( zexpide ) - 0.235 * LOG( zexpide )**2
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zbio.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zbio.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zbio.F90	(revision 13710)
@@ -30,5 +30,5 @@
    USE p4zfechem
    USE p4zligand       !  Prognostic ligand model
-   USE prtctl_trc      !  print control for debugging
+   USE prtctl          !  print control for debugging
    USE iom             !  I/O manager
   
@@ -40,4 +40,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -71,5 +72,5 @@
       xdiss(:,:,:) = 1.
 !!gm the use of nmld should be better here?
-      DO_3D_11_11( 2, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 2, jpkm1 )
 !!gm  :  use nmln  and test on jk ...  less memory acces
          IF( gdepw(ji,jj,jk+1,Kmm) > hmld(ji,jj) )   xdiss(ji,jj,jk) = 0.01
@@ -107,6 +108,6 @@
       IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
          WRITE(charout, FMT="('bio ')")
-         CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+         CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
       ENDIF
       !
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zche.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zche.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zche.F90	(revision 13710)
@@ -132,4 +132,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -178,5 +179,5 @@
       ! 0.04°C relative to an exact computation
       ! ---------------------------------------------------------------------
-      DO_3D_11_11( 1, jpk )
+      DO_3D( 1, 1, 1, 1, 1, jpk )
          zpres = gdept(ji,jj,jk,Kmm) / 1000.
          za1 = 0.04 * ( 1.0 + 0.185 * ts(ji,jj,jk,jp_tem,Kmm) + 0.035 * (salinprac(ji,jj,jk) - 35.0) )
@@ -471,5 +472,5 @@
       IF( ln_timing )  CALL timing_start('ahini_for_at')
       !
-      DO_3D_11_11( 1, jpk )
+      DO_3D( 1, 1, 1, 1, 1, jpk )
       p_alkcb  = tr(ji,jj,jk,jptal,Kbb) * 1000. / (rhop(ji,jj,jk) + rtrn)
       p_dictot = tr(ji,jj,jk,jpdic,Kbb) * 1000. / (rhop(ji,jj,jk) + rtrn)
@@ -569,5 +570,5 @@
 
    ! TOTAL H+ scale: conversion factor for Htot = aphscale * Hfree
-   DO_3D_11_11( 1, jpk )
+   DO_3D( 1, 1, 1, 1, 1, jpk )
       IF (rmask(ji,jj,jk) == 1.) THEN
          p_alktot = tr(ji,jj,jk,jptal,Kbb) * 1000. / (rhop(ji,jj,jk) + rtrn)
@@ -598,5 +599,5 @@
 
    DO jn = 1, jp_maxniter_atgen 
-   DO_3D_11_11( 1, jpk )
+   DO_3D( 1, 1, 1, 1, 1, jpk )
       IF (rmask(ji,jj,jk) == 1.) THEN
          zfact = rhop(ji,jj,jk) / 1000. + rtrn
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zfechem.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zfechem.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zfechem.F90	(revision 13710)
@@ -16,5 +16,5 @@
    USE p4zche          ! chemical model
    USE p4zbc           ! Boundary conditions from sediments
-   USE prtctl_trc      ! print control for debugging
+   USE prtctl          ! print control for debugging
    USE iom             ! I/O manager
 
@@ -33,4 +33,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -91,5 +92,5 @@
       ! Chemistry is supposed to be fast enough to be at equilibrium
       ! ------------------------------------------------------------
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          zTL1(ji,jj,jk)  = ztotlig(ji,jj,jk)
          zkeq            = fekeq(ji,jj,jk)
@@ -106,5 +107,5 @@
 
       zdust = 0.         ! if no dust available
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          ! Scavenging rate of iron. This scavenging rate depends on the load of particles of sea water. 
          ! This parameterization assumes a simple second order kinetics (k[Particles][Fe]).
@@ -117,8 +118,4 @@
          !
          zfeequi = zFe3(ji,jj,jk) * 1E-9
-         zhplus  = max( rtrn, hi(ji,jj,jk) )
-         fe3sol  = fesol(ji,jj,jk,1) * ( zhplus**3 + fesol(ji,jj,jk,2) * zhplus**2  &
-            &         + fesol(ji,jj,jk,3) * zhplus + fesol(ji,jj,jk,4)     &
-            &         + fesol(ji,jj,jk,5) / zhplus )
          zfecoll = 0.5 * zFeL1(ji,jj,jk) * 1E-9
          ! precipitation of Fe3+, creation of nanoparticles
@@ -176,5 +173,5 @@
       IF( ln_ligand ) THEN
          !
-         DO_3D_11_11( 1, jpkm1 )
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
             zlam1a   = ( 0.369  * 0.3 * tr(ji,jj,jk,jpdoc,Kbb) + 102.4  * tr(ji,jj,jk,jppoc,Kbb) ) * xdiss(ji,jj,jk)    &
                 &    + ( 114.   * 0.3 * tr(ji,jj,jk,jpdoc,Kbb) )
@@ -221,6 +218,6 @@
       IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
          WRITE(charout, FMT="('fechem')")
-         CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+         CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
       ENDIF
       !
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zflx.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zflx.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zflx.F90	(revision 13710)
@@ -19,5 +19,5 @@
    USE sms_pisces     !  PISCES Source Minus Sink variables
    USE p4zche         !  Chemical model
-   USE prtctl_trc     !  print control for debugging
+   USE prtctl         !  print control for debugging
    USE iom            !  I/O manager
    USE fldread        !  read input fields
@@ -54,4 +54,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -109,5 +110,5 @@
       IF( l_co2cpl )   satmco2(:,:) = atm_co2(:,:)
 
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          ! DUMMY VARIABLES FOR DIC, H+, AND BORATE
          zfact = rhop(ji,jj,1) / 1000. + rtrn
@@ -125,5 +126,5 @@
       ! -------------------------------------------
 
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          ztc  = MIN( 35., ts(ji,jj,1,jp_tem,Kmm) )
          ztc2 = ztc * ztc
@@ -144,5 +145,5 @@
 
 
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          ztkel = tempis(ji,jj,1) + 273.15
          zsal  = salinprac(ji,jj,1) + ( 1.- tmask(ji,jj,1) ) * 35.
@@ -177,6 +178,6 @@
       IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
          WRITE(charout, FMT="('flx ')")
-         CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+         CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
       ENDIF
 
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zligand.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zligand.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zligand.F90	(revision 13710)
@@ -12,5 +12,5 @@
    USE trc             ! passive tracers common variables 
    USE sms_pisces      ! PISCES Source Minus Sink variables
-   USE prtctl_trc      ! print control for debugging
+   USE prtctl          ! print control for debugging
    USE iom             !  I/O manager
 
@@ -52,5 +52,5 @@
       IF( ln_timing )   CALL timing_start('p4z_ligand')
       !
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          !
          ! ------------------------------------------------------------------
@@ -89,6 +89,6 @@
       IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
          WRITE(charout, FMT="('ligand1')")
-         CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+         CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
       ENDIF
       !
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zlim.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zlim.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zlim.F90	(revision 13710)
@@ -98,5 +98,5 @@
       IF( ln_timing )   CALL timing_start('p4z_lim')
       !
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          
          ! Tuning of the iron concentration to a minimum level that is set to the detection limit
@@ -161,5 +161,5 @@
          zlim1    = xdiatno3(ji,jj,jk) + xdiatnh4(ji,jj,jk)
          zlim2    = tr(ji,jj,jk,jppo4,Kbb) / ( tr(ji,jj,jk,jppo4,Kbb) + zconc1dnh4  )
-         zlim3    = tr(ji,jj,jk,jpsil,Kbb) / ( tr(ji,jj,jk,jpsil,Kbb) + xksi(ji,jj) )
+         zlim3    = tr(ji,jj,jk,jpsil,Kbb) / ( tr(ji,jj,jk,jpsil,Kbb) + xksi(ji,jj) + rtrn )
          zratio   = tr(ji,jj,jk,jpdfe,Kbb) * z1_trbdia
          zironmin = xcoef1 * tr(ji,jj,jk,jpdch,Kbb) * z1_trbdia + xcoef2 * zlim1 + xcoef3 * xdiatno3(ji,jj,jk)
@@ -173,5 +173,5 @@
       ! Compute the fraction of nanophytoplankton that is made of calcifiers
       ! --------------------------------------------------------------------
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          zlim1 =  ( tr(ji,jj,jk,jpno3,Kbb) * concnnh4 + tr(ji,jj,jk,jpnh4,Kbb) * concnno3 )    &
             &   / ( concnno3 * concnnh4 + concnnh4 * tr(ji,jj,jk,jpno3,Kbb) + concnno3 * tr(ji,jj,jk,jpnh4,Kbb) ) 
@@ -193,5 +193,5 @@
       END_3D
       !
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          ! denitrification factor computed from O2 levels
          nitrfac(ji,jj,jk) = MAX(  0.e0, 0.4 * ( 6.e-6  - tr(ji,jj,jk,jpoxy,Kbb) )    &
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zlys.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zlys.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zlys.F90	(revision 13710)
@@ -20,5 +20,5 @@
    USE sms_pisces      !  PISCES Source Minus Sink variables
    USE p4zche          !  Chemical model
-   USE prtctl_trc      !  print control for debugging
+   USE prtctl          !  print control for debugging
    USE iom             !  I/O manager
 
@@ -75,5 +75,5 @@
       CALL solve_at_general( zhinit, zhi, Kbb )
 
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          zco3(ji,jj,jk) = tr(ji,jj,jk,jpdic,Kbb) * ak13(ji,jj,jk) * ak23(ji,jj,jk) / (zhi(ji,jj,jk)**2   &
             &             + ak13(ji,jj,jk) * zhi(ji,jj,jk) + ak13(ji,jj,jk) * ak23(ji,jj,jk) + rtrn )
@@ -87,5 +87,5 @@
       !     ---------------------------------------------------------
 
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
 
          ! DEVIATION OF [CO3--] FROM SATURATION VALUE
@@ -130,6 +130,6 @@
       IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
         WRITE(charout, FMT="('lys ')")
-        CALL prt_ctl_trc_info(charout)
-        CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
+        CALL prt_ctl_info( charout, cdcomp = 'top' )
+        CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
       ENDIF
       !
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zmeso.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zmeso.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zmeso.F90	(revision 13710)
@@ -15,5 +15,5 @@
    USE sms_pisces      ! PISCES Source Minus Sink variables
    USE p4zprod         ! production
-   USE prtctl_trc      ! print control for debugging
+   USE prtctl          ! print control for debugging
    USE iom             ! I/O manager
 
@@ -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
@@ -80,5 +81,5 @@
       IF( ln_timing )   CALL timing_start('p4z_meso')
       !
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          zcompam   = MAX( ( tr(ji,jj,jk,jpmes,Kbb) - 1.e-9 ), 0.e0 )
          zfact     = xstep * tgfunc2(ji,jj,jk) * zcompam
@@ -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
@@ -237,6 +246,6 @@
       IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
         WRITE(charout, FMT="('meso')")
-        CALL prt_ctl_trc_info(charout)
-        CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
+        CALL prt_ctl_info( charout, cdcomp = 'top' )
+        CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
       ENDIF
       !
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zmicro.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zmicro.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zmicro.F90	(revision 13710)
@@ -17,5 +17,5 @@
    USE p4zprod         ! production
    USE iom             ! I/O manager
-   USE prtctl_trc      ! print control for debugging
+   USE prtctl          ! print control for debugging
 
    IMPLICIT NONE
@@ -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
@@ -78,5 +79,5 @@
       IF( ln_timing )   CALL timing_start('p4z_micro')
       !
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          zcompaz = MAX( ( tr(ji,jj,jk,jpzoo,Kbb) - 1.e-9 ), 0.e0 )
          zfact   = xstep * tgfunc2(ji,jj,jk) * zcompaz
@@ -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 ) 
@@ -193,6 +202,6 @@
       IF(sn_cfctl%l_prttrc) THEN      ! print mean trends (used for debugging)
          WRITE(charout, FMT="('micro')")
-         CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+         CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
       ENDIF
       !
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zmort.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zmort.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zmort.F90	(revision 13710)
@@ -15,5 +15,5 @@
    USE p4zprod         ! Primary productivity 
    USE p4zlim          ! Phytoplankton limitation terms
-   USE prtctl_trc      ! print control for debugging
+   USE prtctl          ! print control for debugging
 
    IMPLICIT NONE
@@ -77,5 +77,5 @@
       !
       prodcal(:,:,:) = 0._wp   ! calcite production variable set to zero
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          zcompaph = MAX( ( tr(ji,jj,jk,jpphy,Kbb) - 1e-8 ), 0.e0 )
          !     When highly limited by macronutrients, very small cells 
@@ -120,6 +120,6 @@
        IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
          WRITE(charout, FMT="('nano')")
-         CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+         CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
        ENDIF
       !
@@ -152,5 +152,5 @@
       !     ------------------------------------------------------------
 
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
 
          zcompadi = MAX( ( tr(ji,jj,jk,jpdia,Kbb) - 1e-9), 0. )
@@ -192,6 +192,6 @@
       IF(sn_cfctl%l_prttrc) THEN      ! print mean trends (used for debugging)
          WRITE(charout, FMT="('diat')")
-         CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+         CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
       ENDIF
       !
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zopt.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zopt.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zopt.F90	(revision 13710)
@@ -16,5 +16,5 @@
    USE iom            ! I/O manager
    USE fldread        !  time interpolation
-   USE prtctl_trc     !  print control for debugging
+   USE prtctl         !  print control for debugging
 
    IMPLICIT NONE
@@ -37,11 +37,8 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   par_varsw      ! PAR fraction of shortwave
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ekb, ekg, ekr  ! wavelength (Red-Green-Blue)
-
-   INTEGER  ::   nksrp   ! levels below which the light cannot penetrate ( depth larger than 391 m)
-
-   REAL(wp), DIMENSION(3,61) ::   xkrgb   ! tabulated attenuation coefficients for RGB absorption
    
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -88,12 +85,12 @@
       IF( ln_p5z )   zchl3d(:,:,:) = zchl3d(:,:,:)    + tr(:,:,:,jppch,Kbb)
       !
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          zchl = ( zchl3d(ji,jj,jk) + rtrn ) * 1.e6
          zchl = MIN(  10. , MAX( 0.05, zchl )  )
          irgb = NINT( 41 + 20.* LOG10( zchl ) + rtrn )
          !                                                         
-         ekb(ji,jj,jk) = xkrgb(1,irgb) * e3t(ji,jj,jk,Kmm)
-         ekg(ji,jj,jk) = xkrgb(2,irgb) * e3t(ji,jj,jk,Kmm)
-         ekr(ji,jj,jk) = xkrgb(3,irgb) * e3t(ji,jj,jk,Kmm)
+         ekb(ji,jj,jk) = rkrgb(1,irgb) * e3t(ji,jj,jk,Kmm)
+         ekg(ji,jj,jk) = rkrgb(2,irgb) * e3t(ji,jj,jk,Kmm)
+         ekr(ji,jj,jk) = rkrgb(3,irgb) * e3t(ji,jj,jk,Kmm)
       END_3D
       !                                        !* Photosynthetically Available Radiation (PAR)
@@ -105,5 +102,5 @@
          CALL p4z_opt_par( kt, Kmm, zqsr_corr, ze1, ze2, ze3, pqsr100 = zqsr100 ) 
          !
-         DO jk = 1, nksrp      
+         DO jk = 1, nksr      
             etot_ndcy(:,:,jk) =        ze1(:,:,jk) +        ze2(:,:,jk) +       ze3(:,:,jk)
             enano    (:,:,jk) =  1.85 * ze1(:,:,jk) + 0.69 * ze2(:,:,jk) + 0.46 * ze3(:,:,jk)
@@ -111,5 +108,5 @@
          END DO
          IF( ln_p5z ) THEN
-            DO jk = 1, nksrp      
+            DO jk = 1, nksr      
               epico  (:,:,jk) =  1.94 * ze1(:,:,jk) + 0.66 * ze2(:,:,jk) + 0.4 * ze3(:,:,jk)
             END DO
@@ -120,5 +117,5 @@
          CALL p4z_opt_par( kt, Kmm, zqsr_corr, ze1, ze2, ze3 ) 
          !
-         DO jk = 1, nksrp      
+         DO jk = 1, nksr      
             etot(:,:,jk) =  ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk)
          END DO
@@ -130,5 +127,5 @@
          CALL p4z_opt_par( kt, Kmm, zqsr_corr, ze1, ze2, ze3, pqsr100 = zqsr100  ) 
          !
-         DO jk = 1, nksrp      
+         DO jk = 1, nksr      
             etot (:,:,jk) =        ze1(:,:,jk) +        ze2(:,:,jk) +       ze3(:,:,jk)
             enano(:,:,jk) =  1.85 * ze1(:,:,jk) + 0.69 * ze2(:,:,jk) + 0.46 * ze3(:,:,jk)
@@ -136,5 +133,5 @@
          END DO
          IF( ln_p5z ) THEN
-            DO jk = 1, nksrp      
+            DO jk = 1, nksr      
               epico(:,:,jk) =  1.94 * ze1(:,:,jk) + 0.66 * ze2(:,:,jk) + 0.4 * ze3(:,:,jk)
             END DO
@@ -149,5 +146,5 @@
          !
          etot3(:,:,1) =  qsr(:,:) * tmask(:,:,1)
-         DO jk = 2, nksrp + 1
+         DO jk = 2, nksr + 1
             etot3(:,:,jk) =  ( ze0(:,:,jk) + ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk) ) * tmask(:,:,jk)
          END DO
@@ -159,5 +156,5 @@
       heup_01(:,:) = gdepw(:,:,2,Kmm)
 
-      DO_3D_11_11( 2, nksrp )
+      DO_3D( 1, 1, 1, 1, 2, nksr )
         IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >=  zqsr100(ji,jj) )  THEN
            neln(ji,jj) = jk+1                    ! Euphotic level : 1rst T-level strictly below Euphotic layer
@@ -177,5 +174,5 @@
       zetmp2 (:,:)   = 0.e0
 
-      DO_3D_11_11( 1, nksrp )
+      DO_3D( 1, 1, 1, 1, 1, nksr )
          IF( gdepw(ji,jj,jk+1,Kmm) <= hmld(ji,jj) ) THEN
             zetmp1 (ji,jj) = zetmp1 (ji,jj) + etot     (ji,jj,jk) * e3t(ji,jj,jk,Kmm) ! remineralisation
@@ -188,5 +185,5 @@
       zpar(:,:,:) = etot_ndcy(:,:,:)  ! diagnostic : PAR with no diurnal cycle 
       !
-      DO_3D_11_11( 1, nksrp )
+      DO_3D( 1, 1, 1, 1, 1, nksr )
          IF( gdepw(ji,jj,jk+1,Kmm) <= hmld(ji,jj) ) THEN
             z1_dep = 1. / ( zdepmoy(ji,jj) + rtrn )
@@ -200,5 +197,5 @@
       zetmp4 (:,:)   = 0.e0
       !
-      DO_3D_11_11( 1, nksrp )
+      DO_3D( 1, 1, 1, 1, 1, nksr )
          IF( gdepw(ji,jj,jk+1,Kmm) <= MIN(hmld(ji,jj), heup_01(ji,jj)) ) THEN
             zetmp3 (ji,jj) = zetmp3 (ji,jj) + enano    (ji,jj,jk) * e3t(ji,jj,jk,Kmm) ! production
@@ -210,5 +207,5 @@
       ediatm(:,:,:) = ediat(:,:,:)
       !
-      DO_3D_11_11( 1, nksrp )
+      DO_3D( 1, 1, 1, 1, 1, nksr )
          IF( gdepw(ji,jj,jk+1,Kmm) <= hmld(ji,jj) ) THEN
             z1_dep = 1. / ( zdepmoy(ji,jj) + rtrn )
@@ -220,5 +217,5 @@
       IF( ln_p5z ) THEN
          ALLOCATE( zetmp5(jpi,jpj) )  ;   zetmp5 (:,:) = 0.e0
-         DO_3D_11_11( 1, nksrp )
+         DO_3D( 1, 1, 1, 1, 1, nksr )
             IF( gdepw(ji,jj,jk+1,Kmm) <= MIN(hmld(ji,jj), heup_01(ji,jj)) ) THEN
                zetmp5(ji,jj)  = zetmp5 (ji,jj) + epico(ji,jj,jk) * e3t(ji,jj,jk,Kmm) ! production
@@ -228,5 +225,5 @@
          epicom(:,:,:) = epico(:,:,:)
          !
-         DO_3D_11_11( 1, nksrp )
+         DO_3D( 1, 1, 1, 1, 1, nksr )
             IF( gdepw(ji,jj,jk+1,Kmm) <= hmld(ji,jj) ) THEN
                z1_dep = 1. / ( zdepmoy(ji,jj) + rtrn )
@@ -282,5 +279,5 @@
          pe3(:,:,1) = zqsr(:,:)
          !
-         DO jk = 2, nksrp + 1
+         DO jk = 2, nksr + 1
             DO jj = 1, jpj
                DO ji = 1, jpi
@@ -301,5 +298,5 @@
         pe3(:,:,1) = zqsr(:,:) * EXP( -0.5 * ekr(:,:,1) )
         !
-        DO_3D_11_11( 2, nksrp )
+        DO_3D( 1, 1, 1, 1, 2, nksr )
            pe1(ji,jj,jk) = pe1(ji,jj,jk-1) * EXP( -0.5 * ( ekb(ji,jj,jk-1) + ekb(ji,jj,jk) ) )
            pe2(ji,jj,jk) = pe2(ji,jj,jk-1) * EXP( -0.5 * ( ekg(ji,jj,jk-1) + ekg(ji,jj,jk) ) )
@@ -399,9 +396,4 @@
          ntimes_par = iom_getszuld( numpar )   ! get number of record in file
       ENDIF
-      !
-      CALL trc_oce_rgb( xkrgb )                  ! tabulated attenuation coefficients
-      nksrp = trc_oce_ext_lev( r_si2, 0.33e2 )     ! max level of light extinction (Blue Chl=0.01)
-      !
-      IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksrp, ' ref depth = ', gdepw_1d(nksrp+1), ' m'
       !
                          ekr      (:,:,:) = 0._wp
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zpoc.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zpoc.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zpoc.F90	(revision 13710)
@@ -15,5 +15,5 @@
    USE trc             !  passive tracers common variables 
    USE sms_pisces      !  PISCES Source Minus Sink variables
-   USE prtctl_trc      !  print control for debugging
+   USE prtctl          !  print control for debugging
    USE iom             !  I/O manager
 
@@ -39,4 +39,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -106,5 +107,5 @@
      ! -----------------------------------------------------------------------
      ztremint(:,:,:) = zremigoc(:,:,:)
-     DO_3D_11_11( 2, jpkm1 )
+     DO_3D( 1, 1, 1, 1, 2, jpkm1 )
         IF (tmask(ji,jj,jk) == 1.) THEN
           zdep = hmld(ji,jj)
@@ -191,5 +192,5 @@
 
       IF( ln_p4z ) THEN
-         DO_3D_11_11( 1, jpkm1 )
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
             ! POC disaggregation by turbulence and bacterial activity. 
             ! --------------------------------------------------------
@@ -211,5 +212,5 @@
          END_3D
       ELSE
-         DO_3D_11_11( 1, jpkm1 )
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
              ! POC disaggregation by turbulence and bacterial activity. 
             ! --------------------------------------------------------
@@ -241,6 +242,6 @@
      IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
         WRITE(charout, FMT="('poc1')")
-        CALL prt_ctl_trc_info(charout)
-        CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
+        CALL prt_ctl_info( charout, cdcomp = 'top' )
+        CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
      ENDIF
 
@@ -259,5 +260,5 @@
      ! ----------------------------------------------------------------
      ! 
-     DO_3D_11_11( 1, jpkm1 )
+     DO_3D( 1, 1, 1, 1, 1, jpkm1 )
         zdep = hmld(ji,jj)
         IF (tmask(ji,jj,jk) == 1. .AND. gdept(ji,jj,jk,Kmm) <= zdep ) THEN
@@ -274,5 +275,5 @@
      ! ---------------------------------------------------------------------
      ztremint(:,:,:) = zremipoc(:,:,:)
-     DO_3D_11_11( 1, jpkm1 )
+     DO_3D( 1, 1, 1, 1, 1, jpkm1 )
         IF (tmask(ji,jj,jk) == 1.) THEN
           zdep = hmld(ji,jj)
@@ -309,5 +310,5 @@
      ! -----------------------------------------------------------------------
      !
-     DO_3D_11_11( 2, jpkm1 )
+     DO_3D( 1, 1, 1, 1, 2, jpkm1 )
         IF (tmask(ji,jj,jk) == 1.) THEN
           zdep = hmld(ji,jj)
@@ -383,5 +384,5 @@
 
      IF( ln_p4z ) THEN
-         DO_3D_11_11( 1, jpkm1 )
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
             IF (tmask(ji,jj,jk) == 1.) THEN
               ! POC disaggregation by turbulence and bacterial activity. 
@@ -400,5 +401,5 @@
          END_3D
      ELSE
-       DO_3D_11_11( 1, jpkm1 )
+       DO_3D( 1, 1, 1, 1, 1, jpkm1 )
           ! POC disaggregation by turbulence and bacterial activity. 
           ! --------------------------------------------------------
@@ -433,6 +434,6 @@
       IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
          WRITE(charout, FMT="('poc2')")
-         CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+         CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
       ENDIF
       !
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zprod.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zprod.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zprod.F90	(revision 13710)
@@ -16,5 +16,5 @@
    USE sms_pisces      ! PISCES Source Minus Sink variables
    USE p4zlim          ! Co-limitations of differents nutrients
-   USE prtctl_trc      ! print control for debugging
+   USE prtctl          ! print control for debugging
    USE iom             ! I/O manager
 
@@ -48,4 +48,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -109,5 +110,5 @@
       ! day length in hours
       zstrn(:,:) = 0.
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          zargu = TAN( zcodel ) * TAN( gphit(ji,jj) * rad )
          zargu = MAX( -1., MIN(  1., zargu ) )
@@ -116,5 +117,5 @@
 
       ! Impact of the day duration and light intermittency on phytoplankton growth
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
             zval = MAX( 1., zstrn(ji,jj) )
@@ -134,5 +135,5 @@
 
       ! Computation of the P-I slope for nanos and diatoms
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
             ztn         = MAX( 0., ts(ji,jj,jk,jp_tem,Kmm) - 15. )
@@ -149,5 +150,5 @@
       END_3D
 
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
              ! Computation of production function for Carbon
@@ -170,5 +171,5 @@
       !  Computation of a proxy of the N/C ratio
       !  ---------------------------------------
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
           zval = MIN( xnanopo4(ji,jj,jk), ( xnanonh4(ji,jj,jk) + xnanono3(ji,jj,jk) ) )   &
           &      * zprmaxn(ji,jj,jk) / ( zprbio(ji,jj,jk) + rtrn )
@@ -180,5 +181,5 @@
 
 
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
 
           IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
@@ -204,5 +205,5 @@
       !  Sea-ice effect on production
 
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
          zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
@@ -210,5 +211,5 @@
 
       ! Computation of the various production terms 
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
             !  production terms for nanophyto. (C)
@@ -236,5 +237,5 @@
 
       ! Computation of the chlorophyll production terms
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
             !  production terms for nanophyto. ( chlorophyll )
@@ -259,5 +260,5 @@
 
       !   Update the arrays TRA which contain the biological sources and sinks
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
         IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
            zproreg  = zprorcan(ji,jj,jk) - zpronewn(ji,jj,jk)
@@ -287,5 +288,5 @@
      IF( ln_ligand ) THEN
          zpligprod1(:,:,:) = 0._wp    ;    zpligprod2(:,:,:) = 0._wp
-         DO_3D_11_11( 1, jpkm1 )
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
            IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
               zdocprod = excretd * zprorcad(ji,jj,jk) + excretn * zprorcan(ji,jj,jk)
@@ -330,6 +331,6 @@
      IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
          WRITE(charout, FMT="('prod')")
-         CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+         CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
      ENDIF
       !
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zrem.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zrem.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zrem.F90	(revision 13710)
@@ -18,5 +18,5 @@
    USE p4zprod         !  Growth rate of the 2 phyto groups
    USE p4zlim
-   USE prtctl_trc      !  print control for debugging
+   USE prtctl          !  print control for debugging
    USE iom             !  I/O manager
 
@@ -44,4 +44,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -88,5 +89,5 @@
       ! that was modeling explicitely bacteria
       ! -------------------------------------------------------
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          zdep = MAX( hmld(ji,jj), heup(ji,jj) )
          IF( gdept(ji,jj,jk,Kmm) < zdep ) THEN
@@ -102,5 +103,5 @@
 
       IF( ln_p4z ) THEN
-         DO_3D_11_11( 1, jpkm1 )
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
             ! DOC ammonification. Depends on depth, phytoplankton biomass
             ! and a limitation term which is supposed to be a parameterization of the bacterial activity. 
@@ -133,5 +134,5 @@
          END_3D
       ELSE
-         DO_3D_11_11( 1, jpkm1 )
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
             ! DOC ammonification. Depends on depth, phytoplankton biomass
             ! and a limitation term which is supposed to be a parameterization of the bacterial activity. 
@@ -177,5 +178,5 @@
 
 
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          ! NH4 nitrification to NO3. Ceased for oxygen concentrations
          ! below 2 umol/L. Inhibited at strong light 
@@ -195,9 +196,9 @@
        IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
          WRITE(charout, FMT="('rem1')")
-         CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+         CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
        ENDIF
 
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
 
          ! Bacterial uptake of iron. No iron is available in DOC. So
@@ -217,6 +218,6 @@
        IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
          WRITE(charout, FMT="('rem2')")
-         CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+         CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
        ENDIF
 
@@ -225,5 +226,5 @@
       ! ---------------------------------------------------------------
 
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          zdep     = MAX( hmld(ji,jj), heup_01(ji,jj) )
          zsatur   = MAX( rtrn, ( sio3eq(ji,jj,jk) - tr(ji,jj,jk,jpsil,Kbb) ) / ( sio3eq(ji,jj,jk) + rtrn ) )
@@ -248,6 +249,6 @@
       IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
          WRITE(charout, FMT="('rem3')")
-         CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+         CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
        ENDIF
 
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zsed.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zsed.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zsed.F90	(revision 13710)
@@ -18,5 +18,5 @@
    USE sed             !  Sediment module
    USE iom             !  I/O manager
-   USE prtctl_trc      !  print control for debugging
+   USE prtctl          !  print control for debugging
 
    IMPLICIT NONE
@@ -39,4 +39,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -93,5 +94,5 @@
          ! OA: Warning, the following part is necessary to avoid CFL problems above the sediments
          ! --------------------------------------------------------------------
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             ikt  = mbkt(ji,jj)
             zdep = e3t(ji,jj,ikt,Kmm) / xstep
@@ -103,5 +104,5 @@
          ! Computation of the fraction of organic matter that is permanently buried from Dunne's model
          ! -------------------------------------------------------
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
            IF( tmask(ji,jj,1) == 1 ) THEN
               ikt = mbkt(ji,jj)
@@ -129,5 +130,5 @@
       IF( .NOT.lk_sed )  zrivsil = 1._wp - sedsilfrac
 
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          ikt  = mbkt(ji,jj)
          zdep = xstep / e3t(ji,jj,ikt,Kmm) 
@@ -141,5 +142,5 @@
       !
       IF( .NOT.lk_sed ) THEN
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             ikt  = mbkt(ji,jj)
             zdep = xstep / e3t(ji,jj,ikt,Kmm) 
@@ -159,5 +160,5 @@
       ENDIF
       !
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          ikt  = mbkt(ji,jj)
          zdep = xstep / e3t(ji,jj,ikt,Kmm) 
@@ -171,5 +172,5 @@
       !
       IF( ln_p5z ) THEN
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             ikt  = mbkt(ji,jj)
             zdep = xstep / e3t(ji,jj,ikt,Kmm) 
@@ -186,5 +187,5 @@
          ! The 0.5 factor in zpdenit is to avoid negative NO3 concentration after
          ! denitrification in the sediments. Not very clever, but simpliest option.
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             ikt  = mbkt(ji,jj)
             zdep = xstep / e3t(ji,jj,ikt,Kmm) 
@@ -223,5 +224,5 @@
       ENDDO
       IF( ln_p4z ) THEN
-         DO_3D_11_11( 1, jpkm1 )
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
             !                      ! Potential nitrogen fixation dependant on temperature and iron
             ztemp = ts(ji,jj,jk,jp_tem,Kmm)
@@ -239,5 +240,5 @@
          END_3D
       ELSE       ! p5z
-         DO_3D_11_11( 1, jpkm1 )
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
             !                      ! Potential nitrogen fixation dependant on temperature and iron
             ztemp = ts(ji,jj,jk,jp_tem,Kmm)
@@ -260,5 +261,5 @@
       ! ----------------------------------------
       IF( ln_p4z ) THEN
-         DO_3D_11_11( 1, jpkm1 )
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
             zfact = nitrpot(ji,jj,jk) * nitrfix
             tr(ji,jj,jk,jpnh4,Krhs) = tr(ji,jj,jk,jpnh4,Krhs) + zfact / 3.0
@@ -277,5 +278,5 @@
          END_3D
       ELSE    ! p5z
-         DO_3D_11_11( 1, jpkm1 )
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
             zfact = nitrpot(ji,jj,jk) * nitrfix
             tr(ji,jj,jk,jpnh4,Krhs) = tr(ji,jj,jk,jpnh4,Krhs) + zfact / 3.0
@@ -312,8 +313,8 @@
       ENDIF
       !
-      IF(sn_cfctl%l_prttrc) THEN  ! print mean trends (USEd for debugging)
+      IF(sn_cfctl%l_prttrc) THEN  ! print mean trneds (USEd for debugging)
          WRITE(charout, fmt="('sed ')")
-         CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+         CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
       ENDIF
       !
@@ -365,4 +366,6 @@
       lk_sed = ln_sediment .AND. ln_sed_2way 
       !
+      nitrpot(:,:,jpk) = 0._wp   ! define last level for iom_put
+      !
    END SUBROUTINE p4z_sed_init
 
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zsink.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zsink.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zsink.F90	(revision 13710)
@@ -17,5 +17,5 @@
    USE sms_pisces      !  PISCES Source Minus Sink variables
    USE trcsink         !  General routine to compute sedimentation
-   USE prtctl_trc      !  print control for debugging
+   USE prtctl          !  print control for debugging
    USE iom             !  I/O manager
    USE lib_mpp
@@ -40,4 +40,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -80,5 +81,5 @@
       !    by data and from the coagulation theory
       !    -----------------------------------------------------------
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          zmax  = MAX( heup_01(ji,jj), hmld(ji,jj) )
          zfact = MAX( 0., gdepw(ji,jj,jk+1,Kmm) - zmax ) / wsbio2scale
@@ -143,6 +144,6 @@
       IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
          WRITE(charout, FMT="('sink')")
-         CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+         CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
       ENDIF
       !
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zsms.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zsms.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p4zsms.F90	(revision 13710)
@@ -25,5 +25,5 @@
    USE trdtrc          ! TOP trends variables
    USE sedmodel        ! Sediment model
-   USE prtctl_trc      ! print control for debugging
+   USE prtctl          ! print control for debugging
 
    IMPLICIT NONE
@@ -41,4 +41,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -68,5 +69,5 @@
       REAL(wp), ALLOCATABLE, DIMENSION(:,:    ) :: zw2d
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:  ) :: zw3d
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   ztrdt   ! 4D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jp_pisces) :: ztrbbio
 
       !!---------------------------------------------------------------------
@@ -92,11 +93,4 @@
       rfact = rDt_trc
       !
-      ! trends computation initialisation
-      IF( l_trdtrc )  THEN
-         ALLOCATE( ztrdt(jpi,jpj,jpk,jp_pisces) )  !* store now fields before applying the Asselin filter
-         ztrdt(:,:,:,:)  = tr(:,:,:,:,Kmm)
-      ENDIF
-      !
-
       IF( ( ln_top_euler .AND. kt == nittrc000 )  .OR. ( .NOT.ln_top_euler .AND. kt <= nittrc000 + 1 ) ) THEN
          rfactr  = 1. / rfact
@@ -116,4 +110,9 @@
          END DO
       ENDIF
+
+      DO jn = jp_pcs0, jp_pcs1              !   Store the tracer concentrations before entering PISCES
+         ztrbbio(:,:,:,jn) = tr(:,:,:,jn,Kbb)
+      END DO
+
       !
       IF( ll_bc )    CALL p4z_bc( kt, Kbb, Kmm, Krhs )   ! external sources of nutrients 
@@ -132,5 +131,5 @@
          xnegtr(:,:,:) = 1.e0
          DO jn = jp_pcs0, jp_pcs1
-            DO_3D_11_11( 1, jpk )
+            DO_3D( 1, 1, 1, 1, 1, jpk )
                IF( ( tr(ji,jj,jk,jn,Kbb) + tr(ji,jj,jk,jn,Krhs) ) < 0.e0 ) THEN
                   ztra             = ABS( tr(ji,jj,jk,jn,Kbb) ) / ( ABS( tr(ji,jj,jk,jn,Krhs) ) + rtrn )
@@ -197,32 +196,26 @@
          END DO
          !
-         IF( ln_top_euler ) THEN
-            DO jn = jp_pcs0, jp_pcs1
-               tr(:,:,:,jn,Kmm) = tr(:,:,:,jn,Kbb)
-            END DO
-         ENDIF
+      END DO
+      !
+#endif
+      !
+      IF( ln_sediment ) THEN 
+         !
+         CALL sed_model( kt, Kbb, Kmm, Krhs )     !  Main program of Sediment model
+         !
+      ENDIF
+      !
+      DO jn = jp_pcs0, jp_pcs1
+         tr(:,:,:,jn,Krhs) = ( tr(:,:,:,jn,Kbb) - ztrbbio(:,:,:,jn) ) * rfactr
+         tr(:,:,:,jn,Kbb ) = ztrbbio(:,:,:,jn)
+         ztrbbio(:,:,:,jn) = 0._wp
       END DO
       !
       IF( l_trdtrc ) THEN
          DO jn = jp_pcs0, jp_pcs1
-           ztrdt(:,:,:,jn) = ( tr(:,:,:,jn,Kbb) - ztrdt(:,:,:,jn) ) * rfact2r 
            CALL trd_trc( tr(:,:,:,jn,Krhs), jn, jptra_sms, kt, Kmm )   ! save trends
          END DO
-         DEALLOCATE( ztrdt ) 
       END IF
-#endif
-      !
-      IF( ln_sediment ) THEN 
-         !
-         CALL sed_model( kt, Kbb, Kmm, Krhs )     !  Main program of Sediment model
-         !
-         IF( ln_top_euler ) THEN
-            DO jn = jp_pcs0, jp_pcs1
-               tr(:,:,:,jn,Kmm) = tr(:,:,:,jn,Kbb)
-            END DO
-         ENDIF
-         !
-      ENDIF
-      !
+      !  
       IF( lrst_trc )  CALL p4z_rst( kt, Kbb, Kmm,  'WRITE' )           !* Write PISCES informations in restart file 
       !
@@ -340,12 +333,12 @@
          ! 
          IF( iom_varid( numrtr, 'PH', ldstop = .FALSE. ) > 0 ) THEN
-            CALL iom_get( numrtr, jpdom_autoglo, 'PH' , hi(:,:,:)  )
+            CALL iom_get( numrtr, jpdom_auto, 'PH' , hi(:,:,:)  )
          ELSE
             CALL p4z_che( Kbb, Kmm )                  ! initialize the chemical constants
             CALL ahini_for_at( hi, Kbb )
          ENDIF
-         CALL iom_get( numrtr, jpdom_autoglo, 'Silicalim', xksi(:,:) )
+         CALL iom_get( numrtr, jpdom_auto, 'Silicalim', xksi(:,:) )
          IF( iom_varid( numrtr, 'Silicamax', ldstop = .FALSE. ) > 0 ) THEN
-            CALL iom_get( numrtr, jpdom_autoglo, 'Silicamax' , xksimax(:,:)  )
+            CALL iom_get( numrtr, jpdom_auto, 'Silicamax' , xksimax(:,:)  )
          ELSE
             xksimax(:,:) = xksi(:,:)
@@ -360,7 +353,7 @@
          IF( ln_p5z ) THEN
             IF( iom_varid( numrtr, 'sized', ldstop = .FALSE. ) > 0 ) THEN
-               CALL iom_get( numrtr, jpdom_autoglo, 'sizep' , sizep(:,:,:)  )
-               CALL iom_get( numrtr, jpdom_autoglo, 'sizen' , sizen(:,:,:)  )
-               CALL iom_get( numrtr, jpdom_autoglo, 'sized' , sized(:,:,:)  )
+               CALL iom_get( numrtr, jpdom_auto, 'sizep' , sizep(:,:,:)  )
+               CALL iom_get( numrtr, jpdom_auto, 'sizen' , sizen(:,:,:)  )
+               CALL iom_get( numrtr, jpdom_auto, 'sized' , sized(:,:,:)  )
             ELSE
                sizep(:,:,:) = 1.
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p5zlim.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p5zlim.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p5zlim.F90	(revision 13710)
@@ -131,5 +131,5 @@
       zratchl = 6.0
       !
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          ! 
          ! Tuning of the iron concentration to a minimum level that is set to the detection limit
@@ -306,5 +306,5 @@
          &          / (xqndmax(ji,jj,jk) - 2. * xqndmin(ji,jj,jk) ) )   &
          &          * xqndmax(ji,jj,jk) / (zration + rtrn)
-         zlim3    = tr(ji,jj,jk,jpsil,Kbb) / ( tr(ji,jj,jk,jpsil,Kbb) + xksi(ji,jj) )
+         zlim3    = tr(ji,jj,jk,jpsil,Kbb) / ( tr(ji,jj,jk,jpsil,Kbb) + xksi(ji,jj) + rtrn )
          zlim4    = MAX( 0., ( zratiof - zqfemd ) / qfdopt )
          xlimdfe(ji,jj,jk) = MIN( 1., zlim4 )
@@ -318,5 +318,5 @@
       ! phytoplankton (see Daines et al., 2013). 
       ! --------------------------------------------------------------------------------------------------
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          ! Size estimation of nanophytoplankton
          ! ------------------------------------
@@ -367,5 +367,5 @@
       ! Compute the fraction of nanophytoplankton that is made of calcifiers
       ! --------------------------------------------------------------------
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          zlim1 =  tr(ji,jj,jk,jpnh4,Kbb) / ( tr(ji,jj,jk,jpnh4,Kbb) + concnnh4 ) + tr(ji,jj,jk,jpno3,Kbb)    &
          &        / ( tr(ji,jj,jk,jpno3,Kbb) + concnno3 ) * ( 1.0 - tr(ji,jj,jk,jpnh4,Kbb)   &
@@ -385,5 +385,5 @@
       END_3D
       !
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          ! denitrification factor computed from O2 levels
          nitrfac(ji,jj,jk) = MAX(  0.e0, 0.4 * ( 6.e-6  - tr(ji,jj,jk,jpoxy,Kbb) )    &
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p5zmeso.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p5zmeso.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p5zmeso.F90	(revision 13710)
@@ -15,5 +15,5 @@
    USE trc             !  passive tracers common variables 
    USE sms_pisces      !  PISCES Source Minus Sink variables
-   USE prtctl_trc      !  print control for debugging
+   USE prtctl          !  print control for debugging
    USE iom             !  I/O manager
 
@@ -98,5 +98,5 @@
       IF ( bmetexc2 ) zmetexcess = 1.0
 
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          zcompam   = MAX( ( tr(ji,jj,jk,jpmes,Kbb) - 1.e-9 ), 0.e0 )
          zfact     = xstep * tgfunc2(ji,jj,jk) * zcompam
@@ -359,6 +359,6 @@
       IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
         WRITE(charout, FMT="('meso')")
-        CALL prt_ctl_trc_info(charout)
-        CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
+        CALL prt_ctl_info( charout, cdcomp = 'top' )
+        CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
       ENDIF
       !
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p5zmicro.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p5zmicro.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p5zmicro.F90	(revision 13710)
@@ -18,5 +18,5 @@
    USE p5zlim          !  Phytoplankton limitation terms
    USE iom             !  I/O manager
-   USE prtctl_trc      !  print control for debugging
+   USE prtctl          !  print control for debugging
 
    IMPLICIT NONE
@@ -96,5 +96,5 @@
       IF ( bmetexc ) zmetexcess = 1.0
       !
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          zcompaz = MAX( ( tr(ji,jj,jk,jpzoo,Kbb) - 1.e-9 ), 0.e0 )
          zfact   = xstep * tgfunc2(ji,jj,jk) * zcompaz
@@ -306,6 +306,6 @@
       IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
          WRITE(charout, FMT="('micro')")
-         CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+         CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
       ENDIF
       !
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p5zmort.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p5zmort.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p5zmort.F90	(revision 13710)
@@ -16,5 +16,5 @@
    USE p4zlim
    USE p5zlim          !  Phytoplankton limitation terms
-   USE prtctl_trc      !  print control for debugging
+   USE prtctl          !  print control for debugging
 
    IMPLICIT NONE
@@ -82,5 +82,5 @@
       !
       prodcal(:,:,:) = 0.  !: calcite production variable set to zero
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          zcompaph = MAX( ( tr(ji,jj,jk,jpphy,Kbb) - 1e-9 ), 0.e0 )
          !   Squared mortality of Phyto similar to a sedimentation term during
@@ -121,6 +121,6 @@
        IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
          WRITE(charout, FMT="('nano')")
-         CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+         CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
        ENDIF
       !
@@ -148,5 +148,5 @@
       IF( ln_timing )   CALL timing_start('p5z_pico')
       !
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          zcompaph = MAX( ( tr(ji,jj,jk,jppic,Kbb) - 1e-9 ), 0.e0 )
          !  Squared mortality of Phyto similar to a sedimentation term during
@@ -179,6 +179,6 @@
        IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
          WRITE(charout, FMT="('pico')")
-         CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+         CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
        ENDIF
       !
@@ -207,5 +207,5 @@
       !
 
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
 
          zcompadi = MAX( ( tr(ji,jj,jk,jpdia,Kbb) - 1E-9), 0. )
@@ -254,6 +254,6 @@
       IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
          WRITE(charout, FMT="('diat')")
-         CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+         CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
       ENDIF
       !
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p5zprod.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p5zprod.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/P4Z/p5zprod.F90	(revision 13710)
@@ -18,5 +18,5 @@
    USE p4zlim
    USE p5zlim          !  Co-limitations of differents nutrients
-   USE prtctl_trc      !  print control for debugging
+   USE prtctl          !  print control for debugging
    USE iom             !  I/O manager
 
@@ -52,4 +52,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -124,5 +125,5 @@
       ! day length in hours
       zstrn(:,:) = 0.
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          zargu = TAN( zcodel ) * TAN( gphit(ji,jj) * rad )
          zargu = MAX( -1., MIN(  1., zargu ) )
@@ -131,5 +132,5 @@
 
          ! Impact of the day duration on phytoplankton growth
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
             zval = MAX( 1., zstrn(ji,jj) )
@@ -151,5 +152,5 @@
       WHERE( zstrn(:,:) < 1.e0 ) zstrn(:,:) = 24.
 
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
             ! Computation of the P-I slope for nanos and diatoms
@@ -185,5 +186,5 @@
       END_3D
 
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
 
           IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
@@ -207,5 +208,5 @@
 
       !  Sea-ice effect on production                                                                               
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          zprbio(ji,jj,jk)  = zprbio(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
          zprpic(ji,jj,jk)  = zprpic(ji,jj,jk) * ( 1. - fr_i(ji,jj) ) 
@@ -215,5 +216,5 @@
 
       ! Computation of the various production terms of nanophytoplankton 
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
             !  production terms for nanophyto.
@@ -248,5 +249,5 @@
 
       ! Computation of the various production terms of picophytoplankton 
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
             !  production terms for picophyto.
@@ -281,5 +282,5 @@
 
       ! Computation of the various production terms of diatoms
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
             !  production terms for diatomees
@@ -315,5 +316,5 @@
       END_3D
 
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
          IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
                !  production terms for nanophyto. ( chlorophyll )
@@ -346,5 +347,5 @@
 
       !   Update the arrays TRA which contain the biological sources and sinks
-      DO_3D_11_11( 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
         zprontot = zpronewn(ji,jj,jk) + zproregn(ji,jj,jk)
         zproptot = zpronewp(ji,jj,jk) + zproregp(ji,jj,jk)
@@ -409,5 +410,5 @@
      IF( ln_ligand ) THEN
          zpligprod1(:,:,:) = 0._wp    ;    zpligprod2(:,:,:) = 0._wp             
-         DO_3D_11_11( 1, jpkm1 )
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
            zdocprod = excretd * zprorcad(ji,jj,jk) + excretn * zprorcan(ji,jj,jk) + excretp * zprorcap(ji,jj,jk)
            zfeup    = texcretn * zprofen(ji,jj,jk) + texcretd * zprofed(ji,jj,jk) + texcretp * zprofep(ji,jj,jk)
@@ -460,6 +461,6 @@
       IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
          WRITE(charout, FMT="('prod')")
-         CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+         CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
       ENDIF
       !
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/SED/oce_sed.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/SED/oce_sed.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/SED/oce_sed.F90	(revision 13710)
@@ -13,5 +13,8 @@
    USE dom_oce , ONLY :   glamt     =>   glamt          !: longitude of t-point (degre)
    USE dom_oce , ONLY :   gphit     =>   gphit          !: latitude  of t-point (degre)
+!!st 
+#if ! defined key_qco
    USE dom_oce , ONLY :   e3t       =>   e3t            !: latitude  of t-point (degre)
+#endif
    USE dom_oce , ONLY :   e3t_1d    =>   e3t_1d         !: reference depth of t-points (m)
    USE dom_oce , ONLY :   gdepw_0   =>   gdepw_0        !: reference depth of t-points (m)
@@ -53,4 +56,2 @@
 
 END MODULE oce_sed
-
-
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/SED/sedchem.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/SED/sedchem.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/SED/sedchem.F90	(revision 13710)
@@ -138,5 +138,5 @@
          CALL sed_chem_cst
       ELSE
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             ikt = mbkt(ji,jj) 
             IF ( tmask(ji,jj,ikt) == 1 ) THEN
@@ -577,5 +577,5 @@
          saltprac(:) = salt(:) * 35.0 / 35.16504
       ELSE
-         saltprac(:) = temp(:)
+         saltprac(:) = salt(:)
       ENDIF
 
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/SED/seddta.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/SED/seddta.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/SED/seddta.F90	(revision 13710)
@@ -24,4 +24,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !! $Id$
 CONTAINS
@@ -95,5 +96,5 @@
       !    -----------------------------------------------------------
       IF (ln_sediment_offline) THEN
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             ikt = mbkt(ji,jj)
             zwsbio4(ji,jj) = wsbio2 / rday
@@ -101,5 +102,5 @@
          END_2D
       ELSE
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             ikt = mbkt(ji,jj)
             zdep = e3t(ji,jj,ikt,Kmm) / rDt_trc
@@ -110,5 +111,5 @@
 
       trc_data(:,:,:) = 0.
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          ikt = mbkt(ji,jj)
          IF ( tmask(ji,jj,ikt) == 1 ) THEN
@@ -164,5 +165,5 @@
       CALL pack_arr ( jpoce,  rainrm_dta(1:jpoce,jscal), trc_data(1:jpi,1:jpj,14), iarroce(1:jpoce) )
       rainrm_dta(1:jpoce,jscal) = rainrm_dta(1:jpoce,jscal) * 1e-4
-      ! vector temperature [°C] and salinity 
+      ! vector temperature [�C] and salinity 
       CALL pack_arr ( jpoce,  temp(1:jpoce), trc_data(1:jpi,1:jpj,15), iarroce(1:jpoce) )
       CALL pack_arr ( jpoce,  salt(1:jpoce), trc_data(1:jpi,1:jpj,16), iarroce(1:jpoce) )
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/SED/sedini.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/SED/sedini.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/SED/sedini.F90	(revision 13710)
@@ -135,5 +135,5 @@
       ! Determination of sediments number of points and allocate global variables
       epkbot(:,:) = 0.
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          ikt = mbkt(ji,jj) 
          IF( tmask(ji,jj,ikt) == 1 ) epkbot(ji,jj) = e3t_1d(ikt)
@@ -247,5 +247,5 @@
       ! Computation of 1D array of sediments points
       indoce = 0
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          IF (  epkbot(ji,jj) > 0. ) THEN
             indoce          = indoce + 1
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/SED/sedinorg.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/SED/sedinorg.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/SED/sedinorg.F90	(revision 13710)
@@ -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_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/SED/sedrst.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/SED/sedrst.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/SED/sedrst.F90	(revision 13710)
@@ -123,5 +123,5 @@
          cltra = TRIM(sedtrcd(jn))
          IF( iom_varid( numrsr, TRIM(cltra) , ldstop = .FALSE. ) > 0 ) THEN
-            CALL iom_get( numrsr, jpdom_autoglo, TRIM(cltra), zdta(:,:,:,jn) )
+            CALL iom_get( numrsr, jpdom_auto, TRIM(cltra), zdta(:,:,:,jn) )
          ELSE
             zdta(:,:,:,jn) = 0.0
@@ -142,5 +142,5 @@
          cltra = TRIM(seddia3d(jn))
          IF( iom_varid( numrsr, TRIM(cltra) , ldstop = .FALSE. ) > 0 ) THEN
-            CALL iom_get( numrsr, jpdom_autoglo, TRIM(cltra), zdta1(:,:,:,jn) )
+            CALL iom_get( numrsr, jpdom_auto, TRIM(cltra), zdta1(:,:,:,jn) )
          ELSE
             zdta1(:,:,:,jn) = 0.0
@@ -169,5 +169,5 @@
       cltra = "dbioturb"
       IF( iom_varid( numrsr, TRIM(cltra) , ldstop = .FALSE. ) > 0 ) THEN
-         CALL iom_get( numrsr, jpdom_autoglo, TRIM(cltra), zdta2(:,:,:) )
+         CALL iom_get( numrsr, jpdom_auto, TRIM(cltra), zdta2(:,:,:) )
       ELSE
          zdta2(:,:,:) = 0.0
@@ -179,5 +179,5 @@
       cltra = "irrig"
       IF( iom_varid( numrsr, TRIM(cltra) , ldstop = .FALSE. ) > 0 ) THEN
-         CALL iom_get( numrsr, jpdom_autoglo, TRIM(cltra), zdta2(:,:,:) )
+         CALL iom_get( numrsr, jpdom_auto, TRIM(cltra), zdta2(:,:,:) )
       ELSE
          zdta2(:,:,:) = 0.0
@@ -189,5 +189,5 @@
       cltra = "sedligand"
       IF( iom_varid( numrsr, TRIM(cltra) , ldstop = .FALSE. ) > 0 ) THEN
-         CALL iom_get( numrsr, jpdom_autoglo, TRIM(cltra), zdta2(:,:,:) )
+         CALL iom_get( numrsr, jpdom_auto, TRIM(cltra), zdta2(:,:,:) )
       ELSE
          zdta2(:,:,:) = 0.0
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/SED/sedsfc.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/SED/sedsfc.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/SED/sedsfc.F90	(revision 13710)
@@ -48,5 +48,5 @@
 
 
-      DO_2D_11_11
+      DO_2D( 1, 1, 1, 1 )
          ikt = mbkt(ji,jj)
          IF ( tmask(ji,jj,ikt) == 1 ) THEN
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/SED/trcdmp_sed.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/SED/trcdmp_sed.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/SED/trcdmp_sed.F90	(revision 13710)
@@ -21,5 +21,5 @@
    USE trc             ! ocean passive tracers variables
    USE trcdta
-   USE prtctl_trc      ! Print control for debbuging
+   USE prtctl          ! Print control for debbuging
    USE iom
 
@@ -93,5 +93,5 @@
                CALL trc_dta( kt, Kmm, sf_trcdta(jl), rf_trfac(jl), ztrcdta )   ! read tracer data at nit000
                !
-               DO_2D_11_11
+               DO_2D( 1, 1, 1, 1 )
                   ikt = mbkt(ji,jj)
                   tr(ji,jj,ikt,jn,Kbb) = ztrcdta(ji,jj,ikt) + ( tr(ji,jj,ikt,jn,Kbb) -  ztrcdta(ji,jj,ikt) )     &
@@ -107,6 +107,6 @@
       IF( sn_cfctl%l_prttrc ) THEN
          WRITE(charout, FMT="('dmp ')")
-         CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc( tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+         CALL prt_ctl( tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm, clinfo3='trd' )
       ENDIF
       !
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/trcwri_pisces.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/trcwri_pisces.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/PISCES/trcwri_pisces.F90	(revision 13710)
@@ -21,4 +21,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -68,5 +69,5 @@
             zo2min   (:,:) = tr(:,:,1,jpoxy,Kmm) * tmask(:,:,1)
             zdepo2min(:,:) = gdepw(:,:,1,Kmm)   * tmask(:,:,1)
-            DO_3D_11_11( 2, jpkm1 )
+            DO_3D( 1, 1, 1, 1, 2, jpkm1 )
                IF( tmask(ji,jj,jk) == 1 ) then
                   IF( tr(ji,jj,jk,jpoxy,Kmm) < zo2min(ji,jj) ) then
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/TRP/trcadv.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/TRP/trcadv.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/TRP/trcadv.F90	(revision 13710)
@@ -29,5 +29,5 @@
    USE ldfslp         ! Lateral diffusion: slopes of neutral surfaces
    !
-   USE prtctl_trc     ! control print
+   USE prtctl         ! control print
    USE timing         ! Timing
 
@@ -59,4 +59,5 @@
    INTEGER, PARAMETER ::   np_QCK     = 5   ! QUICK scheme
    
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -137,6 +138,6 @@
       IF( sn_cfctl%l_prttrc ) THEN        !== print mean trends (used for debugging)
          WRITE(charout, FMT="('adv ')")
-         CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc( tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+         CALL prt_ctl( tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm, clinfo3='trd' )
       END IF
       !
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/TRP/trcatf.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/TRP/trcatf.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/TRP/trcatf.F90	(revision 13710)
@@ -31,5 +31,9 @@
    USE trd_oce
    USE trdtra
+# if defined key_qco
+   USE traatfqco
+# else
    USE traatf
+# endif
    USE bdy_oce   , ONLY: ln_bdy
    USE trcbdy          ! BDY open boundaries
@@ -39,5 +43,5 @@
    !
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE prtctl_trc      ! Print control for debbuging
+   USE prtctl          ! Print control for debbuging
 
    IMPLICIT NONE
@@ -50,4 +54,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -112,5 +117,5 @@
          ! total trend for the non-time-filtered variables. 
          zfact = 1.0 / rn_Dt
-         ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn; e3tn cancel from ts(Kmm) terms
+         ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3ta*Ta)/e3tn; e3tn cancel from ts(Kmm) terms
          IF( ln_linssh ) THEN       ! linear sea surface height only
             DO jn = 1, jptra
@@ -151,6 +156,11 @@
       ELSE     
          IF( .NOT. l_offline ) THEN ! Leap-Frog + Asselin filter time stepping
+# if defined key_qco
+            IF( ln_linssh ) THEN   ;   CALL tra_atf_fix_lf( kt, Kbb, Kmm, Kaa, nittrc000,        'TRC', ptr, jptra )                     !     linear ssh
+            ELSE                   ;   CALL tra_atf_qco_lf( kt, Kbb, Kmm, Kaa, nittrc000, rn_Dt, 'TRC', ptr, sbc_trc, sbc_trc_b, jptra ) ! non-linear ssh
+# else
             IF( ln_linssh ) THEN   ;   CALL tra_atf_fix( kt, Kbb, Kmm, Kaa, nittrc000,         'TRC', ptr, jptra )                     !     linear ssh
             ELSE                   ;   CALL tra_atf_vvl( kt, Kbb, Kmm, Kaa, nittrc000, rn_Dt, 'TRC', ptr, sbc_trc, sbc_trc_b, jptra ) ! non-linear ssh
+# endif
             ENDIF
          ELSE
@@ -174,6 +184,6 @@
       IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
          WRITE(charout, FMT="('nxt')")
-         CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc(tab4d=ptr(:,:,:,:,Kmm), mask=tmask, clinfo=ctrcnm)
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+         CALL prt_ctl(tab4d_1=ptr(:,:,:,:,Kmm), mask1=tmask, clinfo=ctrcnm)
       ENDIF
       !
@@ -182,5 +192,5 @@
    END SUBROUTINE trc_atf
 
-
+# if ! defined key_qco
    SUBROUTINE trc_atf_off( kt, Kbb, Kmm, Kaa, ptr )
       !!----------------------------------------------------------------------
@@ -198,8 +208,8 @@
       !!                This can be summurized for tempearture as:
       !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T
-      !!                  /( e3t(:,:,:,Kmm)    + rbcp*[ e3t(:,:,:,Kbb)    - 2 e3t(:,:,:,Kmm)    + e3t(:,:,:,Kaa)    ] )   
+      !!                  /( e3t(:,:,jk,Kmm)    + rbcp*[ e3t(:,:,jk,Kbb)    - 2 e3t(:,:,jk,Kmm)    + e3t(:,:,jk,Kaa)    ] )   
       !!             ztm = 0                                                       otherwise
       !!             tb  = ( e3t_n*tn + rn_atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
-      !!                  /( e3t(:,:,:,Kmm)    + rn_atfp*[ e3t(:,:,:,Kbb)    - 2 e3t(:,:,:,Kmm)    + e3t(:,:,:,Kaa)    ] )
+      !!                  /( e3t(:,:,jk,Kmm)    + rn_atfp*[ e3t(:,:,jk,Kbb)    - 2 e3t(:,:,jk,Kmm)    + e3t(:,:,jk,Kaa)    ] )
       !!             tn  = ta 
       !!             ta  = zt        (NB: reset to 0 after eos_bn2 call)
@@ -229,5 +239,5 @@
       !
       DO jn = 1, jptra      
-         DO_3D_11_11( 1, jpkm1 )
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
             ze3t_b = e3t(ji,jj,jk,Kbb)
             ze3t_n = e3t(ji,jj,jk,Kmm)
@@ -257,4 +267,78 @@
       !
    END SUBROUTINE trc_atf_off
+# else
+   SUBROUTINE trc_atf_off( kt, Kbb, Kmm, Kaa, ptr )
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE tra_atf_off  ***
+      !!
+      !!          !!!!!!!!!!!!!!!!! REWRITE HEADER COMMENTS !!!!!!!!!!!!!!
+      !!
+      !! ** Purpose :   Time varying volume: apply the Asselin time filter  
+      !! 
+      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
+      !!              - save in (ta,sa) a thickness weighted average over the three 
+      !!             time levels which will be used to compute rdn and thus the semi-
+      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
+      !!              - swap tracer fields to prepare the next time_step.
+      !!                This can be summurized for tempearture as:
+      !!             ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )   ln_dynhpg_imp = T
+      !!                  /( e3t(:,:,jk,Kmm)    + rbcp*[ e3t(:,:,jk,Kbb)    - 2 e3t(:,:,jk,Kmm)    + e3t(:,:,jk,Kaa)    ] )   
+      !!             ztm = 0                                                       otherwise
+      !!             tb  = ( e3t_n*tn + rn_atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
+      !!                  /( e3t(:,:,jk,Kmm)    + rn_atfp*[ e3t(:,:,jk,Kbb)    - 2 e3t(:,:,jk,Kmm)    + e3t(:,:,jk,Kaa)    ] )
+      !!             tn  = ta 
+      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call)
+      !!
+      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
+      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
+      !!----------------------------------------------------------------------
+      INTEGER                                   , INTENT(in   ) ::  kt            ! ocean time-step index
+      INTEGER                                   , INTENT(in   ) ::  Kbb, Kmm, Kaa ! time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) ::  ptr           ! passive tracers
+      !!     
+      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
+      REAL(wp) ::   ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
+      REAL(wp) ::   ze3t_b, ze3t_n, ze3t_a, ze3t_f           !   -      -
+      !!----------------------------------------------------------------------
+      !
+      IF( kt == nittrc000 )  THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'trc_atf_off : Asselin time filtering'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
+         IF( .NOT. ln_linssh ) THEN
+            rfact1 = rn_atfp * rn_Dt
+            rfact2 = rfact1 / rho0
+         ENDIF
+        !  
+      ENDIF
+      !
+      DO jn = 1, jptra      
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+            ze3t_b = 1._wp + r3t(ji,jj,Kbb) * tmask(ji,jj,jk)
+            ze3t_n = 1._wp + r3t(ji,jj,Kmm) * tmask(ji,jj,jk)
+            ze3t_a = 1._wp + r3t(ji,jj,Kaa) * tmask(ji,jj,jk)
+            !                                         ! tracer content at Before, now and after
+            ztc_b  = ptr(ji,jj,jk,jn,Kbb) * ze3t_b
+            ztc_n  = ptr(ji,jj,jk,jn,Kmm) * ze3t_n
+            ztc_a  = ptr(ji,jj,jk,jn,Kaa) * ze3t_a
+            !
+            ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
+            !
+            ze3t_f = 1._wp + r3t_f(ji,jj)*tmask(ji,jj,jk)
+            ztc_f  = ztc_n  + rn_atfp * ztc_d
+            !
+            IF( .NOT. ln_linssh .AND. jk == mikt(ji,jj) ) THEN           ! first level 
+               ztc_f  = ztc_f  - rfact1 * ( sbc_trc(ji,jj,jn) - sbc_trc_b(ji,jj,jn) )
+            ENDIF
+
+            ze3t_f = 1.e0 / ze3t_f
+            ptr(ji,jj,jk,jn,Kmm) = ztc_f * ze3t_f     ! time filtered "now" field
+            !
+         END_3D
+         ! 
+      END DO
+      !
+   END SUBROUTINE trc_atf_off
+# endif
 
 #else
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/TRP/trcbbl.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/TRP/trcbbl.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/TRP/trcbbl.F90	(revision 13710)
@@ -25,5 +25,5 @@
    USE trdtra         ! tracer trends
    USE trabbl         ! bottom boundary layer 
-   USE prtctl_trc     ! Print control for debbuging
+   USE prtctl         ! Print control for debbuging
 
    PUBLIC   trc_bbl   !  routine called by trctrp.F90
@@ -70,6 +70,6 @@
          CALL tra_bbl_dif( ptr(:,:,:,:,Kbb), ptr(:,:,:,:,Krhs), jptra, Kmm )  
          IF( sn_cfctl%l_prttrc )   THEN
-            WRITE(charout, FMT="(' bbl_dif')")  ;  CALL prt_ctl_trc_info(charout)
-            CALL prt_ctl_trc( tab4d=ptr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
+            WRITE(charout, FMT="(' bbl_dif')")  ;  CALL prt_ctl_info( charout, cdcomp = 'top' )
+            CALL prt_ctl( tab4d_1=ptr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm, clinfo3='trd' )
          ENDIF
          !
@@ -81,6 +81,6 @@
          CALL tra_bbl_adv( ptr(:,:,:,:,Kbb), ptr(:,:,:,:,Krhs), jptra, Kmm )  
          IF( sn_cfctl%l_prttrc )   THEN
-            WRITE(charout, FMT="(' bbl_adv')")  ;  CALL prt_ctl_trc_info(charout)
-            CALL prt_ctl_trc( tab4d=ptr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
+            WRITE(charout, FMT="(' bbl_adv')")  ;  CALL prt_ctl_info( charout, cdcomp = 'top' )
+            CALL prt_ctl( tab4d_1=ptr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm, clinfo3='trd' )
          ENDIF
          !
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/TRP/trcdmp.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/TRP/trcdmp.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/TRP/trcdmp.F90	(revision 13710)
@@ -24,5 +24,5 @@
    !
    USE iom
-   USE prtctl_trc      ! Print control for debbuging
+   USE prtctl          ! Print control for debbuging
 
    IMPLICIT NONE
@@ -45,4 +45,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -112,10 +113,10 @@
                !
                CASE( 0 )                !==  newtonian damping throughout the water column  ==!
-                  DO_3D_00_00( 1, jpkm1 )
+                  DO_3D( 0, 0, 0, 0, 1, jpkm1 )
                      ptr(ji,jj,jk,jn,Krhs) = ptr(ji,jj,jk,jn,Krhs) + restotr(ji,jj,jk) * ( ztrcdta(ji,jj,jk) - ptr(ji,jj,jk,jn,Kbb) )
                   END_3D
                   !
                CASE ( 1 )                !==  no damping in the turbocline (avt > 5 cm2/s)  ==!
-                  DO_3D_00_00( 1, jpkm1 )
+                  DO_3D( 0, 0, 0, 0, 1, jpkm1 )
                      IF( avt(ji,jj,jk) <= avt_c )  THEN 
                         ptr(ji,jj,jk,jn,Krhs) = ptr(ji,jj,jk,jn,Krhs) + restotr(ji,jj,jk) * ( ztrcdta(ji,jj,jk) - ptr(ji,jj,jk,jn,Kbb) )
@@ -124,5 +125,5 @@
                   !
                CASE ( 2 )               !==  no damping in the mixed layer   ==! 
-                  DO_3D_00_00( 1, jpkm1 )
+                  DO_3D( 0, 0, 0, 0, 1, jpkm1 )
                      IF( gdept(ji,jj,jk,Kmm) >= hmlp (ji,jj) ) THEN
                         ptr(ji,jj,jk,jn,Krhs) = ptr(ji,jj,jk,jn,Krhs) + restotr(ji,jj,jk) * ( ztrcdta(ji,jj,jk) - ptr(ji,jj,jk,jn,Kbb) )
@@ -148,6 +149,6 @@
       IF( sn_cfctl%l_prttrc ) THEN
          WRITE(charout, FMT="('dmp ')")
-         CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc( tab4d=ptr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+         CALL prt_ctl( tab4d_1=ptr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm, clinfo3='trd' )
       ENDIF
       !
@@ -204,5 +205,5 @@
          !Read in mask from file
          CALL iom_open ( cn_resto_tr, imask)
-         CALL iom_get  ( imask, jpdom_autoglo, 'resto', restotr)
+         CALL iom_get  ( imask, jpdom_auto, 'resto', restotr)
          CALL iom_close( imask )
          !
@@ -245,32 +246,33 @@
             !                                           ! =======================
             CASE ( 1 )                                  ! eORCA_R1 configuration
-            !                                           ! =======================
-            isrow = 332 - jpjglo
-            !
-            nctsi1(1)   = 333  ; nctsj1(1)   = 243 - isrow   ! Caspian Sea
-            nctsi2(1)   = 342  ; nctsj2(1)   = 274 - isrow
-            !                                        
-            nctsi1(2)   = 198  ; nctsj1(2)   = 258 - isrow   ! Lake Superior
-            nctsi2(2)   = 204  ; nctsj2(2)   = 262 - isrow
-            !                                         
-            nctsi1(3)   = 201  ; nctsj1(3)   = 250 - isrow   ! Lake Michigan
-            nctsi2(3)   = 203  ; nctsj2(3)   = 256 - isrow
-            !                                        
-            nctsi1(4)   = 204  ; nctsj1(4)   = 252 - isrow   ! Lake Huron
-            nctsi2(4)   = 209  ; nctsj2(4)   = 256 - isrow
-            !                                        
-            nctsi1(5)   = 206  ; nctsj1(5)   = 249 - isrow   ! Lake Erie
-            nctsi2(5)   = 209  ; nctsj2(5)   = 251 - isrow
-            !                                        
-            nctsi1(6)   = 210  ; nctsj1(6)   = 252 - isrow   ! Lake Ontario
-            nctsi2(6)   = 212  ; nctsj2(6)   = 252 - isrow
-            !                                        
-            nctsi1(7)   = 321  ; nctsj1(7)   = 180 - isrow   ! Victoria Lake
-            nctsi2(7)   = 322  ; nctsj2(7)   = 189 - isrow
-            !                                        
-            nctsi1(8)   = 297  ; nctsj1(8)   = 270 - isrow   ! Baltic Sea
-            nctsi2(8)   = 308  ; nctsj2(8)   = 293 - isrow
-            !                                        
-            !                                           ! =======================
+               !                                        ! =======================
+               !
+               isrow = 332 - (Nj0glo + 1)   ! was 332 - jpjglo -> jpjglo_old_version = Nj0glo + 1
+               !
+               nctsi1(1)   = 333  ; nctsj1(1)   = 243 - isrow   ! Caspian Sea
+               nctsi2(1)   = 342  ; nctsj2(1)   = 274 - isrow
+               !                                        
+               nctsi1(2)   = 198  ; nctsj1(2)   = 258 - isrow   ! Lake Superior
+               nctsi2(2)   = 204  ; nctsj2(2)   = 262 - isrow
+               !                                         
+               nctsi1(3)   = 201  ; nctsj1(3)   = 250 - isrow   ! Lake Michigan
+               nctsi2(3)   = 203  ; nctsj2(3)   = 256 - isrow
+               !                                        
+               nctsi1(4)   = 204  ; nctsj1(4)   = 252 - isrow   ! Lake Huron
+               nctsi2(4)   = 209  ; nctsj2(4)   = 256 - isrow
+               !                                        
+               nctsi1(5)   = 206  ; nctsj1(5)   = 249 - isrow   ! Lake Erie
+               nctsi2(5)   = 209  ; nctsj2(5)   = 251 - isrow
+               !                                        
+               nctsi1(6)   = 210  ; nctsj1(6)   = 252 - isrow   ! Lake Ontario
+               nctsi2(6)   = 212  ; nctsj2(6)   = 252 - isrow
+               !                                        
+               nctsi1(7)   = 321  ; nctsj1(7)   = 180 - isrow   ! Victoria Lake
+               nctsi2(7)   = 322  ; nctsj2(7)   = 189 - isrow
+               !                                        
+               nctsi1(8)   = 297  ; nctsj1(8)   = 270 - isrow   ! Baltic Sea
+               nctsi2(8)   = 308  ; nctsj2(8)   = 293 - isrow
+               !
+               !                                        ! =======================
             CASE ( 2 )                                  !  ORCA_R2 configuration
                !                                        ! =======================
@@ -285,9 +287,10 @@
                nctsi2(3)   = 181  ;  nctsj2(3)   = 112
               !                                      
-               nctsi1(4)   =   2  ;  nctsj1(4)   = 107      ! Black Sea 2 : est part of the Black Sea
+               nctsi1(4)   =   2  ;  nctsj1(4)   = 107       ! Black Sea 2 : est part of the Black Sea
                nctsi2(4)   =   6  ;  nctsj2(4)   = 112
                !                                     
                nctsi1(5)   =  145 ;  nctsj1(5)   = 116       ! Baltic Sea
                nctsi2(5)   =  150 ;  nctsj2(5)   = 126
+               !
                !                                        ! =======================
             CASE ( 4 )                                  !  ORCA_R4 configuration
@@ -305,4 +308,5 @@
                nctsi1(4)   = 75  ;  nctsj1(4)   = 59         ! Baltic Sea
                nctsi2(4)   = 76  ;  nctsj2(4)   = 61
+               !
                !                                        ! =======================
             CASE ( 025 )                                ! ORCA_R025 configuration
@@ -318,4 +322,7 @@
             !
          ENDIF
+         !
+         nctsi1(:) = nctsi1(:) + nn_hls - 1   ;   nctsi2(:) = nctsi2(:) + nn_hls - 1   ! -1 as x-perio included in old input files
+         nctsj1(:) = nctsj1(:) + nn_hls       ;   nctsj2(:) = nctsj2(:) + nn_hls
          !
          ! convert the position in local domain indices
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/TRP/trcldf.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/TRP/trcldf.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/TRP/trcldf.F90	(revision 13710)
@@ -25,5 +25,5 @@
    USE trdtra         ! trends manager: tracers
    !
-   USE prtctl_trc     ! Print control
+   USE prtctl         ! Print control
 
    IMPLICIT NONE
@@ -44,4 +44,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -81,5 +82,5 @@
       zahv(:,:,:) = rldf * ahtv(:,:,:)
       !                                  !* Enhanced zonal diffusivity coefficent in the equatorial domain
-      DO_3D_11_11( 1, jpk )
+      DO_3D( 1, 1, 1, 1, 1, jpk )
          IF( gdept(ji,jj,jk,Kmm) > 200. .AND. gphit(ji,jj) < 5. .AND. gphit(ji,jj) > -5. ) THEN
             zdep = MAX( gdept(ji,jj,jk,Kmm) - 1000., 0. ) / 1000.
@@ -114,6 +115,6 @@
       IF( sn_cfctl%l_prttrc ) THEN ! print mean trends (used for debugging)
          WRITE(charout, FMT="('ldf ')")
-         CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc( tab4d=ptr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+         CALL prt_ctl( tab4d_1=ptr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm, clinfo3='trd' )
       ENDIF
       !
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/TRP/trcrad.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/TRP/trcrad.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/TRP/trcrad.F90	(revision 13710)
@@ -19,5 +19,5 @@
    USE trd_oce
    USE trdtra
-   USE prtctl_trc          ! Print control for debbuging
+   USE prtctl              ! Print control for debbuging
    USE lib_fortran
 
@@ -72,6 +72,6 @@
       IF(sn_cfctl%l_prttrc) THEN      ! print mean trends (used for debugging)
          WRITE(charout, FMT="('rad')")
-         CALL prt_ctl_trc_info( charout )
-         CALL prt_ctl_trc( tab4d=ptr(:,:,:,:,Kbb), mask=tmask, clinfo=ctrcnm )
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+         CALL prt_ctl( tab4d_1=ptr(:,:,:,:,Kbb), mask1=tmask, clinfo=ctrcnm )
       ENDIF
       !
@@ -168,5 +168,5 @@
               IF( l_trdtrc )   ztrtrd(:,:,:) = ptr(:,:,:,jn,itime)                       ! save input tr(:,:,:,:,Kbb) for trend computation           
               !
-              DO_3D_11_11( 1, jpkm1 )
+              DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
                  IF( ztrneg(ji,jj,jn) /= 0. ) THEN                                 ! if negative values over the 3x3 box
                     !
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/TRP/trcsbc.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/TRP/trcsbc.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/TRP/trcsbc.F90	(revision 13710)
@@ -18,5 +18,5 @@
    USE oce_trc         ! ocean dynamics and active tracers variables
    USE trc             ! ocean  passive tracers variables
-   USE prtctl_trc      ! Print control for debbuging
+   USE prtctl          ! Print control for debbuging
    USE iom
    USE trd_oce
@@ -30,4 +30,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -49,5 +50,5 @@
       !!            The surface freshwater flux modify the ocean volume
       !!         and thus the concentration of a tracer as :
-      !!            tr(Krhs) = tr(Krhs) + emp * tr(Kmm) / e3t   for k=1
+      !!            tr(Krhs) = tr(Krhs) + emp * tr(Kmm) / e3t_   for k=1
       !!         where emp, the surface freshwater budget (evaporation minus
       !!         precipitation ) given in kg/m2/s is divided
@@ -87,5 +88,5 @@
             zfact = 0.5_wp
             DO jn = 1, jptra
-               CALL iom_get( numrtr, jpdom_autoglo, 'sbc_'//TRIM(ctrcnm(jn))//'_b', sbc_trc_b(:,:,jn) )   ! before tracer content sbc
+               CALL iom_get( numrtr, jpdom_auto, 'sbc_'//TRIM(ctrcnm(jn))//'_b', sbc_trc_b(:,:,jn) )   ! before tracer content sbc
             END DO
          ELSE                                         ! No restart or restart not found: Euler forward time stepping
@@ -120,5 +121,5 @@
          !
          DO jn = 1, jptra
-            DO_2D_01_00
+            DO_2D( 0, 1, 0, 0 )
                sbc_trc(ji,jj,jn) = zsfx(ji,jj) * r1_rho0 * ptr(ji,jj,1,jn,Kmm)
             END_2D
@@ -128,5 +129,5 @@
          !
          DO jn = 1, jptra
-            DO_2D_01_00
+            DO_2D( 0, 1, 0, 0 )
                sbc_trc(ji,jj,jn) = ( zsfx(ji,jj) + fmmflx(ji,jj) ) * r1_rho0 * ptr(ji,jj,1,jn,Kmm)
             END_2D
@@ -136,5 +137,5 @@
          !
          DO jn = 1, jptra
-            DO_2D_01_00
+            DO_2D( 0, 1, 0, 0 )
                zse3t = 1. / e3t(ji,jj,1,Kmm)
                ! tracer flux at the ice/ocean interface (tracer/m2/s)
@@ -154,5 +155,5 @@
       END SELECT
       !
-      CALL lbc_lnk( 'trcsbc', sbc_trc(:,:,:), 'T', 1. )
+      CALL lbc_lnk( 'trcsbc', sbc_trc(:,:,:), 'T', 1.0_wp )
       !                                       Concentration dilution effect on tracers due to evaporation & precipitation 
       DO jn = 1, jptra
@@ -160,5 +161,5 @@
          IF( l_trdtrc )   ztrtrd(:,:,:) = ptr(:,:,:,jn,Krhs)  ! save trends
          !
-         DO_2D_01_00
+         DO_2D( 0, 1, 0, 0 )
             zse3t = zfact / e3t(ji,jj,1,Kmm)
             ptr(ji,jj,1,jn,Krhs) = ptr(ji,jj,1,jn,Krhs) + ( sbc_trc_b(ji,jj,jn) + sbc_trc(ji,jj,jn) ) * zse3t
@@ -186,6 +187,6 @@
       !
       IF( sn_cfctl%l_prttrc )   THEN
-         WRITE(charout, FMT="('sbc ')") ;  CALL prt_ctl_trc_info(charout)
-                                           CALL prt_ctl_trc( tab4d=ptr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
+         WRITE(charout, FMT="('sbc ')") ;  CALL prt_ctl_info( charout, cdcomp = 'top' )
+                                           CALL prt_ctl( tab4d_1=ptr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm, clinfo3='trd' )
       ENDIF
       IF( l_trdtrc )  DEALLOCATE( ztrtrd )
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/TRP/trcsink.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/TRP/trcsink.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/TRP/trcsink.F90	(revision 13710)
@@ -26,4 +26,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -73,5 +74,5 @@
          iiter(:,:) = 1
       ELSE
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             iiter(ji,jj) = 1
             DO jk = 1, jpkm1
@@ -85,5 +86,5 @@
       ENDIF
 
-      DO_3D_11_11( 1,jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1,jpkm1 )
          IF( tmask(ji,jj,jk) == 1.0 ) THEN
            zwsmax = 0.5 * e3t(ji,jj,jk,Kmm) * rday / rsfact
@@ -145,5 +146,5 @@
       DO jn = 1, 2
          !  first guess of the slopes interior values
-         DO_2D_11_11
+         DO_2D( 1, 1, 1, 1 )
             !
             zstep = rsfact / REAL( kiter(ji,jj), wp ) / 2.
@@ -157,5 +158,5 @@
             ! slopes
             DO jk = 2, jpkm1
-               zign = 0.25 + SIGN( 0.25, ztraz(ji,jj,jk) * ztraz(ji,jj,jk+1) )
+               zign = 0.25 + SIGN( 0.25_wp, ztraz(ji,jj,jk) * ztraz(ji,jj,jk+1) )
                zakz(ji,jj,jk) = ( ztraz(ji,jj,jk) + ztraz(ji,jj,jk+1) ) * zign
             END DO
@@ -163,5 +164,5 @@
             ! Slopes limitation
             DO jk = 2, jpkm1
-               zakz(ji,jj,jk) = SIGN( 1., zakz(ji,jj,jk) ) *        &
+               zakz(ji,jj,jk) = SIGN( 1.0_wp, zakz(ji,jj,jk) ) *        &
                   &             MIN( ABS( zakz(ji,jj,jk) ), 2. * ABS(ztraz(ji,jj,jk+1)), 2. * ABS(ztraz(ji,jj,jk) ) )
             END DO
@@ -185,5 +186,5 @@
       END DO
 
-      DO_3D_11_11( 1,jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1,jpkm1 )
          zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm)
          ztrb(ji,jj,jk) = ztrb(ji,jj,jk) + 2. * zflx
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/TRP/trczdf.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/TRP/trczdf.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/TRP/trczdf.F90	(revision 13710)
@@ -22,5 +22,5 @@
 !!gm
    USE trdtra        ! trends manager: tracers 
-   USE prtctl_trc    ! Print control
+   USE prtctl        ! Print control
 
    IMPLICIT NONE
@@ -69,6 +69,6 @@
       IF( sn_cfctl%l_prttrc )   THEN
          WRITE(charout, FMT="('zdf ')")
-         CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc( tab4d=tr(:,:,:,:,Kaa), mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+         CALL prt_ctl( tab4d_1=tr(:,:,:,:,Kaa), mask1=tmask, clinfo=ctrcnm, clinfo3='trd' )
       END IF
       !
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/TRP/trdmxl_trc.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/TRP/trdmxl_trc.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/TRP/trdmxl_trc.F90	(revision 13710)
@@ -51,4 +51,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -124,5 +125,5 @@
 
             IF( jpktrd_trc < jpk ) THEN                           ! description ???
-               DO_2D_11_11
+               DO_2D( 1, 1, 1, 1 )
                   IF( nmld_trc(ji,jj) <= jpktrd_trc ) THEN
                      zvlmsk(ji,jj) = tmask(ji,jj,1)
@@ -147,5 +148,5 @@
          ! ... Weights for vertical averaging
          wkx_trc(:,:,:) = 0.e0
-         DO_3D_11_11( 1, jpktrd_trc )
+         DO_3D( 1, 1, 1, 1, 1, jpktrd_trc )                       ! initialize wkx_trc with vertical scale factor in mixed-layer
             IF( jk - nmld_trc(ji,jj) < 0 )   wkx_trc(ji,jj,jk) = e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
          END_3D
@@ -258,5 +259,5 @@
          !
          DO jn = 1, jptra
-            DO_2D_11_11
+            DO_2D( 1, 1, 1, 1 )
                ik = nmld_trc(ji,jj)
                IF( ln_trdtrc(jn) )    &
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/TRP/trdmxl_trc_rst.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/TRP/trdmxl_trc_rst.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/TRP/trdmxl_trc_rst.F90	(revision 13710)
@@ -144,22 +144,22 @@
          
          DO jn = 1, jptra
-            CALL iom_get( inum, jpdom_autoglo, 'tmlbb_trc_'  //ctrcnm(jn), tmlbb_trc  (:,:,jn) )
-            CALL iom_get( inum, jpdom_autoglo, 'tmlbn_trc_'  //ctrcnm(jn), tmlbn_trc  (:,:,jn) )
-            CALL iom_get( inum, jpdom_autoglo, 'tmlatfb_trc_'//ctrcnm(jn), tmlatfb_trc(:,:,jn) )
-            CALL iom_get( inum, jpdom_autoglo, 'tmlradb_trc_'//ctrcnm(jn), tmlradb_trc(:,:,jn) )
+            CALL iom_get( inum, jpdom_auto, 'tmlbb_trc_'  //ctrcnm(jn), tmlbb_trc  (:,:,jn) )
+            CALL iom_get( inum, jpdom_auto, 'tmlbn_trc_'  //ctrcnm(jn), tmlbn_trc  (:,:,jn) )
+            CALL iom_get( inum, jpdom_auto, 'tmlatfb_trc_'//ctrcnm(jn), tmlatfb_trc(:,:,jn) )
+            CALL iom_get( inum, jpdom_auto, 'tmlradb_trc_'//ctrcnm(jn), tmlradb_trc(:,:,jn) )
          END DO
          
       ELSE
-         CALL iom_get( inum, jpdom_autoglo, 'rmldbn_trc', rmldbn_trc ) ! needed for rmld_sum
+         CALL iom_get( inum, jpdom_auto, 'rmldbn_trc', rmldbn_trc ) ! needed for rmld_sum
          
          !                                                          ! ===========
          DO jn = 1, jptra                                           ! tracer loop
             !                                                       ! ===========
-            CALL iom_get( inum, jpdom_autoglo, 'tmlatfb_trc_' //ctrcnm(jn), tmlatfb_trc(:,:,jn) )
-            CALL iom_get( inum, jpdom_autoglo, 'tmlbb_trc_'   //ctrcnm(jn), tmlbb_trc  (:,:,jn) )
-            CALL iom_get( inum, jpdom_autoglo, 'tmlradb_trc_' //ctrcnm(jn), tmlradb_trc(:,:,jn) )
-
-            CALL iom_get( inum, jpdom_autoglo, 'tmlbn_trc_'   //ctrcnm(jn), tmlbn_trc   (:,:,jn) ) ! needed for tml_sum
-            CALL iom_get( inum, jpdom_autoglo, 'tml_sumb_trc_'//ctrcnm(jn), tml_sumb_trc(:,:,jn) )
+            CALL iom_get( inum, jpdom_auto, 'tmlatfb_trc_' //ctrcnm(jn), tmlatfb_trc(:,:,jn) )
+            CALL iom_get( inum, jpdom_auto, 'tmlbb_trc_'   //ctrcnm(jn), tmlbb_trc  (:,:,jn) )
+            CALL iom_get( inum, jpdom_auto, 'tmlradb_trc_' //ctrcnm(jn), tmlradb_trc(:,:,jn) )
+
+            CALL iom_get( inum, jpdom_auto, 'tmlbn_trc_'   //ctrcnm(jn), tmlbn_trc   (:,:,jn) ) ! needed for tml_sum
+            CALL iom_get( inum, jpdom_auto, 'tml_sumb_trc_'//ctrcnm(jn), tml_sumb_trc(:,:,jn) )
             
             DO jk = 1, jpltrd_trc
@@ -169,11 +169,11 @@
                   WRITE(charout,FMT="('tmltrd_csum_ub_trc_', A3, '_', I2)") ctrcnm(jn), jk
                ENDIF
-               CALL iom_get( inum, jpdom_autoglo, charout, tmltrd_csum_ub_trc(:,:,jk,jn) )
+               CALL iom_get( inum, jpdom_auto, charout, tmltrd_csum_ub_trc(:,:,jk,jn) )
             END DO
             
-            CALL iom_get( inum, jpdom_autoglo, 'tmltrd_atf_sumb_trc_'//ctrcnm(jn) , &
+            CALL iom_get( inum, jpdom_auto, 'tmltrd_atf_sumb_trc_'//ctrcnm(jn) , &
                  &        tmltrd_atf_sumb_trc(:,:,jn) )
 
-            CALL iom_get( inum, jpdom_autoglo, 'tmltrd_rad_sumb_trc_'//ctrcnm(jn) , &
+            CALL iom_get( inum, jpdom_auto, 'tmltrd_rad_sumb_trc_'//ctrcnm(jn) , &
                  &        tmltrd_rad_sumb_trc(:,:,jn) )
             !                                                       ! ===========
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/TRP/trdtrc.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/TRP/trdtrc.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/TRP/trdtrc.F90	(revision 13710)
@@ -18,4 +18,5 @@
    USE trdmxl_trc        ! Mixed layer trends diag.
    USE iom               ! I/O library
+   USE par_kind
 
    IMPLICIT NONE
@@ -107,4 +108,6 @@
    !!----------------------------------------------------------------------
 
+   USE par_kind
+
    PUBLIC trd_trc
 
@@ -116,5 +119,5 @@
       INTEGER               , INTENT( in )     ::   kjn     ! tracer index
       INTEGER               , INTENT( in )     ::   ktrd    ! tracer trend index
-      REAL, DIMENSION(:,:,:), INTENT( inout )  ::   ptrtrd  ! Temperature or U trend
+      REAL(wp), DIMENSION(:,:,:), INTENT( inout )  ::   ptrtrd  ! Temperature or U trend
       WRITE(*,*) 'trd_trc : You should not have seen this print! error?', ptrtrd(1,1,1)
       WRITE(*,*) '  "      "      : You should not have seen this print! error?', kjn
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/oce_trc.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/oce_trc.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/oce_trc.F90	(revision 13710)
@@ -18,4 +18,27 @@
    USE par_oce , ONLY :   jp_tem   =>   jp_tem     !: indice for temperature
    USE par_oce , ONLY :   jp_sal   =>   jp_sal     !: indice for salinity
+   USE par_oce , ONLY :   nn_hls   =>   nn_hls     !: 
+   USE par_oce , ONLY :   Nis0    =>   Nis0      !: 
+   USE par_oce , ONLY :   Njs0    =>   Njs0      !: 
+   USE par_oce , ONLY :   Nie0    =>   Nie0      !: 
+   USE par_oce , ONLY :   Nje0    =>   Nje0      !: 
+   USE par_oce , ONLY :   Nis1    =>   Nis1      !: 
+   USE par_oce , ONLY :   Njs1    =>   Njs1      !: 
+   USE par_oce , ONLY :   Nie1    =>   Nie1      !: 
+   USE par_oce , ONLY :   Nje1    =>   Nje1      !: 
+   USE par_oce , ONLY :   Nis1nxt2    =>   Nis1nxt2      !: 
+   USE par_oce , ONLY :   Njs1nxt2    =>   Njs1nxt2      !: 
+   USE par_oce , ONLY :   Nie1nxt2    =>   Nie1nxt2      !: 
+   USE par_oce , ONLY :   Nje1nxt2    =>   Nje1nxt2      !: 
+   USE par_oce , ONLY :   Nis2    =>   Nis2      !: 
+   USE par_oce , ONLY :   Njs2    =>   Njs2      !: 
+   USE par_oce , ONLY :   Nie2    =>   Nie2      !: 
+   USE par_oce , ONLY :   Nje2    =>   Nje2      !: 
+   USE par_oce , ONLY :   Ni_0    =>   Ni_0      !: 
+   USE par_oce , ONLY :   Nj_0    =>   Nj_0      !: 
+   USE par_oce , ONLY :   Ni_1    =>   Ni_1      !: 
+   USE par_oce , ONLY :   Nj_1    =>   Nj_1      !: 
+   USE par_oce , ONLY :   Ni_2    =>   Ni_2      !: 
+   USE par_oce , ONLY :   Nj_2    =>   Nj_2      !: 
 
    USE in_out_manager                           !* IO manager *
@@ -62,4 +85,6 @@
    USE traqsr  , ONLY :   rn_abs     =>    rn_abs     !: fraction absorbed in the very near surface
    USE traqsr  , ONLY :   rn_si0     =>    rn_si0     !: very near surface depth of extinction
+   USE traqsr  , ONLY :   nksr       =>    nksr       !: levels below which the light cannot penetrate (depth larger than 391 m)
+   USE traqsr  , ONLY :   rkrgb      =>    rkrgb      !: tabulated attenuation coefficients for RGB absorption
    USE traqsr  , ONLY :   ln_qsr_bio =>    ln_qsr_bio !: flag to use or not the biological fluxes for light
    USE sbcrnf  , ONLY :   rnfmsk     =>    rnfmsk     !: mixed adv scheme in runoffs vicinity (hori.) 
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/prtctl_trc.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/prtctl_trc.F90	(revision 12990)
+++ 	(revision )
@@ -1,295 +1,0 @@
-MODULE prtctl_trc
-   !!======================================================================
-   !!                         ***  MODULE prtctl_trc  ***
-   !! TOP :   print all SUM trends for each processor domain
-   !!======================================================================
-   !! History :    -   !  2005-07  (C. Talandier) original code for OPA
-   !!             1.0  !  2005-10  (C. Ethe     ) adapted to passive tracer
-   !!----------------------------------------------------------------------
-   !!   prt_ctl_trc      :   control print in mpp for passive tracers
-   !!   prt_ctl_trc_info :   ???
-   !!   prt_ctl_trc_init :   ???
-   !!----------------------------------------------------------------------
-   USE par_trc          ! TOP parameters
-   USE oce_trc          ! ocean space and time domain variables
-   USE prtctl           ! print control for OPA
-
-   IMPLICIT NONE
-   PRIVATE
-
-   INTEGER , DIMENSION(:), ALLOCATABLE ::   numid_trc          !: logical unit 
-   INTEGER , DIMENSION(:), ALLOCATABLE ::   nlditl , nldjtl    !: first, last indoor index for each i-domain
-   INTEGER , DIMENSION(:), ALLOCATABLE ::   nleitl , nlejtl    !: first, last indoor index for each j-domain
-   INTEGER , DIMENSION(:), ALLOCATABLE ::   nimpptl, njmpptl   !: i-, j-indexes for each processor
-   INTEGER , DIMENSION(:), ALLOCATABLE ::   nlcitl , nlcjtl    !: dimensions of every subdomain
-   INTEGER , DIMENSION(:), ALLOCATABLE ::   ibonitl, ibonjtl
-
-   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   tra_ctl         !: previous trend values
-
-   PUBLIC prt_ctl_trc         ! called by all subroutines
-   PUBLIC prt_ctl_trc_info    !
-   PUBLIC prt_ctl_trc_init    ! called by opa.F90
-
-CONTAINS
-
-   SUBROUTINE prt_ctl_trc( tab4d, mask, clinfo, ovlap, kdim, clinfo2 )
-      !!----------------------------------------------------------------------
-      !!                     ***  ROUTINE prt_ctl_trc  ***
-      !!
-      !! ** Purpose : - print sum control 3D arrays over the same area 
-      !!                in mono and mpp case. This way can be usefull when
-      !!                debugging a new parametrization in mono or mpp. 
-      !!
-      !! ** Method  : 2 possibilities exist when setting the sn_cfctl%prttrc parameter to
-      !!                .true. in the ocean namelist:
-      !!              - to debug a MPI run .vs. a mono-processor one; 
-      !!                the control print will be done over each sub-domain.
-      !!                The nictl[se] and njctl[se] parameters in the namelist must 
-      !!                be set to zero and [ij]splt to the corresponding splitted
-      !!                domain in MPI along respectively i-, j- directions.
-      !!              - to debug a mono-processor run over the whole domain/a specific area; 
-      !!                in the first case the nictl[se] and njctl[se] parameters must be set
-      !!                to zero else to the indices of the area to be controled. In both cases
-      !!                isplt and jsplt must be set to 1.
-      !!              - All arguments of the above calling sequence are optional so their
-      !!                name must be explicitly typed if used. For instance if the mask
-      !!                array tmask(:,:,:) must be passed through the prt_ctl_trc subroutine, 
-      !!                it must look like: CALL prt_ctl_trc( mask=tmask ).
-      !!----------------------------------------------------------------------
-      REAL(wp)         , DIMENSION(:,:,:,:), INTENT(in), OPTIONAL ::   tab4d     ! 4D array
-      REAL(wp)         , DIMENSION(:,:,:)  , INTENT(in), OPTIONAL ::   mask      ! 3D mask to apply to the tab4d array
-      CHARACTER (len=*), DIMENSION(:)      , INTENT(in), OPTIONAL ::   clinfo    ! information about the tab3d array
-      CHARACTER (len=*)                    , INTENT(in), OPTIONAL ::   clinfo2   ! ???
-      INTEGER                              , INTENT(in), OPTIONAL ::   ovlap     ! overlap value
-      INTEGER                              , INTENT(in), OPTIONAL ::   kdim      ! k- direction for 4D arrays
-      !!
-      INTEGER  ::   overlap, jn, js, sind, eind, kdir, j_id
-      REAL(wp) ::   zsum, zvctl
-      CHARACTER (len=20), ALLOCATABLE, DIMENSION(:) ::   cl
-      CHARACTER (len=10) ::   cl2
-      REAL(wp), DIMENSION(jpi,jpj,jpk)  :: zmask, ztab3d 
-      !!----------------------------------------------------------------------
-
-      ALLOCATE( cl(jptra) )
-      !                                      ! Arrays, scalars initialization 
-      overlap       = 0
-      kdir          = jpkm1
-      zsum          = 0.e0
-      zvctl         = 0.e0
-      cl(:)         = ''
-      cl2           = ''
-      ztab3d(:,:,:) = 0.e0
-      zmask (:,:,:) = 1.e0
-
-      !                                      ! Control of optional arguments
-      IF( PRESENT(ovlap)   )   overlap       = ovlap
-      IF( PRESENT(kdim)    )   kdir          = kdim
-      IF( PRESENT(clinfo ) )   cl(:)         = clinfo(:)
-      IF( PRESENT(clinfo2) )   cl2           = clinfo2
-      IF( PRESENT(mask)    )   zmask (:,:,:) = mask(:,:,:)
-
-      IF( lk_mpp )   THEN      ! processor number
-         sind = narea
-         eind = narea
-      ELSE                     ! processors total number
-         sind = 1
-         eind = ijsplt
-      ENDIF
-
-      ! Loop over each sub-domain, i.e. the total number of processors ijsplt
-      DO js = sind, eind
-         !
-         ! Set logical unit
-         j_id = numid_trc( js - narea + 1 )
-         ! Set indices for the SUM control
-         IF( .NOT. lsp_area ) THEN
-            IF (lk_mpp )   THEN
-               nictls = MAX( 1, nlditl(js) - overlap )
-               nictle = nleitl(js) + overlap * MIN( 1, nlcitl(js) - nleitl(js)) 
-               njctls = MAX( 1, nldjtl(js) - overlap )
-               njctle = nlejtl(js) + overlap * MIN( 1, nlcjtl(js) - nlejtl(js))
-               ! Do not take into account the bound of the domain
-               IF( ibonitl(js) == -1 .OR. ibonitl(js) == 2 )   nictls = MAX( 2, nictls )
-               IF( ibonitl(js) ==  1 .OR. ibonitl(js) == 2 )   nictle = MIN( nictle, nleitl(js) - 1 )
-               IF( ibonjtl(js) == -1 .OR. ibonjtl(js) == 2 )   njctls = MAX( 2, njctls )
-               IF( ibonjtl(js) ==  1 .OR. ibonjtl(js) == 2 )   njctle = MIN( njctle, nlejtl(js) - 1 )
-            ELSE
-               nictls = MAX( 1, nimpptl(js) + nlditl(js) - 1 - overlap )
-               nictle = nimpptl(js) + nleitl(js) - 1 + overlap * MIN( 1, nlcitl(js) - nleitl(js) ) 
-               njctls = MAX( 1, njmpptl(js) + nldjtl(js) - 1 - overlap )
-               njctle = njmpptl(js) + nlejtl(js) - 1 + overlap * MIN( 1, nlcjtl(js) - nlejtl(js) ) 
-               ! Do not take into account the bound of the domain
-               IF( ibonitl(js) == -1 .OR. ibonitl(js) == 2 )   nictls = MAX( 2, nictls )
-               IF( ibonjtl(js) == -1 .OR. ibonjtl(js) == 2 )   njctls = MAX( 2, njctls )
-               IF( ibonitl(js) ==  1 .OR. ibonitl(js) == 2 )   nictle = MIN( nictle, nimpptl(js) + nleitl(js) - 2 )
-               IF( ibonjtl(js) ==  1 .OR. ibonjtl(js) == 2 )   njctle = MIN( njctle, njmpptl(js) + nlejtl(js) - 2 )
-            ENDIF
-         ENDIF
-         !
-         IF( PRESENT(clinfo2) ) THEN
-            DO jn = 1, jptra
-               zvctl  = tra_ctl(jn,js)
-               ztab3d(:,:,:) = tab4d(:,:,:,jn)
-               zsum          = SUM( ztab3d(nictls:nictle,njctls:njctle,1:kdir)   &
-                  &                * zmask(nictls:nictle,njctls:njctle,1:kdir) )
-               WRITE(j_id,FMT="(3x,a,' : ',D23.16)") cl(jn), zsum-zvctl
-               tra_ctl(jn,js) = zsum
-            END DO
-         ELSE
-            DO jn = 1, jptra
-               ztab3d(:,:,:) = tab4d(:,:,:,jn)
-               zsum          = SUM( ztab3d(nictls:nictle,njctls:njctle,1:kdir)   &
-                  &               *  zmask(nictls:nictle,njctls:njctle,1:kdir) )
-               WRITE(j_id,FMT="(3x,a,' : ',D23.16)") cl(jn), zsum
-            END DO
-         ENDIF
-         !
-      END DO
-      !
-      DEALLOCATE( cl )
-      !
-   END SUBROUTINE prt_ctl_trc
-
-
-   SUBROUTINE prt_ctl_trc_info( clinfo )
-      !!----------------------------------------------------------------------
-      !!                     ***  ROUTINE prt_ctl_trc_info  ***
-      !!
-      !! ** Purpose : - print information without any computation
-      !!----------------------------------------------------------------------
-      CHARACTER (len=*), INTENT(in) ::   clinfo      ! information to print
-      !! 
-      INTEGER ::   js, sind, eind, j_id
-      !!----------------------------------------------------------------------
-
-      IF( lk_mpp ) THEN      ! processor number
-         sind = narea
-         eind = narea
-      ELSE                   ! total number of processors
-         sind = 1
-         eind = ijsplt
-      ENDIF
-
-      ! Loop over each sub-domain, i.e. number of processors ijsplt
-      DO js = sind, eind
-         j_id = numid_trc(js - narea + 1)
-         WRITE(j_id,*) clinfo
-      END DO
-      !
-   END SUBROUTINE prt_ctl_trc_info
-
-
-   SUBROUTINE prt_ctl_trc_init
-      !!----------------------------------------------------------------------
-      !!                     ***  ROUTINE prt_ctl_trc_init  ***
-      !!
-      !! ** Purpose :   open ASCII files & compute indices
-      !!----------------------------------------------------------------------
-      INTEGER            ::   js, sind, eind, j_id
-      CHARACTER (len=31) :: clfile_out
-      CHARACTER (len=27) :: clb_name
-      CHARACTER (len=19) :: cl_run
-      !!----------------------------------------------------------------------
-
-      !                             ! Allocate arrays
-      ALLOCATE( nlditl (ijsplt) )
-      ALLOCATE( nldjtl (ijsplt) )
-      ALLOCATE( nleitl (ijsplt) )
-      ALLOCATE( nlejtl (ijsplt) )
-      ALLOCATE( nimpptl(ijsplt) )
-      ALLOCATE( njmpptl(ijsplt) )
-      ALLOCATE( nlcitl (ijsplt) )
-      ALLOCATE( nlcjtl (ijsplt) )
-      ALLOCATE( tra_ctl(jptra,ijsplt) )
-      ALLOCATE( ibonitl(ijsplt) )
-      ALLOCATE( ibonjtl(ijsplt) )
-
-      tra_ctl(:,:) = 0.e0           ! Initialization to zero
-
-      IF( lk_mpp ) THEN
-         sind = narea
-         eind = narea
-         clb_name = "('mpp.top.output_',I4.4)"
-         cl_run = 'MULTI processor run'
-         ! use indices for each area computed by mpp_init subroutine
-         nlditl(1:jpnij) = nldit(:) 
-         nleitl(1:jpnij) = nleit(:) 
-         nldjtl(1:jpnij) = nldjt(:) 
-         nlejtl(1:jpnij) = nlejt(:) 
-         !
-         nimpptl(1:jpnij) = nimppt(:)
-         njmpptl(1:jpnij) = njmppt(:)
-         !
-         nlcitl(1:jpnij) = nlcit(:)
-         nlcjtl(1:jpnij) = nlcjt(:)
-         !
-         ibonitl(1:jpnij) = ibonit(:)
-         ibonjtl(1:jpnij) = ibonjt(:)
-      ELSE
-         sind = 1
-         eind = ijsplt
-         clb_name = "('mono.top.output_',I4.4)"
-         cl_run   = 'MONO processor run '
-         ! compute indices for each area as done in mpp_init subroutine
-         CALL sub_dom
-      ENDIF
-
-      ALLOCATE( numid_trc(eind-sind+1) )
-
-      DO js = sind, eind
-         WRITE(clfile_out,FMT=clb_name) js-1
-         CALL ctl_opn( numid_trc(js -narea + 1), clfile_out, 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
-         j_id = numid_trc(js -narea + 1)
-         WRITE(j_id,*)
-         WRITE(j_id,*) '                 L O D Y C - I P S L'
-         WRITE(j_id,*) '                       N E M 0 '
-         WRITE(j_id,*) '            Ocean General Circulation Model'
-         WRITE(j_id,*) '               version TOP 1.0  (2005) '
-         WRITE(j_id,*)
-         WRITE(j_id,*) '                   PROC number: ', js
-         WRITE(j_id,*)
-         WRITE(j_id,FMT="(19x,a20)") cl_run
-
-         ! Print the SUM control indices
-         IF( .NOT. lsp_area )   THEN
-            IF ( lk_mpp )   THEN
-               nictls = nlditl(js) 
-               nictle = nleitl(js)
-               njctls = nldjtl(js)
-               njctle = nlejtl(js)
-            ELSE
-               nictls = nimpptl(js) + nlditl(js) - 1
-               nictle = nimpptl(js) + nleitl(js) - 1
-               njctls = njmpptl(js) + nldjtl(js) - 1
-               njctle = njmpptl(js) + nlejtl(js) - 1
-            ENDIF
-         ENDIF
-         WRITE(j_id,*) 
-         WRITE(j_id,*) 'prt_tra_ctl :  Sum control indices'
-         WRITE(j_id,*) '~~~~~~~'
-         WRITE(j_id,*)
-         WRITE(j_id,9000)'                                nlej   = ', nlejtl(js), '              '
-         WRITE(j_id,9000)'                  ------------- njctle = ', njctle, ' -------------'
-         WRITE(j_id,9001)'                  |                                       |'
-         WRITE(j_id,9001)'                  |                                       |'
-         WRITE(j_id,9001)'                  |                                       |'
-         WRITE(j_id,9002)'           nictls = ', nictls,  '                           nictle = ', nictle
-         WRITE(j_id,9002)'           nldi   = ', nlditl(js),  '                           nlei   = ', nleitl(js)
-         WRITE(j_id,9001)'                  |                                       |'
-         WRITE(j_id,9001)'                  |                                       |'
-         WRITE(j_id,9001)'                  |                                       |'
-         WRITE(j_id,9004)'  njmpp  = ',njmpptl(js),'   ------------- njctls = ', njctls, ' -------------'
-         WRITE(j_id,9003)'           nimpp  = ', nimpptl(js), '        nldj   = ', nldjtl(js), '              '
-         WRITE(j_id,*)
-         WRITE(j_id,*)
-
-9000     FORMAT(a41,i4.4,a14)
-9001     FORMAT(a59)
-9002     FORMAT(a20,i4.4,a36,i3.3)
-9003     FORMAT(a20,i4.4,a17,i4.4)
-9004     FORMAT(a11,i4.4,a26,i4.4,a14)
-      END DO
-      !
-   END SUBROUTINE prt_ctl_trc_init
-
-END MODULE prtctl_trc
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/trc.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/trc.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/trc.F90	(revision 13710)
@@ -21,5 +21,5 @@
    INTEGER, PUBLIC ::   numonr     = -1   !: reference passive tracer namelist output output.namelist.top
    INTEGER, PUBLIC ::   numstr            !: tracer statistics
-   INTEGER, PUBLIC ::   numrtr            !: trc restart (read )
+   INTEGER, PUBLIC ::   numrtr     = -1   !: trc restart (read )
    INTEGER, PUBLIC ::   numrtw            !: trc restart ( write )
    CHARACTER(:), ALLOCATABLE, PUBLIC ::   numnat_ref   !: character buffer for reference passive tracer namelist_top_ref
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/trcbc.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/trcbc.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/trcbc.F90	(revision 13710)
@@ -48,4 +48,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -151,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
@@ -414,5 +415,5 @@
          ! Remove river dilution for tracers with absent river load
          IF( ln_rnf_ctl .AND. .NOT.ln_trc_cbc(jn) ) THEN
-            DO_2D_01_00
+            DO_2D( 0, 1, 0, 0 )
                DO jk = 1, nk_rnf(ji,jj)
                   zrnf = (rnf(ji,jj) + rnf_b(ji,jj)) * 0.5_wp * r1_rho0 / h_rnf(ji,jj)
@@ -428,5 +429,5 @@
             jl = n_trc_indsbc(jn)
             sf_trcsbc(jl)%fnow(:,:,1) = MAX( rtrn, sf_trcsbc(jl)%fnow(:,:,1) ) ! avoid nedgative value due to interpolation
-            DO_2D_01_00
+            DO_2D( 0, 1, 0, 0 )
                zfact = 1. / ( e3t(ji,jj,1,Kmm) * rn_sbc_time )
                ptr(ji,jj,1,jn,Krhs) = ptr(ji,jj,1,jn,Krhs) + rf_trsfac(jl) * sf_trcsbc(jl)%fnow(ji,jj,1) * zfact
@@ -438,5 +439,5 @@
             IF( l_offline )   rn_rfact = 1._wp
             jl = n_trc_indcbc(jn)
-            DO_2D_01_00
+            DO_2D( 0, 1, 0, 0 )
                DO jk = 1, nk_rnf(ji,jj)
                   zfact = rn_rfact / ( e1e2t(ji,jj) * h_rnf(ji,jj) * rn_cbc_time ) * tmask(ji,jj,1)
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/trcbdy.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/trcbdy.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/trcbdy.F90	(revision 13710)
@@ -49,5 +49,4 @@
       INTEGER                           :: ib_bdy ,ir, jn ,igrd ! Loop indices
       REAL(wp), POINTER, DIMENSION(:,:) ::  ztrc
-      REAL(wp), POINTER                 ::  zfac
       LOGICAL                           :: llrim0               ! indicate if rim 0 is treated
       LOGICAL, DIMENSION(4)             :: llsend1, llrecv1     ! indicate how communications are to be carried out
@@ -61,30 +60,38 @@
          IF( ir == 0 ) THEN   ;   llrim0 = .TRUE.
          ELSE                 ;   llrim0 = .FALSE.
-         END IF
+         ENDIF
          DO ib_bdy=1, nb_bdy
+            !
             DO jn = 1, jptra
                !
-               ztrc => trcdta_bdy(jn,ib_bdy)%trc 
-               zfac => trcdta_bdy(jn,ib_bdy)%rn_fac
+               IF( ASSOCIATED(trcdta_bdy(jn,ib_bdy)%trc) .AND. trcdta_bdy(jn,ib_bdy)%cn_obc /= 'neumann' ) THEN
+                  IF( .NOT. ASSOCIATED(ztrc) )   ALLOCATE( ztrc(idx_bdy(ib_bdy)%nblen(igrd),jpk) )
+                  ztrc(:,:) = trcdta_bdy(jn,ib_bdy)%trc(:,:) * trcdta_bdy(jn,ib_bdy)%rn_fac
+               ENDIF
                !
-               SELECT CASE( TRIM(trcdta_bdy(jn,ib_bdy)%cn_obc) )
+               SELECT CASE( trcdta_bdy(jn,ib_bdy)%cn_obc )
                CASE('none'        )   ;   CYCLE
                CASE('frs'         )   ! treat the whole boundary at once
-                  IF( ir == 0 ) CALL bdy_frs( idx_bdy(ib_bdy),                tr(:,:,:,jn,Krhs), ztrc*zfac )
+                  IF( ir == 0 )           CALL bdy_frs( idx_bdy(ib_bdy),                   tr(:,:,:,jn,Krhs), ztrc )
                CASE('specified'   )   ! treat the whole rim      at once
-                  IF( ir == 0 ) CALL bdy_spe( idx_bdy(ib_bdy),                tr(:,:,:,jn,Krhs), ztrc*zfac )
-               CASE('neumann'     )   ;   CALL bdy_nmn( idx_bdy(ib_bdy), igrd         , tr(:,:,:,jn,Krhs) )   ! tra masked
-               CASE('orlanski'    )   ;   CALL bdy_orl( idx_bdy(ib_bdy), tr(:,:,:,jn,Kbb), tr(:,:,:,jn,Krhs), ztrc*zfac, ll_npo=.false. )
-               CASE('orlanski_npo')   ;   CALL bdy_orl( idx_bdy(ib_bdy), tr(:,:,:,jn,Kbb), tr(:,:,:,jn,Krhs), ztrc*zfac, ll_npo=.true. )
+                  IF( ir == 0 )           CALL bdy_spe( idx_bdy(ib_bdy),                   tr(:,:,:,jn,Krhs), ztrc )
+               CASE('neumann'     )   ;   CALL bdy_nmn( idx_bdy(ib_bdy), igrd            , tr(:,:,:,jn,Krhs),       llrim0 )   ! tra masked
+               CASE('orlanski'    )   ;   CALL bdy_orl( idx_bdy(ib_bdy), tr(:,:,:,jn,Kbb), tr(:,:,:,jn,Krhs), ztrc, llrim0,   &
+                  &                                     ll_npo=.FALSE. )
+               CASE('orlanski_npo')   ;   CALL bdy_orl( idx_bdy(ib_bdy), tr(:,:,:,jn,Kbb), tr(:,:,:,jn,Krhs), ztrc, llrim0,   &
+                  &                                     ll_npo=.TRUE.  )
                CASE DEFAULT           ;   CALL ctl_stop( 'trc_bdy : unrecognised option for open boundaries for passive tracers' )
                END SELECT
                !
             END DO
+            !
+            IF( ASSOCIATED(ztrc) )   DEALLOCATE(ztrc)
+            !
          END DO
          !
          IF( nn_hls > 1 .AND. ir == 1 ) CYCLE   ! at least 2 halos will be corrected -> no need to correct rim 1 before rim 0
-         IF( nn_hls == 1 ) THEN   ;   llsend1(:) = .false.   ;   llrecv1(:) = .false.   ;   END IF
+         IF( nn_hls == 1 ) THEN   ;   llsend1(:) = .false.   ;   llrecv1(:) = .false.   ;   ENDIF
          DO ib_bdy=1, nb_bdy
-            SELECT CASE( TRIM(cn_tra(ib_bdy)) )
+            SELECT CASE( cn_tra(ib_bdy) )
             CASE('neumann')
                llsend1(:) = llsend1(:) .OR. lsend_bdyint(ib_bdy,1,:,ir)   ! possibly every direction, T points
@@ -96,6 +103,6 @@
          END DO
          IF( ANY(llsend1) .OR. ANY(llrecv1) ) THEN   ! if need to send/recv in at least one direction
-            CALL lbc_lnk( 'trcbdy', tr(:,:,:,:,Krhs), 'T',  1., kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 )
-         END IF
+            CALL lbc_lnk( 'trcbdy', tr(:,:,:,:,Krhs), 'T',  1.0_wp, kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 )
+         ENDIF
          !
       END DO   ! ir
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/trcdta.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/trcdta.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/trcdta.F90	(revision 13710)
@@ -41,4 +41,5 @@
    !! Substitutions
 #include "do_loop_substitute.h90"
+#include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -198,5 +199,5 @@
                WRITE(numout,*) 'trc_dta: interpolates passive tracer data onto the s- or mixed s-z-coordinate mesh'
             ENDIF
-            DO_2D_11_11
+            DO_2D( 1, 1, 1, 1 )                 ! vertical interpolation of T & S
                DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points
                   zl = gdept(ji,jj,jk,Kmm)
@@ -206,5 +207,5 @@
                      ztp(jk) = ptrcdta(ji,jj,jpkm1)
                   ELSE                                      ! inbetween : vertical interpolation between jkk & jkk+1
-                     DO jkk = 1, jpkm1                                  ! when  gdept(jkk) < zl < gdept(jkk+1)
+                     DO jkk = 1, jpkm1                                  ! when  gdept_1d(jkk) < zl < gdept_1d(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))
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/trcini.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/trcini.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/trcini.F90	(revision 13710)
@@ -20,5 +20,5 @@
    USE trcnam          ! Namelist read
    USE daymod          ! calendar manager
-   USE prtctl_trc      ! Print control passive tracers (prt_ctl_trc_init routine)
+   USE prtctl          ! Print control passive tracers (prt_ctl_init routine)
    USE trcrst
    USE lib_mpp         ! distribued memory computing library
@@ -31,4 +31,5 @@
    PUBLIC   trc_init   ! called by opa
 
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -93,4 +94,6 @@
       INTEGER             ::  jk, jn  ! dummy loop indices
       CHARACTER (len=25) :: charout
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jptra) :: zzmsk
+      CHARACTER (len=25), DIMENSION(jptra) :: clseb  
       !!----------------------------------------------------------------------
       !
@@ -124,8 +127,13 @@
       IF(lwp) WRITE(numout,*)
       IF(sn_cfctl%l_prttrc) THEN            ! print mean trends (used for debugging)
-         CALL prt_ctl_trc_init
+         CALL prt_ctl_init( 'top', jptra )
          WRITE(charout, FMT="('ini ')")
-         CALL prt_ctl_trc_info( charout )
-         CALL prt_ctl_trc( tab4d=tr(:,:,:,:,Kmm), mask=tmask, clinfo=ctrcnm )
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+         CALL prt_ctl( tab4d_1=tr(:,:,:,:,Kmm), mask1=tmask, clinfo=ctrcnm )
+         DO jn = 1, jptra
+            zzmsk(:,:,:,jn) = tmask(:,:,:)
+            WRITE(clseb(jn),'(a,i2.2)') 'seb ', jn
+         END DO
+         CALL prt_ctl( tab4d_1=zzmsk, mask1=tmask, clinfo=clseb )
       ENDIF
 9000  FORMAT('      tracer nb : ',i2,'      name :',a10,'      initial content :',e18.10)
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/trcrst.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/trcrst.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/trcrst.F90	(revision 13710)
@@ -33,4 +33,5 @@
    PUBLIC   trc_rst_cal
 
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -113,9 +114,9 @@
       ! READ prognostic variables and computes diagnostic variable
       DO jn = 1, jptra
-         CALL iom_get( numrtr, jpdom_autoglo, 'TRN'//ctrcnm(jn), tr(:,:,:,jn,Kmm) )
-      END DO
-
-      DO jn = 1, jptra
-         CALL iom_get( numrtr, jpdom_autoglo, 'TRB'//ctrcnm(jn), tr(:,:,:,jn,Kbb) )
+         CALL iom_get( numrtr, jpdom_auto, 'TRN'//ctrcnm(jn), tr(:,:,:,jn,Kmm) )
+      END DO
+
+      DO jn = 1, jptra
+         CALL iom_get( numrtr, jpdom_auto, 'TRB'//ctrcnm(jn), tr(:,:,:,jn,Kbb) )
       END DO
       !
@@ -236,5 +237,5 @@
                ! calculate start time in hours and minutes
                zdayfrac=adatrj-INT(adatrj)
-               ksecs = NINT(zdayfrac*86400)            ! Nearest second to catch rounding errors in adatrj              
+               ksecs = NINT(zdayfrac*86400)            ! Nearest second to catch rounding errors in adatrj
                ihour = INT(ksecs/3600)
                iminute = ksecs/60-ihour*60
@@ -257,7 +258,8 @@
                adatrj = INT(adatrj)                    ! adatrj set to integer as nn_time0 updated            
              ELSE
+               ndt05 = NINT( 0.5 * rn_Dt  )   !  --- WARNING --- not defined yet are we did not go through day_init
                ! parameters corresponding to nit000 - 1 (as we start the step
                ! loop with a call to day)
-               ndastp = ndate0 - 1       ! ndate0 read in the namelist in dom_nam
+               ndastp = ndate0        ! ndate0 read in the namelist in dom_nam
                nhour   =   nn_time0 / 100
                nminute = ( nn_time0 - nhour * 100 )
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/trcsms.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/trcsms.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/trcsms.F90	(revision 13710)
@@ -20,5 +20,5 @@
    USE trcsms_age         ! AGE
    USE trcsms_my_trc      ! MY_TRC  tracers
-   USE prtctl_trc         ! Print control for debbuging
+   USE prtctl             ! Print control for debbuging
 
    IMPLICIT NONE
@@ -58,6 +58,6 @@
       IF(sn_cfctl%l_prttrc) THEN                       ! print mean trends (used for debugging)
          WRITE(charout, FMT="('sms ')")
-         CALL prt_ctl_trc_info( charout )
-         CALL prt_ctl_trc( tab4d=tr(:,:,:,:,Kmm), mask=tmask, clinfo=ctrcnm )
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+         CALL prt_ctl( tab4d_1=tr(:,:,:,:,Kmm), mask1=tmask, clinfo=ctrcnm )
       ENDIF
       !
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/trcstp.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/trcstp.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/trcstp.F90	(revision 13710)
@@ -22,5 +22,5 @@
    USE sms_pisces,  ONLY : ln_check_mass
    !
-   USE prtctl_trc     ! Print control for debbuging
+   USE prtctl         ! Print control for debbuging
    USE iom            !
    USE in_out_manager !
@@ -37,4 +37,5 @@
    REAL(wp), DIMENSION(:,:,:), SAVE, ALLOCATABLE ::   qsr_arr   ! save qsr during TOP time-step
 
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -91,5 +92,5 @@
       IF(sn_cfctl%l_prttrc) THEN
          WRITE(charout,FMT="('kt =', I4,'  d/m/y =',I2,I2,I4)") kt, nday, nmonth, nyear
-         CALL prt_ctl_trc_info(charout)
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
       ENDIF
       !
@@ -199,5 +200,5 @@
             rsecfst = INT( zkt ) * rn_Dt
             IF(lwp) WRITE(numout,*) 'trc_qsr_mean:   qsr_mean read in the restart file at time-step rsecfst =', rsecfst, ' s '
-            CALL iom_get( numrtr, jpdom_autoglo, 'qsr_mean', qsr_mean )   !  A mean of qsr
+            CALL iom_get( numrtr, jpdom_auto, 'qsr_mean', qsr_mean )   !  A mean of qsr
             CALL iom_get( numrtr, 'nrdcy', zrec )   !  Number of record per days
             IF( INT( zrec ) == nb_rec_per_day ) THEN
@@ -205,8 +206,8 @@
                   IF( jn <= 9 )  THEN
                     WRITE(cl1,'(i1)') jn
-                    CALL iom_get( numrtr, jpdom_autoglo, 'qsr_arr_'//cl1, qsr_arr(:,:,jn) )   !  A mean of qsr
+                    CALL iom_get( numrtr, jpdom_auto, 'qsr_arr_'//cl1, qsr_arr(:,:,jn) )   !  A mean of qsr
                   ELSE
                     WRITE(cl2,'(i2.2)') jn
-                    CALL iom_get( numrtr, jpdom_autoglo, 'qsr_arr_'//cl2, qsr_arr(:,:,jn) )   !  A mean of qsr
+                    CALL iom_get( numrtr, jpdom_auto, 'qsr_arr_'//cl2, qsr_arr(:,:,jn) )   !  A mean of qsr
                   ENDIF
               END DO
Index: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/trcwri.F90
===================================================================
--- NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/trcwri.F90	(revision 12990)
+++ NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/TOP/trcwri.F90	(revision 13710)
@@ -60,7 +60,9 @@
        CALL iom_put( "e3v_0", e3v_0(:,:,:) )
        !
+#if ! defined key_qco
        CALL iom_put( "e3t" , e3t(:,:,:,Kmm) )
        CALL iom_put( "e3u" , e3u(:,:,:,Kmm) )
        CALL iom_put( "e3v" , e3v(:,:,:,Kmm) )
+#endif 
        !
       ENDIF
