Index: NEMO/trunk/src/TOP/PISCES/P4Z/p4zmicro.F90
===================================================================
--- NEMO/trunk/src/TOP/PISCES/P4Z/p4zmicro.F90	(revision 12276)
+++ NEMO/trunk/src/TOP/PISCES/P4Z/p4zmicro.F90	(revision 12377)
@@ -42,4 +42,6 @@
    REAL(wp), PUBLIC ::   epshermin   !: minimum growth efficiency for grazing 1
 
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -49,5 +51,5 @@
 CONTAINS
 
-   SUBROUTINE p4z_micro( kt, knt )
+   SUBROUTINE p4z_micro( kt, knt, Kbb, Krhs )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_micro  ***
@@ -59,4 +61,5 @@
       INTEGER, INTENT(in) ::   kt    ! ocean time step
       INTEGER, INTENT(in) ::   knt   ! ??? 
+      INTEGER, INTENT(in) ::   Kbb, Krhs  ! time level indices
       !
       INTEGER  :: ji, jj, jk
@@ -75,115 +78,111 @@
       IF( ln_timing )   CALL timing_start('p4z_micro')
       !
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zcompaz = MAX( ( trb(ji,jj,jk,jpzoo) - 1.e-9 ), 0.e0 )
-               zfact   = xstep * tgfunc2(ji,jj,jk) * zcompaz
-
-               !  Respiration rates of both zooplankton
-               !  -------------------------------------
-               zrespz = resrat * zfact * trb(ji,jj,jk,jpzoo) / ( xkmort + trb(ji,jj,jk,jpzoo) )  &
-                  &   + resrat * zfact * 3. * nitrfac(ji,jj,jk)
-
-               !  Zooplankton mortality. A square function has been selected with
-               !  no real reason except that it seems to be more stable and may mimic predation.
-               !  ---------------------------------------------------------------
-               ztortz = mzrat * 1.e6 * zfact * trb(ji,jj,jk,jpzoo) * (1. - nitrfac(ji,jj,jk))
-
-               zcompadi  = MIN( MAX( ( trb(ji,jj,jk,jpdia) - xthreshdia ), 0.e0 ), xsizedia )
-               zcompaph  = MAX( ( trb(ji,jj,jk,jpphy) - xthreshphy ), 0.e0 )
-               zcompapoc = MAX( ( trb(ji,jj,jk,jppoc) - xthreshpoc ), 0.e0 )
-               
-               !     Microzooplankton grazing
-               !     ------------------------
-               zfood     = xprefn * zcompaph + xprefc * zcompapoc + xprefd * zcompadi
-               zfoodlim  = MAX( 0. , zfood - min(xthresh,0.5*zfood) )
-               zdenom    = zfoodlim / ( xkgraz + zfoodlim )
-               zdenom2   = zdenom / ( zfood + rtrn )
-               zgraze    = grazrat * xstep * tgfunc2(ji,jj,jk) * trb(ji,jj,jk,jpzoo) * (1. - nitrfac(ji,jj,jk))
-
-               zgrazp    = zgraze  * xprefn * zcompaph  * zdenom2 
-               zgrazm    = zgraze  * xprefc * zcompapoc * zdenom2 
-               zgrazsd   = zgraze  * xprefd * zcompadi  * zdenom2 
-
-               zgrazpf   = zgrazp  * trb(ji,jj,jk,jpnfe) / (trb(ji,jj,jk,jpphy) + rtrn)
-               zgrazmf   = zgrazm  * trb(ji,jj,jk,jpsfe) / (trb(ji,jj,jk,jppoc) + rtrn)
-               zgrazsf   = zgrazsd * trb(ji,jj,jk,jpdfe) / (trb(ji,jj,jk,jpdia) + rtrn)
-               !
-               zgraztotc = zgrazp  + zgrazm  + zgrazsd 
-               zgraztotf = zgrazpf + zgrazsf + zgrazmf 
-               zgraztotn = zgrazp * quotan(ji,jj,jk) + zgrazm + zgrazsd * quotad(ji,jj,jk)
-
-               ! Grazing by microzooplankton
-               zgrazing(ji,jj,jk) = zgraztotc
-
-               !    Various remineralization and excretion terms
-               !    --------------------------------------------
-               zgrasrat  = ( zgraztotf + rtrn ) / ( zgraztotc + rtrn )
-               zgrasratn = ( zgraztotn + rtrn ) / ( zgraztotc + rtrn )
-               zepshert  =  MIN( 1., zgrasratn, zgrasrat / ferat3)
-               zbeta     = MAX(0., (epsher - epshermin) )
-               zepsherf  = epshermin + zbeta / ( 1.0 + 0.04E6 * 12. * zfood * zbeta )
-               zepsherv  = zepsherf * zepshert 
-
-               zgrafer   = zgraztotc * MAX( 0. , ( 1. - unass ) * zgrasrat - ferat3 * zepsherv ) 
-               zgrarem   = zgraztotc * ( 1. - zepsherv - unass )
-               zgrapoc   = zgraztotc * unass
-
-               !  Update of the TRA arrays
-               !  ------------------------
-               zgrarsig  = zgrarem * sigma1
-               tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) + zgrarsig
-               tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zgrarsig
-               tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zgrarem - zgrarsig
-               !
-               IF( ln_ligand ) THEN
-                  tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) + (zgrarem - zgrarsig) * ldocz
-                  zzligprod(ji,jj,jk) = (zgrarem - zgrarsig) * ldocz
-               ENDIF
-               !
-               tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) - o2ut * zgrarsig
-               tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + zgrafer
-               zfezoo(ji,jj,jk)    = zgrafer
-               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zgrapoc
-               prodpoc(ji,jj,jk)   = prodpoc(ji,jj,jk) + zgrapoc
-               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zgraztotf * unass
-               tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) + zgrarsig
-               tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * zgrarsig
-               !   Update the arrays TRA which contain the biological sources and sinks
-               !   --------------------------------------------------------------------
-               zmortz = ztortz + zrespz
-               tra(ji,jj,jk,jpzoo) = tra(ji,jj,jk,jpzoo) - zmortz + zepsherv * zgraztotc 
-               tra(ji,jj,jk,jpphy) = tra(ji,jj,jk,jpphy) - zgrazp
-               tra(ji,jj,jk,jpdia) = tra(ji,jj,jk,jpdia) - zgrazsd
-               tra(ji,jj,jk,jpnch) = tra(ji,jj,jk,jpnch) - zgrazp  * trb(ji,jj,jk,jpnch)/(trb(ji,jj,jk,jpphy)+rtrn)
-               tra(ji,jj,jk,jpdch) = tra(ji,jj,jk,jpdch) - zgrazsd * trb(ji,jj,jk,jpdch)/(trb(ji,jj,jk,jpdia)+rtrn)
-               tra(ji,jj,jk,jpdsi) = tra(ji,jj,jk,jpdsi) - zgrazsd * trb(ji,jj,jk,jpdsi)/(trb(ji,jj,jk,jpdia)+rtrn)
-               tra(ji,jj,jk,jpgsi) = tra(ji,jj,jk,jpgsi) + zgrazsd * trb(ji,jj,jk,jpdsi)/(trb(ji,jj,jk,jpdia)+rtrn)
-               tra(ji,jj,jk,jpnfe) = tra(ji,jj,jk,jpnfe) - zgrazpf
-               tra(ji,jj,jk,jpdfe) = tra(ji,jj,jk,jpdfe) - zgrazsf
-               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zmortz - zgrazm
-               prodpoc(ji,jj,jk) = prodpoc(ji,jj,jk) + zmortz
-               conspoc(ji,jj,jk) = conspoc(ji,jj,jk) - zgrazm
-               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + ferat3 * zmortz - zgrazmf
-               !
-               ! calcite production
-               zprcaca = xfracal(ji,jj,jk) * zgrazp
-               prodcal(ji,jj,jk) = prodcal(ji,jj,jk) + zprcaca  ! prodcal=prodcal(nanophy)+prodcal(microzoo)+prodcal(mesozoo)
-               !
-               zprcaca = part * zprcaca
-               tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) - zprcaca
-               tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) - 2. * zprcaca
-               tra(ji,jj,jk,jpcal) = tra(ji,jj,jk,jpcal) + zprcaca
-            END DO
-         END DO
-      END DO
+      DO_3D_11_11( 1, jpkm1 )
+         zcompaz = MAX( ( tr(ji,jj,jk,jpzoo,Kbb) - 1.e-9 ), 0.e0 )
+         zfact   = xstep * tgfunc2(ji,jj,jk) * zcompaz
+
+         !  Respiration rates of both zooplankton
+         !  -------------------------------------
+         zrespz = resrat * zfact * tr(ji,jj,jk,jpzoo,Kbb) / ( xkmort + tr(ji,jj,jk,jpzoo,Kbb) )  &
+            &   + resrat * zfact * 3. * nitrfac(ji,jj,jk)
+
+         !  Zooplankton mortality. A square function has been selected with
+         !  no real reason except that it seems to be more stable and may mimic predation.
+         !  ---------------------------------------------------------------
+         ztortz = mzrat * 1.e6 * zfact * tr(ji,jj,jk,jpzoo,Kbb) * (1. - nitrfac(ji,jj,jk))
+
+         zcompadi  = MIN( MAX( ( tr(ji,jj,jk,jpdia,Kbb) - xthreshdia ), 0.e0 ), xsizedia )
+         zcompaph  = MAX( ( tr(ji,jj,jk,jpphy,Kbb) - xthreshphy ), 0.e0 )
+         zcompapoc = MAX( ( tr(ji,jj,jk,jppoc,Kbb) - xthreshpoc ), 0.e0 )
+         
+         !     Microzooplankton grazing
+         !     ------------------------
+         zfood     = xprefn * zcompaph + xprefc * zcompapoc + xprefd * zcompadi
+         zfoodlim  = MAX( 0. , zfood - min(xthresh,0.5*zfood) )
+         zdenom    = zfoodlim / ( xkgraz + zfoodlim )
+         zdenom2   = zdenom / ( zfood + rtrn )
+         zgraze    = grazrat * xstep * tgfunc2(ji,jj,jk) * tr(ji,jj,jk,jpzoo,Kbb) * (1. - nitrfac(ji,jj,jk))
+
+         zgrazp    = zgraze  * xprefn * zcompaph  * zdenom2 
+         zgrazm    = zgraze  * xprefc * zcompapoc * zdenom2 
+         zgrazsd   = zgraze  * xprefd * zcompadi  * zdenom2 
+
+         zgrazpf   = zgrazp  * tr(ji,jj,jk,jpnfe,Kbb) / (tr(ji,jj,jk,jpphy,Kbb) + rtrn)
+         zgrazmf   = zgrazm  * tr(ji,jj,jk,jpsfe,Kbb) / (tr(ji,jj,jk,jppoc,Kbb) + rtrn)
+         zgrazsf   = zgrazsd * tr(ji,jj,jk,jpdfe,Kbb) / (tr(ji,jj,jk,jpdia,Kbb) + rtrn)
+         !
+         zgraztotc = zgrazp  + zgrazm  + zgrazsd 
+         zgraztotf = zgrazpf + zgrazsf + zgrazmf 
+         zgraztotn = zgrazp * quotan(ji,jj,jk) + zgrazm + zgrazsd * quotad(ji,jj,jk)
+
+         ! Grazing by microzooplankton
+         zgrazing(ji,jj,jk) = zgraztotc
+
+         !    Various remineralization and excretion terms
+         !    --------------------------------------------
+         zgrasrat  = ( zgraztotf + rtrn ) / ( zgraztotc + rtrn )
+         zgrasratn = ( zgraztotn + rtrn ) / ( zgraztotc + rtrn )
+         zepshert  =  MIN( 1., zgrasratn, zgrasrat / ferat3)
+         zbeta     = MAX(0., (epsher - epshermin) )
+         zepsherf  = epshermin + zbeta / ( 1.0 + 0.04E6 * 12. * zfood * zbeta )
+         zepsherv  = zepsherf * zepshert 
+
+         zgrafer   = zgraztotc * MAX( 0. , ( 1. - unass ) * zgrasrat - ferat3 * zepsherv ) 
+         zgrarem   = zgraztotc * ( 1. - zepsherv - unass )
+         zgrapoc   = zgraztotc * unass
+
+         !  Update of the TRA arrays
+         !  ------------------------
+         zgrarsig  = zgrarem * sigma1
+         tr(ji,jj,jk,jppo4,Krhs) = tr(ji,jj,jk,jppo4,Krhs) + zgrarsig
+         tr(ji,jj,jk,jpnh4,Krhs) = tr(ji,jj,jk,jpnh4,Krhs) + zgrarsig
+         tr(ji,jj,jk,jpdoc,Krhs) = tr(ji,jj,jk,jpdoc,Krhs) + zgrarem - zgrarsig
+         !
+         IF( ln_ligand ) THEN
+            tr(ji,jj,jk,jplgw,Krhs) = tr(ji,jj,jk,jplgw,Krhs) + (zgrarem - zgrarsig) * ldocz
+            zzligprod(ji,jj,jk) = (zgrarem - zgrarsig) * ldocz
+         ENDIF
+         !
+         tr(ji,jj,jk,jpoxy,Krhs) = tr(ji,jj,jk,jpoxy,Krhs) - o2ut * zgrarsig
+         tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + zgrafer
+         zfezoo(ji,jj,jk)    = zgrafer
+         tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) + zgrapoc
+         prodpoc(ji,jj,jk)   = prodpoc(ji,jj,jk) + zgrapoc
+         tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) + zgraztotf * unass
+         tr(ji,jj,jk,jpdic,Krhs) = tr(ji,jj,jk,jpdic,Krhs) + zgrarsig
+         tr(ji,jj,jk,jptal,Krhs) = tr(ji,jj,jk,jptal,Krhs) + rno3 * zgrarsig
+         !   Update the arrays TRA which contain the biological sources and sinks
+         !   --------------------------------------------------------------------
+         zmortz = ztortz + zrespz
+         tr(ji,jj,jk,jpzoo,Krhs) = tr(ji,jj,jk,jpzoo,Krhs) - zmortz + zepsherv * zgraztotc 
+         tr(ji,jj,jk,jpphy,Krhs) = tr(ji,jj,jk,jpphy,Krhs) - zgrazp
+         tr(ji,jj,jk,jpdia,Krhs) = tr(ji,jj,jk,jpdia,Krhs) - zgrazsd
+         tr(ji,jj,jk,jpnch,Krhs) = tr(ji,jj,jk,jpnch,Krhs) - zgrazp  * tr(ji,jj,jk,jpnch,Kbb)/(tr(ji,jj,jk,jpphy,Kbb)+rtrn)
+         tr(ji,jj,jk,jpdch,Krhs) = tr(ji,jj,jk,jpdch,Krhs) - zgrazsd * tr(ji,jj,jk,jpdch,Kbb)/(tr(ji,jj,jk,jpdia,Kbb)+rtrn)
+         tr(ji,jj,jk,jpdsi,Krhs) = tr(ji,jj,jk,jpdsi,Krhs) - zgrazsd * tr(ji,jj,jk,jpdsi,Kbb)/(tr(ji,jj,jk,jpdia,Kbb)+rtrn)
+         tr(ji,jj,jk,jpgsi,Krhs) = tr(ji,jj,jk,jpgsi,Krhs) + zgrazsd * tr(ji,jj,jk,jpdsi,Kbb)/(tr(ji,jj,jk,jpdia,Kbb)+rtrn)
+         tr(ji,jj,jk,jpnfe,Krhs) = tr(ji,jj,jk,jpnfe,Krhs) - zgrazpf
+         tr(ji,jj,jk,jpdfe,Krhs) = tr(ji,jj,jk,jpdfe,Krhs) - zgrazsf
+         tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) + zmortz - zgrazm
+         prodpoc(ji,jj,jk) = prodpoc(ji,jj,jk) + zmortz
+         conspoc(ji,jj,jk) = conspoc(ji,jj,jk) - zgrazm
+         tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) + ferat3 * zmortz - zgrazmf
+         !
+         ! calcite production
+         zprcaca = xfracal(ji,jj,jk) * zgrazp
+         prodcal(ji,jj,jk) = prodcal(ji,jj,jk) + zprcaca  ! prodcal=prodcal(nanophy)+prodcal(microzoo)+prodcal(mesozoo)
+         !
+         zprcaca = part * zprcaca
+         tr(ji,jj,jk,jpdic,Krhs) = tr(ji,jj,jk,jpdic,Krhs) - zprcaca
+         tr(ji,jj,jk,jptal,Krhs) = tr(ji,jj,jk,jptal,Krhs) - 2. * zprcaca
+         tr(ji,jj,jk,jpcal,Krhs) = tr(ji,jj,jk,jpcal,Krhs) + zprcaca
+      END_3D
       !
       IF( lk_iomput .AND. knt == nrdttrc ) THEN
-       IF( iom_use("GRAZ1") ) THEN  !   Total grazing of phyto by zooplankton
+        IF( iom_use("GRAZ1") ) THEN  !   Total grazing of phyto by zooplankton
            zgrazing(:,:,jpk) = 0._wp   ; CALL iom_put( "GRAZ1" , zgrazing(:,:,:) * 1.e+3  * rfact2r * tmask(:,:,:) ) 
          ENDIF
          IF( iom_use("FEZOO") ) THEN  
-           zfezoo (:,:,jpk) = 0._wp    ; CALL iom_put( "FEZOO" , zfezoo(:,:,:) * 1e9 * 1.e+3 * rfact2r * tmask(:,:,:) )
+           zfezoo (:,:,jpk) = 0._wp    ; CALL iom_put( "FEZOO", zfezoo(:,:,:) * 1e9 * 1.e+3 * rfact2r * tmask(:,:,:) )
          ENDIF
          IF( ln_ligand ) THEN
@@ -192,8 +191,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="('micro')")
          CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
+         CALL prt_ctl_trc(tab4d=tr(:,:,:,:,Krhs), mask=tmask, clinfo=ctrcnm)
       ENDIF
       !
@@ -228,9 +227,6 @@
       ENDIF
       !
-      REWIND( numnatp_ref )
       READ  ( numnatp_ref, namp4zzoo, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namp4zzoo in reference namelist' )
-
-      REWIND( numnatp_cfg )
       READ  ( numnatp_cfg, namp4zzoo, IOSTAT = ios, ERR = 902 )
 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namp4zzoo in configuration namelist' )
