Index: NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zagg.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zagg.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zagg.F90	(revision 13540)
@@ -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/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zbc.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zbc.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zbc.F90	(revision 13540)
@@ -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/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zbio.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zbio.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zbio.F90	(revision 13540)
@@ -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/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zche.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zche.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zche.F90	(revision 13540)
@@ -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/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zfechem.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zfechem.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zfechem.F90	(revision 13540)
@@ -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/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zflx.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zflx.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zflx.F90	(revision 13540)
@@ -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/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zligand.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zligand.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zligand.F90	(revision 13540)
@@ -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/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zlim.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zlim.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zlim.F90	(revision 13540)
@@ -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/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zlys.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zlys.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zlys.F90	(revision 13540)
@@ -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/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zmeso.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zmeso.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zmeso.F90	(revision 13540)
@@ -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/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zmicro.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zmicro.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zmicro.F90	(revision 13540)
@@ -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/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zmort.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zmort.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zmort.F90	(revision 13540)
@@ -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/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zopt.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zopt.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zopt.F90	(revision 13540)
@@ -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/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zpoc.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zpoc.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zpoc.F90	(revision 13540)
@@ -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/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zprod.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zprod.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zprod.F90	(revision 13540)
@@ -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/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zrem.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zrem.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zrem.F90	(revision 13540)
@@ -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/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zsed.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zsed.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zsed.F90	(revision 13540)
@@ -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
@@ -314,6 +315,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/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zsink.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zsink.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zsink.F90	(revision 13540)
@@ -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/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zsms.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zsms.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p4zsms.F90	(revision 13540)
@@ -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/r12377_ticket2386/src/TOP/PISCES/P4Z/p5zlim.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p5zlim.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p5zlim.F90	(revision 13540)
@@ -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/r12377_ticket2386/src/TOP/PISCES/P4Z/p5zmeso.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p5zmeso.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p5zmeso.F90	(revision 13540)
@@ -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/r12377_ticket2386/src/TOP/PISCES/P4Z/p5zmicro.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p5zmicro.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p5zmicro.F90	(revision 13540)
@@ -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/r12377_ticket2386/src/TOP/PISCES/P4Z/p5zmort.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p5zmort.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p5zmort.F90	(revision 13540)
@@ -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/r12377_ticket2386/src/TOP/PISCES/P4Z/p5zprod.F90
===================================================================
--- NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p5zprod.F90	(revision 12511)
+++ NEMO/branches/2020/r12377_ticket2386/src/TOP/PISCES/P4Z/p5zprod.F90	(revision 13540)
@@ -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
       !
