Index: NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/TOP/PISCES/P4Z/p4zpoc.F90
===================================================================
--- NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/TOP/PISCES/P4Z/p4zpoc.F90	(revision 11396)
+++ NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/TOP/PISCES/P4Z/p4zpoc.F90	(revision 13463)
@@ -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
 
@@ -37,4 +37,7 @@
 
 
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -44,5 +47,5 @@
 CONTAINS
 
-   SUBROUTINE p4z_poc( kt, knt )
+   SUBROUTINE p4z_poc( kt, knt, Kbb, Kmm, Krhs )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_poc  ***
@@ -52,5 +55,6 @@
       !! ** Method  : - ???
       !!---------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt, knt   ! ocean time step and ???
+      INTEGER, INTENT(in) ::   kt, knt         ! ocean time step and ???
+      INTEGER, INTENT(in) ::   Kbb, Kmm, Krhs  ! time level indices
       !
       INTEGER  ::   ji, jj, jk, jn
@@ -103,87 +107,83 @@
      ! -----------------------------------------------------------------------
      ztremint(:,:,:) = zremigoc(:,:,:)
-     DO jk = 2, jpkm1
-        DO jj = 1, jpj
-           DO ji = 1, jpi
-              IF (tmask(ji,jj,jk) == 1.) THEN
-                zdep = hmld(ji,jj)
-                !
-                ! In the case of GOC, lability is constant in the mixed layer 
-                ! It is computed only below the mixed layer depth
-                ! ------------------------------------------------------------
-                !
-                IF( gdept_n(ji,jj,jk) > zdep ) THEN
-                  alphat = 0.
-                  remint = 0.
-                  !
-                  zsizek1  = e3t_n(ji,jj,jk-1) / 2. / (wsbio4(ji,jj,jk-1) + rtrn) * tgfunc(ji,jj,jk-1)
-                  zsizek = e3t_n(ji,jj,jk) / 2. / (wsbio4(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk)
-                  !
-                  IF ( gdept_n(ji,jj,jk-1) <= zdep ) THEN
-                    ! 
-                    ! The first level just below the mixed layer needs a 
-                    ! specific treatment because lability is supposed constant
-                    ! everywhere within the mixed layer. This means that 
-                    ! change in lability in the bottom part of the previous cell
-                    ! should not be computed
-                    ! ----------------------------------------------------------
-                    !
-                    ! POC concentration is computed using the lagrangian 
-                    ! framework. It is only used for the lability param
-                    zpoc = trb(ji,jj,jk-1,jpgoc) + consgoc(ji,jj,jk) * rday / rfact2               &
-                    &   * e3t_n(ji,jj,jk) / 2. / (wsbio4(ji,jj,jk) + rtrn)
-                    zpoc = MAX(0., zpoc)
-                    !
-                    DO jn = 1, jcpoc
-                       !
-                       ! Lagrangian based algorithm. The fraction of each 
-                       ! lability class is computed starting from the previous
-                       ! level
-                       ! -----------------------------------------------------
-                       !
-                       ! the concentration of each lability class is calculated
-                       ! as the sum of the different sources and sinks
-                       ! Please note that production of new GOC experiences
-                       ! degradation 
-                       alphag(ji,jj,jk,jn) = alphag(ji,jj,jk-1,jn) * exp( -reminp(jn) * zsizek ) * zpoc &
-                       &   + prodgoc(ji,jj,jk) * alphan(jn) / tgfunc(ji,jj,jk) / reminp(jn)             &
-                       &   * ( 1. - exp( -reminp(jn) * zsizek ) ) * rday / rfact2 
-                       alphat = alphat + alphag(ji,jj,jk,jn)
-                       remint = remint + alphag(ji,jj,jk,jn) * reminp(jn)
-                    END DO
-                  ELSE
-                    !
-                    ! standard algorithm in the rest of the water column
-                    ! See the comments in the previous block.
-                    ! ---------------------------------------------------
-                    !
-                    zpoc = trb(ji,jj,jk-1,jpgoc) + consgoc(ji,jj,jk-1) * rday / rfact2               &
-                    &   * e3t_n(ji,jj,jk-1) / 2. / (wsbio4(ji,jj,jk-1) + rtrn) + consgoc(ji,jj,jk)   &
-                    &   * rday / rfact2 * e3t_n(ji,jj,jk) / 2. / (wsbio4(ji,jj,jk) + rtrn)
-                    zpoc = max(0., zpoc)
-                    !
-                    DO jn = 1, jcpoc
-                       alphag(ji,jj,jk,jn) = alphag(ji,jj,jk-1,jn) * exp( -reminp(jn) * ( zsizek              &
-                       &   + zsizek1 ) ) * zpoc + ( prodgoc(ji,jj,jk-1) / tgfunc(ji,jj,jk-1) * ( 1.           &
-                       &   - exp( -reminp(jn) * zsizek1 ) ) * exp( -reminp(jn) * zsizek ) + prodgoc(ji,jj,jk) &
-                       &   / tgfunc(ji,jj,jk) * ( 1. - exp( -reminp(jn) * zsizek ) ) ) * rday / rfact2 / reminp(jn) * alphan(jn) 
-                       alphat = alphat + alphag(ji,jj,jk,jn)
-                       remint = remint + alphag(ji,jj,jk,jn) * reminp(jn)
-                    END DO
-                  ENDIF
-                  !
-                  DO jn = 1, jcpoc
-                     ! The contribution of each lability class at the current
-                     ! level is computed
-                     alphag(ji,jj,jk,jn) = alphag(ji,jj,jk,jn) / ( alphat + rtrn)
-                  END DO
-                  ! Computation of the mean remineralisation rate
-                  ztremint(ji,jj,jk) =  MAX(0., remint / ( alphat + rtrn) )
-                  !
-                ENDIF
-              ENDIF
+     DO_3D( 1, 1, 1, 1, 2, jpkm1 )
+        IF (tmask(ji,jj,jk) == 1.) THEN
+          zdep = hmld(ji,jj)
+          !
+          ! In the case of GOC, lability is constant in the mixed layer 
+          ! It is computed only below the mixed layer depth
+          ! ------------------------------------------------------------
+          !
+          IF( gdept(ji,jj,jk,Kmm) > zdep ) THEN
+            alphat = 0.
+            remint = 0.
+            !
+            zsizek1  = e3t(ji,jj,jk-1,Kmm) / 2. / (wsbio4(ji,jj,jk-1) + rtrn) * tgfunc(ji,jj,jk-1)
+            zsizek = e3t(ji,jj,jk,Kmm) / 2. / (wsbio4(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk)
+            !
+            IF ( gdept(ji,jj,jk-1,Kmm) <= zdep ) THEN
+              ! 
+              ! The first level just below the mixed layer needs a 
+              ! specific treatment because lability is supposed constant
+              ! everywhere within the mixed layer. This means that 
+              ! change in lability in the bottom part of the previous cell
+              ! should not be computed
+              ! ----------------------------------------------------------
+              !
+              ! POC concentration is computed using the lagrangian 
+              ! framework. It is only used for the lability param
+              zpoc = tr(ji,jj,jk-1,jpgoc,Kbb) + consgoc(ji,jj,jk) * rday / rfact2               &
+              &   * e3t(ji,jj,jk,Kmm) / 2. / (wsbio4(ji,jj,jk) + rtrn)
+              zpoc = MAX(0., zpoc)
+              !
+              DO jn = 1, jcpoc
+                 !
+                 ! Lagrangian based algorithm. The fraction of each 
+                 ! lability class is computed starting from the previous
+                 ! level
+                 ! -----------------------------------------------------
+                 !
+                 ! the concentration of each lability class is calculated
+                 ! as the sum of the different sources and sinks
+                 ! Please note that production of new GOC experiences
+                 ! degradation 
+                 alphag(ji,jj,jk,jn) = alphag(ji,jj,jk-1,jn) * exp( -reminp(jn) * zsizek ) * zpoc &
+                 &   + prodgoc(ji,jj,jk) * alphan(jn) / tgfunc(ji,jj,jk) / reminp(jn)             &
+                 &   * ( 1. - exp( -reminp(jn) * zsizek ) ) * rday / rfact2 
+                 alphat = alphat + alphag(ji,jj,jk,jn)
+                 remint = remint + alphag(ji,jj,jk,jn) * reminp(jn)
+              END DO
+            ELSE
+              !
+              ! standard algorithm in the rest of the water column
+              ! See the comments in the previous block.
+              ! ---------------------------------------------------
+              !
+              zpoc = tr(ji,jj,jk-1,jpgoc,Kbb) + consgoc(ji,jj,jk-1) * rday / rfact2               &
+              &   * e3t(ji,jj,jk-1,Kmm) / 2. / (wsbio4(ji,jj,jk-1) + rtrn) + consgoc(ji,jj,jk)   &
+              &   * rday / rfact2 * e3t(ji,jj,jk,Kmm) / 2. / (wsbio4(ji,jj,jk) + rtrn)
+              zpoc = max(0., zpoc)
+              !
+              DO jn = 1, jcpoc
+                 alphag(ji,jj,jk,jn) = alphag(ji,jj,jk-1,jn) * exp( -reminp(jn) * ( zsizek              &
+                 &   + zsizek1 ) ) * zpoc + ( prodgoc(ji,jj,jk-1) / tgfunc(ji,jj,jk-1) * ( 1.           &
+                 &   - exp( -reminp(jn) * zsizek1 ) ) * exp( -reminp(jn) * zsizek ) + prodgoc(ji,jj,jk) &
+                 &   / tgfunc(ji,jj,jk) * ( 1. - exp( -reminp(jn) * zsizek ) ) ) * rday / rfact2 / reminp(jn) * alphan(jn) 
+                 alphat = alphat + alphag(ji,jj,jk,jn)
+                 remint = remint + alphag(ji,jj,jk,jn) * reminp(jn)
+              END DO
+            ENDIF
+            !
+            DO jn = 1, jcpoc
+               ! The contribution of each lability class at the current
+               ! level is computed
+               alphag(ji,jj,jk,jn) = alphag(ji,jj,jk,jn) / ( alphat + rtrn)
             END DO
-         END DO
-      END DO
+            ! Computation of the mean remineralisation rate
+            ztremint(ji,jj,jk) =  MAX(0., remint / ( alphat + rtrn) )
+            !
+          ENDIF
+        ENDIF
+     END_3D
 
       IF( ln_p4z ) THEN   ;   zremigoc(:,:,:) = MIN( xremip , ztremint(:,:,:) )
@@ -192,64 +192,56 @@
 
       IF( ln_p4z ) THEN
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ! POC disaggregation by turbulence and bacterial activity. 
-                  ! --------------------------------------------------------
-                  zremig = zremigoc(ji,jj,jk) * xstep * tgfunc(ji,jj,jk)
-                  zorem2  = zremig * trb(ji,jj,jk,jpgoc)
-                  orem(ji,jj,jk)      = zorem2
-                  zorem3(ji,jj,jk) = zremig * solgoc * trb(ji,jj,jk,jpgoc)
-                  zofer2 = zremig * trb(ji,jj,jk,jpbfe)
-                  zofer3 = zremig * solgoc * trb(ji,jj,jk,jpbfe)
-
-                  ! -------------------------------------
-                  tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zorem3(ji,jj,jk)
-                  tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) - zorem2 - zorem3(ji,jj,jk)
-                  tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zofer3
-                  tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) - zofer2 - zofer3
-                  tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zorem2
-                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + zofer2
-                  zfolimi(ji,jj,jk)   = zofer2
-               END DO
-            END DO
-         END DO
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+            ! POC disaggregation by turbulence and bacterial activity. 
+            ! --------------------------------------------------------
+            zremig = zremigoc(ji,jj,jk) * xstep * tgfunc(ji,jj,jk)
+            zorem2  = zremig * tr(ji,jj,jk,jpgoc,Kbb)
+            orem(ji,jj,jk)      = zorem2
+            zorem3(ji,jj,jk) = zremig * solgoc * tr(ji,jj,jk,jpgoc,Kbb)
+            zofer2 = zremig * tr(ji,jj,jk,jpbfe,Kbb)
+            zofer3 = zremig * solgoc * tr(ji,jj,jk,jpbfe,Kbb)
+
+            ! -------------------------------------
+            tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) + zorem3(ji,jj,jk)
+            tr(ji,jj,jk,jpgoc,Krhs) = tr(ji,jj,jk,jpgoc,Krhs) - zorem2 - zorem3(ji,jj,jk)
+            tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) + zofer3
+            tr(ji,jj,jk,jpbfe,Krhs) = tr(ji,jj,jk,jpbfe,Krhs) - zofer2 - zofer3
+            tr(ji,jj,jk,jpdoc,Krhs) = tr(ji,jj,jk,jpdoc,Krhs) + zorem2
+            tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + zofer2
+            zfolimi(ji,jj,jk)   = zofer2
+         END_3D
       ELSE
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                   ! POC disaggregation by turbulence and bacterial activity. 
-                  ! --------------------------------------------------------
-                  zremig = zremigoc(ji,jj,jk) * xstep * tgfunc(ji,jj,jk)
-                  zopoc2 = zremig  * trb(ji,jj,jk,jpgoc)
-                  orem(ji,jj,jk) = zopoc2
-                  zorem3(ji,jj,jk) = zremig * solgoc * trb(ji,jj,jk,jpgoc)
-                  zopon2 = xremipn / xremipc * zremig * trb(ji,jj,jk,jpgon)
-                  zopop2 = xremipp / xremipc * zremig * trb(ji,jj,jk,jpgop)
-                  zofer2 = xremipn / xremipc * zremig * trb(ji,jj,jk,jpbfe)
-
-                  ! -------------------------------------
-                  tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zorem3(ji,jj,jk)
-                  tra(ji,jj,jk,jppon) = tra(ji,jj,jk,jppon) + solgoc * zopon2 
-                  tra(ji,jj,jk,jppop) = tra(ji,jj,jk,jppop) + solgoc * zopop2
-                  tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + solgoc * zofer2
-                  tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zopoc2
-                  tra(ji,jj,jk,jpdon) = tra(ji,jj,jk,jpdon) + zopon2
-                  tra(ji,jj,jk,jpdop) = tra(ji,jj,jk,jpdop) + zopop2
-                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + zofer2
-                  tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) - zopoc2 - zorem3(ji,jj,jk)
-                  tra(ji,jj,jk,jpgon) = tra(ji,jj,jk,jpgon) - zopon2 * (1. + solgoc)
-                  tra(ji,jj,jk,jpgop) = tra(ji,jj,jk,jpgop) - zopop2 * (1. + solgoc)
-                  tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) - zofer2 * (1. + solgoc)
-                  zfolimi(ji,jj,jk)   = zofer2
-               END DO
-            END DO
-         END DO
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+             ! POC disaggregation by turbulence and bacterial activity. 
+            ! --------------------------------------------------------
+            zremig = zremigoc(ji,jj,jk) * xstep * tgfunc(ji,jj,jk)
+            zopoc2 = zremig  * tr(ji,jj,jk,jpgoc,Kbb)
+            orem(ji,jj,jk) = zopoc2
+            zorem3(ji,jj,jk) = zremig * solgoc * tr(ji,jj,jk,jpgoc,Kbb)
+            zopon2 = xremipn / xremipc * zremig * tr(ji,jj,jk,jpgon,Kbb)
+            zopop2 = xremipp / xremipc * zremig * tr(ji,jj,jk,jpgop,Kbb)
+            zofer2 = xremipn / xremipc * zremig * tr(ji,jj,jk,jpbfe,Kbb)
+
+            ! -------------------------------------
+            tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) + zorem3(ji,jj,jk)
+            tr(ji,jj,jk,jppon,Krhs) = tr(ji,jj,jk,jppon,Krhs) + solgoc * zopon2 
+            tr(ji,jj,jk,jppop,Krhs) = tr(ji,jj,jk,jppop,Krhs) + solgoc * zopop2
+            tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) + solgoc * zofer2
+            tr(ji,jj,jk,jpdoc,Krhs) = tr(ji,jj,jk,jpdoc,Krhs) + zopoc2
+            tr(ji,jj,jk,jpdon,Krhs) = tr(ji,jj,jk,jpdon,Krhs) + zopon2
+            tr(ji,jj,jk,jpdop,Krhs) = tr(ji,jj,jk,jpdop,Krhs) + zopop2
+            tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + zofer2
+            tr(ji,jj,jk,jpgoc,Krhs) = tr(ji,jj,jk,jpgoc,Krhs) - zopoc2 - zorem3(ji,jj,jk)
+            tr(ji,jj,jk,jpgon,Krhs) = tr(ji,jj,jk,jpgon,Krhs) - zopon2 * (1. + solgoc)
+            tr(ji,jj,jk,jpgop,Krhs) = tr(ji,jj,jk,jpgop,Krhs) - zopop2 * (1. + solgoc)
+            tr(ji,jj,jk,jpbfe,Krhs) = tr(ji,jj,jk,jpbfe,Krhs) - zofer2 * (1. + solgoc)
+            zfolimi(ji,jj,jk)   = zofer2
+         END_3D
       ENDIF
 
-     IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
+     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=tra, mask=tmask, clinfo=ctrcnm)
+        CALL prt_ctl_info( charout, cdcomp = 'top' )
+        CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
      ENDIF
 
@@ -268,18 +260,14 @@
      ! ----------------------------------------------------------------
      ! 
-     DO jk = 1, jpkm1
-        DO jj = 1, jpj
-           DO ji = 1, jpi
-              zdep = hmld(ji,jj)
-              IF (tmask(ji,jj,jk) == 1. .AND. gdept_n(ji,jj,jk) <= zdep ) THEN
-                totprod(ji,jj) = totprod(ji,jj) + prodpoc(ji,jj,jk) * e3t_n(ji,jj,jk) * rday/ rfact2
-                ! The temperature effect is included here
-                totthick(ji,jj) = totthick(ji,jj) + e3t_n(ji,jj,jk)* tgfunc(ji,jj,jk)
-                totcons(ji,jj) = totcons(ji,jj) - conspoc(ji,jj,jk) * e3t_n(ji,jj,jk) * rday/ rfact2    &
-                &                / ( trb(ji,jj,jk,jppoc) + rtrn )
-              ENDIF
-           END DO
-        END DO
-     END DO
+     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
+          totprod(ji,jj) = totprod(ji,jj) + prodpoc(ji,jj,jk) * e3t(ji,jj,jk,Kmm) * rday/ rfact2
+          ! The temperature effect is included here
+          totthick(ji,jj) = totthick(ji,jj) + e3t(ji,jj,jk,Kmm)* tgfunc(ji,jj,jk)
+          totcons(ji,jj) = totcons(ji,jj) - conspoc(ji,jj,jk) * e3t(ji,jj,jk,Kmm) * rday/ rfact2    &
+          &                / ( tr(ji,jj,jk,jppoc,Kbb) + rtrn )
+        ENDIF
+     END_3D
 
      ! Computation of the lability spectrum in the mixed layer. In the mixed 
@@ -287,30 +275,26 @@
      ! ---------------------------------------------------------------------
      ztremint(:,:,:) = zremipoc(:,:,:)
-     DO jk = 1, jpkm1
-        DO jj = 1, jpj
-           DO ji = 1, jpi
-              IF (tmask(ji,jj,jk) == 1.) THEN
-                zdep = hmld(ji,jj)
-                alphat = 0.0
-                remint = 0.0
-                IF( gdept_n(ji,jj,jk) <= zdep ) THEN
-                   DO jn = 1, jcpoc
-                      ! For each lability class, the system is supposed to be 
-                      ! at equilibrium: Prod - Sink - w alphap = 0.
-                      alphap(ji,jj,jk,jn) = totprod(ji,jj) * alphan(jn) / ( reminp(jn)    &
-                      &                     * totthick(ji,jj) + totcons(ji,jj) + wsbio + rtrn )
-                      alphat = alphat + alphap(ji,jj,jk,jn)
-                   END DO
-                   DO jn = 1, jcpoc
-                      alphap(ji,jj,jk,jn) = alphap(ji,jj,jk,jn) / ( alphat + rtrn)
-                      remint = remint + alphap(ji,jj,jk,jn) * reminp(jn)
-                   END DO
-                   ! Mean remineralization rate in the mixed layer
-                   ztremint(ji,jj,jk) =  MAX( 0., remint )
-                ENDIF
-              ENDIF
-           END DO
-        END DO
-     END DO
+     DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+        IF (tmask(ji,jj,jk) == 1.) THEN
+          zdep = hmld(ji,jj)
+          alphat = 0.0
+          remint = 0.0
+          IF( gdept(ji,jj,jk,Kmm) <= zdep ) THEN
+             DO jn = 1, jcpoc
+                ! For each lability class, the system is supposed to be 
+                ! at equilibrium: Prod - Sink - w alphap = 0.
+                alphap(ji,jj,jk,jn) = totprod(ji,jj) * alphan(jn) / ( reminp(jn)    &
+                &                     * totthick(ji,jj) + totcons(ji,jj) + wsbio + rtrn )
+                alphat = alphat + alphap(ji,jj,jk,jn)
+             END DO
+             DO jn = 1, jcpoc
+                alphap(ji,jj,jk,jn) = alphap(ji,jj,jk,jn) / ( alphat + rtrn)
+                remint = remint + alphap(ji,jj,jk,jn) * reminp(jn)
+             END DO
+             ! Mean remineralization rate in the mixed layer
+             ztremint(ji,jj,jk) =  MAX( 0., remint )
+          ENDIF
+        ENDIF
+     END_3D
      !
      IF( ln_p4z ) THEN   ;  zremipoc(:,:,:) = MIN( xremip , ztremint(:,:,:) )
@@ -326,76 +310,72 @@
      ! -----------------------------------------------------------------------
      !
-     DO jk = 2, jpkm1
-        DO jj = 1, jpj
-           DO ji = 1, jpi
-              IF (tmask(ji,jj,jk) == 1.) THEN
-                zdep = hmld(ji,jj)
-                IF( gdept_n(ji,jj,jk) > zdep ) THEN
-                  alphat = 0.
-                  remint = 0.
-                  !
-                  ! the scale factors are corrected with temperature
-                  zsizek1  = e3t_n(ji,jj,jk-1) / 2. / (wsbio3(ji,jj,jk-1) + rtrn) * tgfunc(ji,jj,jk-1)
-                  zsizek = e3t_n(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk)
-                  !
-                  ! Special treatment of the level just below the MXL
-                  ! See the comments in the GOC section
-                  ! ---------------------------------------------------
-                  !
-                  IF ( gdept_n(ji,jj,jk-1) <= zdep ) THEN
-                    !
-                    ! Computation of the POC concentration using the 
-                    ! lagrangian algorithm
-                    zpoc = trb(ji,jj,jk-1,jppoc) + conspoc(ji,jj,jk) * rday / rfact2               &
-                    &   * e3t_n(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn)
-                    zpoc = max(0., zpoc)
-                    ! 
-                    DO jn = 1, jcpoc
-                       ! computation of the lability spectrum applying the 
-                       ! different sources and sinks
-                       alphap(ji,jj,jk,jn) = alphap(ji,jj,jk-1,jn) * exp( -reminp(jn) * zsizek ) * zpoc  &
-                       &   + ( prodpoc(ji,jj,jk) * alphan(jn) + zorem3(ji,jj,jk) * alphag(ji,jj,jk,jn) ) &
-                       &   / tgfunc(ji,jj,jk) / reminp(jn) * rday / rfact2 * ( 1. - exp( -reminp(jn)     &
-                       &   * zsizek ) )
-                       alphap(ji,jj,jk,jn) = MAX( 0., alphap(ji,jj,jk,jn) )
-                       alphat = alphat + alphap(ji,jj,jk,jn)
-                    END DO
-                  ELSE
-                    !
-                    ! Lability parameterization for the interior of the ocean
-                    ! This is very similar to what is done in the previous 
-                    ! block
-                    ! --------------------------------------------------------
-                    !
-                    zpoc = trb(ji,jj,jk-1,jppoc) + conspoc(ji,jj,jk-1) * rday / rfact2               &
-                    &   * e3t_n(ji,jj,jk-1) / 2. / (wsbio3(ji,jj,jk-1) + rtrn) + conspoc(ji,jj,jk)   &
-                    &   * rday / rfact2 * e3t_n(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn)
-                    zpoc = max(0., zpoc)
-                    !
-                    DO jn = 1, jcpoc
-                       alphap(ji,jj,jk,jn) = alphap(ji,jj,jk-1,jn) * exp( -reminp(jn)                       &
-                       &   * ( zsizek + zsizek1 ) ) * zpoc + ( prodpoc(ji,jj,jk-1) * alphan(jn)             & 
-                       &   + zorem3(ji,jj,jk-1) * alphag(ji,jj,jk-1,jn) ) * rday / rfact2 / reminp(jn)      &
-                       &   / tgfunc(ji,jj,jk-1) * ( 1. - exp( -reminp(jn) * zsizek1 ) ) * exp( -reminp(jn)  &
-                       &   * zsizek ) + ( prodpoc(ji,jj,jk) * alphan(jn) + zorem3(ji,jj,jk)                 &
-                       &   * alphag(ji,jj,jk,jn) ) * rday / rfact2 / reminp(jn) / tgfunc(ji,jj,jk) * ( 1.   &
-                       &   - exp( -reminp(jn) * zsizek ) )
-                       alphap(ji,jj,jk,jn) = max(0., alphap(ji,jj,jk,jn) )
-                       alphat = alphat + alphap(ji,jj,jk,jn)
-                    END DO
-                  ENDIF
-                  ! Normalization of the lability spectrum so that the 
-                  ! integral is equal to 1
-                  DO jn = 1, jcpoc
-                     alphap(ji,jj,jk,jn) = alphap(ji,jj,jk,jn) / ( alphat + rtrn)
-                     remint = remint + alphap(ji,jj,jk,jn) * reminp(jn)
-                  END DO
-                  ! Mean remineralization rate in the water column
-                  ztremint(ji,jj,jk) =  MAX( 0., remint )
-                ENDIF
-              ENDIF
+     DO_3D( 1, 1, 1, 1, 2, jpkm1 )
+        IF (tmask(ji,jj,jk) == 1.) THEN
+          zdep = hmld(ji,jj)
+          IF( gdept(ji,jj,jk,Kmm) > zdep ) THEN
+            alphat = 0.
+            remint = 0.
+            !
+            ! the scale factors are corrected with temperature
+            zsizek1  = e3t(ji,jj,jk-1,Kmm) / 2. / (wsbio3(ji,jj,jk-1) + rtrn) * tgfunc(ji,jj,jk-1)
+            zsizek = e3t(ji,jj,jk,Kmm) / 2. / (wsbio3(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk)
+            !
+            ! Special treatment of the level just below the MXL
+            ! See the comments in the GOC section
+            ! ---------------------------------------------------
+            !
+            IF ( gdept(ji,jj,jk-1,Kmm) <= zdep ) THEN
+              !
+              ! Computation of the POC concentration using the 
+              ! lagrangian algorithm
+              zpoc = tr(ji,jj,jk-1,jppoc,Kbb) + conspoc(ji,jj,jk) * rday / rfact2               &
+              &   * e3t(ji,jj,jk,Kmm) / 2. / (wsbio3(ji,jj,jk) + rtrn)
+              zpoc = max(0., zpoc)
+              ! 
+              DO jn = 1, jcpoc
+                 ! computation of the lability spectrum applying the 
+                 ! different sources and sinks
+                 alphap(ji,jj,jk,jn) = alphap(ji,jj,jk-1,jn) * exp( -reminp(jn) * zsizek ) * zpoc  &
+                 &   + ( prodpoc(ji,jj,jk) * alphan(jn) + zorem3(ji,jj,jk) * alphag(ji,jj,jk,jn) ) &
+                 &   / tgfunc(ji,jj,jk) / reminp(jn) * rday / rfact2 * ( 1. - exp( -reminp(jn)     &
+                 &   * zsizek ) )
+                 alphap(ji,jj,jk,jn) = MAX( 0., alphap(ji,jj,jk,jn) )
+                 alphat = alphat + alphap(ji,jj,jk,jn)
+              END DO
+            ELSE
+              !
+              ! Lability parameterization for the interior of the ocean
+              ! This is very similar to what is done in the previous 
+              ! block
+              ! --------------------------------------------------------
+              !
+              zpoc = tr(ji,jj,jk-1,jppoc,Kbb) + conspoc(ji,jj,jk-1) * rday / rfact2               &
+              &   * e3t(ji,jj,jk-1,Kmm) / 2. / (wsbio3(ji,jj,jk-1) + rtrn) + conspoc(ji,jj,jk)   &
+              &   * rday / rfact2 * e3t(ji,jj,jk,Kmm) / 2. / (wsbio3(ji,jj,jk) + rtrn)
+              zpoc = max(0., zpoc)
+              !
+              DO jn = 1, jcpoc
+                 alphap(ji,jj,jk,jn) = alphap(ji,jj,jk-1,jn) * exp( -reminp(jn)                       &
+                 &   * ( zsizek + zsizek1 ) ) * zpoc + ( prodpoc(ji,jj,jk-1) * alphan(jn)             & 
+                 &   + zorem3(ji,jj,jk-1) * alphag(ji,jj,jk-1,jn) ) * rday / rfact2 / reminp(jn)      &
+                 &   / tgfunc(ji,jj,jk-1) * ( 1. - exp( -reminp(jn) * zsizek1 ) ) * exp( -reminp(jn)  &
+                 &   * zsizek ) + ( prodpoc(ji,jj,jk) * alphan(jn) + zorem3(ji,jj,jk)                 &
+                 &   * alphag(ji,jj,jk,jn) ) * rday / rfact2 / reminp(jn) / tgfunc(ji,jj,jk) * ( 1.   &
+                 &   - exp( -reminp(jn) * zsizek ) )
+                 alphap(ji,jj,jk,jn) = max(0., alphap(ji,jj,jk,jn) )
+                 alphat = alphat + alphap(ji,jj,jk,jn)
+              END DO
+            ENDIF
+            ! Normalization of the lability spectrum so that the 
+            ! integral is equal to 1
+            DO jn = 1, jcpoc
+               alphap(ji,jj,jk,jn) = alphap(ji,jj,jk,jn) / ( alphat + rtrn)
+               remint = remint + alphap(ji,jj,jk,jn) * reminp(jn)
             END DO
-         END DO
-      END DO
+            ! Mean remineralization rate in the water column
+            ztremint(ji,jj,jk) =  MAX( 0., remint )
+          ENDIF
+        ENDIF
+     END_3D
 
      IF( ln_p4z ) THEN   ;   zremipoc(:,:,:) = MIN( xremip , ztremint(:,:,:) )
@@ -404,49 +384,41 @@
 
      IF( ln_p4z ) THEN
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  IF (tmask(ji,jj,jk) == 1.) THEN
-                    ! POC disaggregation by turbulence and bacterial activity. 
-                    ! --------------------------------------------------------
-                    zremip          = zremipoc(ji,jj,jk) * xstep * tgfunc(ji,jj,jk)
-                    zorem           = zremip * trb(ji,jj,jk,jppoc)
-                    zofer           = zremip * trb(ji,jj,jk,jpsfe)
-
-                    tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zorem
-                    orem(ji,jj,jk)      = orem(ji,jj,jk) + zorem
-                    tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + zofer
-                    tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) - zorem
-                    tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) - zofer
-                    zfolimi(ji,jj,jk)   = zfolimi(ji,jj,jk) + zofer
-                  ENDIF
-               END DO
-            END DO
-         END DO
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+            IF (tmask(ji,jj,jk) == 1.) THEN
+              ! POC disaggregation by turbulence and bacterial activity. 
+              ! --------------------------------------------------------
+              zremip          = zremipoc(ji,jj,jk) * xstep * tgfunc(ji,jj,jk)
+              zorem           = zremip * tr(ji,jj,jk,jppoc,Kbb)
+              zofer           = zremip * tr(ji,jj,jk,jpsfe,Kbb)
+
+              tr(ji,jj,jk,jpdoc,Krhs) = tr(ji,jj,jk,jpdoc,Krhs) + zorem
+              orem(ji,jj,jk)      = orem(ji,jj,jk) + zorem
+              tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + zofer
+              tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) - zorem
+              tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) - zofer
+              zfolimi(ji,jj,jk)   = zfolimi(ji,jj,jk) + zofer
+            ENDIF
+         END_3D
      ELSE
-       DO jk = 1, jpkm1
-          DO jj = 1, jpj
-             DO ji = 1, jpi
-                ! POC disaggregation by turbulence and bacterial activity. 
-                ! --------------------------------------------------------
-                zremip = zremipoc(ji,jj,jk) * xstep * tgfunc(ji,jj,jk)
-                zopoc  = zremip * trb(ji,jj,jk,jppoc)
-                orem(ji,jj,jk)  = orem(ji,jj,jk) + zopoc
-                zopon  = xremipn / xremipc * zremip * trb(ji,jj,jk,jppon)
-                zopop  = xremipp / xremipc * zremip * trb(ji,jj,jk,jppop)
-                zofer  = xremipn / xremipc * zremip * trb(ji,jj,jk,jpsfe)
-
-                tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) - zopoc
-                tra(ji,jj,jk,jppon) = tra(ji,jj,jk,jppon) - zopon
-                tra(ji,jj,jk,jppop) = tra(ji,jj,jk,jppop) - zopop
-                tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) - zofer
-                tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zopoc
-                tra(ji,jj,jk,jpdon) = tra(ji,jj,jk,jpdon) + zopon 
-                tra(ji,jj,jk,jpdop) = tra(ji,jj,jk,jpdop) + zopop 
-                tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + zofer 
-                zfolimi(ji,jj,jk)   = zfolimi(ji,jj,jk) + zofer
-             END DO
-           END DO
-        END DO
+       DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+          ! POC disaggregation by turbulence and bacterial activity. 
+          ! --------------------------------------------------------
+          zremip = zremipoc(ji,jj,jk) * xstep * tgfunc(ji,jj,jk)
+          zopoc  = zremip * tr(ji,jj,jk,jppoc,Kbb)
+          orem(ji,jj,jk)  = orem(ji,jj,jk) + zopoc
+          zopon  = xremipn / xremipc * zremip * tr(ji,jj,jk,jppon,Kbb)
+          zopop  = xremipp / xremipc * zremip * tr(ji,jj,jk,jppop,Kbb)
+          zofer  = xremipn / xremipc * zremip * tr(ji,jj,jk,jpsfe,Kbb)
+
+          tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) - zopoc
+          tr(ji,jj,jk,jppon,Krhs) = tr(ji,jj,jk,jppon,Krhs) - zopon
+          tr(ji,jj,jk,jppop,Krhs) = tr(ji,jj,jk,jppop,Krhs) - zopop
+          tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) - zofer
+          tr(ji,jj,jk,jpdoc,Krhs) = tr(ji,jj,jk,jpdoc,Krhs) + zopoc
+          tr(ji,jj,jk,jpdon,Krhs) = tr(ji,jj,jk,jpdon,Krhs) + zopon 
+          tr(ji,jj,jk,jpdop,Krhs) = tr(ji,jj,jk,jpdop,Krhs) + zopop 
+          tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + zofer 
+          zfolimi(ji,jj,jk)   = zfolimi(ji,jj,jk) + zofer
+       END_3D
      ENDIF
 
@@ -460,8 +432,8 @@
      ENDIF
 
-      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
+      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=tra, mask=tmask, clinfo=ctrcnm)
+         CALL prt_ctl_info( charout, cdcomp = 'top' )
+         CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
       ENDIF
       !
@@ -497,10 +469,8 @@
       ENDIF
       !
-      REWIND( numnatp_ref )              ! Namelist nampisrem in reference namelist : Pisces remineralization
       READ  ( numnatp_ref, nampispoc, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampispoc in reference namelist', lwp )
-      REWIND( numnatp_cfg )              ! Namelist nampisrem in configuration namelist : Pisces remineralization
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampispoc in reference namelist' )
       READ  ( numnatp_cfg, nampispoc, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampispoc in configuration namelist', lwp )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampispoc in configuration namelist' )
       IF(lwm) WRITE( numonp, nampispoc )
 
