Index: NEMO/trunk/src/TOP/PISCES/P4Z/p4zligand.F90
===================================================================
--- NEMO/trunk/src/TOP/PISCES/P4Z/p4zligand.F90	(revision 12276)
+++ NEMO/trunk/src/TOP/PISCES/P4Z/p4zligand.F90	(revision 12377)
@@ -26,4 +26,6 @@
    REAL(wp), PUBLIC ::  prlgw    !: Photochemical of weak ligand
 
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -33,5 +35,5 @@
 CONTAINS
 
-   SUBROUTINE p4z_ligand( kt, knt )
+   SUBROUTINE p4z_ligand( kt, knt, Kbb, Krhs )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_ligand  ***
@@ -39,5 +41,6 @@
       !! ** Purpose :   Compute remineralization/scavenging of organic ligands
       !!---------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt, knt ! ocean time step
+      INTEGER, INTENT(in) ::   kt, knt   ! ocean time step
+      INTEGER, INTENT(in)  ::  Kbb, Krhs ! time level indices
       !
       INTEGER  ::   ji, jj, jk
@@ -49,28 +52,24 @@
       IF( ln_timing )   CALL timing_start('p4z_ligand')
       !
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               !
-               ! ------------------------------------------------------------------
-               ! Remineralization of iron ligands
-               ! ------------------------------------------------------------------
-               ! production from remineralisation of organic matter
-               zlgwp = orem(ji,jj,jk) * rlig
-               ! decay of weak ligand
-               ! This is based on the idea that as LGW is lower
-               ! there is a larger fraction of refractory OM
-               zlgwr = max( rlgs , rlgw * exp( -2 * (trb(ji,jj,jk,jplgw)*1e9) ) ) ! years
-               zlgwr = 1. / zlgwr * tgfunc(ji,jj,jk) * ( xstep / nyear_len(1) ) * blim(ji,jj,jk) * trb(ji,jj,jk,jplgw)
-               ! photochem loss of weak ligand
-               zlgwpr = prlgw * xstep * etot(ji,jj,jk) * trb(ji,jj,jk,jplgw) * (1. - fr_i(ji,jj))
-               tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) + zlgwp - zlgwr - zlgwpr
-               zligrem(ji,jj,jk)   = zlgwr
-               zligpr(ji,jj,jk)    = zlgwpr
-               zligprod(ji,jj,jk)  = zlgwp
-               !
-            END DO
-         END DO
-      END DO
+      DO_3D_11_11( 1, jpkm1 )
+         !
+         ! ------------------------------------------------------------------
+         ! Remineralization of iron ligands
+         ! ------------------------------------------------------------------
+         ! production from remineralisation of organic matter
+         zlgwp = orem(ji,jj,jk) * rlig
+         ! decay of weak ligand
+         ! This is based on the idea that as LGW is lower
+         ! there is a larger fraction of refractory OM
+         zlgwr = max( rlgs , rlgw * exp( -2 * (tr(ji,jj,jk,jplgw,Kbb)*1e9) ) ) ! years
+         zlgwr = 1. / zlgwr * tgfunc(ji,jj,jk) * ( xstep / nyear_len(1) ) * blim(ji,jj,jk) * tr(ji,jj,jk,jplgw,Kbb)
+         ! photochem loss of weak ligand
+         zlgwpr = prlgw * xstep * etot(ji,jj,jk) * tr(ji,jj,jk,jplgw,Kbb) * (1. - fr_i(ji,jj))
+         tr(ji,jj,jk,jplgw,Krhs) = tr(ji,jj,jk,jplgw,Krhs) + zlgwp - zlgwr - zlgwpr
+         zligrem(ji,jj,jk)   = zlgwr
+         zligpr(ji,jj,jk)    = zlgwpr
+         zligprod(ji,jj,jk) = zlgwp
+         !
+      END_3D
       !
       !  Output of some diagnostics variables
@@ -88,8 +87,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="('ligand1')")
          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
       !
@@ -119,10 +118,6 @@
          WRITE(numout,*) '~~~~~~~~~~~~~~~'
       ENDIF
-
-      REWIND( numnatp_ref )
       READ  ( numnatp_ref, nampislig, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nampislig in reference namelist' )
-
-      REWIND( numnatp_cfg )
       READ  ( numnatp_cfg, nampislig, IOSTAT = ios, ERR = 902 )
 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nampislig in configuration namelist' )
