Index: NEMO/trunk/src/TOP/PISCES/P4Z/p4zsink.F90
===================================================================
--- NEMO/trunk/src/TOP/PISCES/P4Z/p4zsink.F90	(revision 12276)
+++ NEMO/trunk/src/TOP/PISCES/P4Z/p4zsink.F90	(revision 12377)
@@ -38,4 +38,6 @@
    INTEGER  :: ik100
 
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -49,5 +51,5 @@
    !!----------------------------------------------------------------------
 
-   SUBROUTINE p4z_sink ( kt, knt )
+   SUBROUTINE p4z_sink ( kt, knt, Kbb, Kmm, Krhs )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_sink  ***
@@ -59,4 +61,5 @@
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) :: kt, knt
+      INTEGER, INTENT(in) :: Kbb, Kmm, Krhs  ! time level indices
       INTEGER  ::   ji, jj, jk
       CHARACTER (len=25) :: charout
@@ -77,13 +80,9 @@
       !    by data and from the coagulation theory
       !    -----------------------------------------------------------
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1,jpi
-               zmax  = MAX( heup_01(ji,jj), hmld(ji,jj) )
-               zfact = MAX( 0., gdepw_n(ji,jj,jk+1) - zmax ) / wsbio2scale
-               wsbio4(ji,jj,jk) = wsbio2 + MAX(0., ( wsbio2max - wsbio2 )) * zfact
-            END DO
-         END DO
-      END DO
+      DO_3D_11_11( 1, jpkm1 )
+         zmax  = MAX( heup_01(ji,jj), hmld(ji,jj) )
+         zfact = MAX( 0., gdepw(ji,jj,jk+1,Kmm) - zmax ) / wsbio2scale
+         wsbio4(ji,jj,jk) = wsbio2 + MAX(0., ( wsbio2max - wsbio2 )) * zfact
+      END_3D
 
       ! limit the values of the sinking speeds to avoid numerical instabilities  
@@ -102,10 +101,10 @@
       !   Compute the sedimentation term using p4zsink2 for all the sinking particles
       !   -----------------------------------------------------
-      CALL trc_sink( kt, wsbio3, sinking , jppoc, rfact2 )
-      CALL trc_sink( kt, wsbio3, sinkfer , jpsfe, rfact2 )
-      CALL trc_sink( kt, wsbio4, sinking2, jpgoc, rfact2 )
-      CALL trc_sink( kt, wsbio4, sinkfer2, jpbfe, rfact2 )
-      CALL trc_sink( kt, wsbio4, sinksil , jpgsi, rfact2 )
-      CALL trc_sink( kt, wsbio4, sinkcal , jpcal, rfact2 )
+      CALL trc_sink( kt, Kbb, Kmm, wsbio3, sinking , jppoc, rfact2 )
+      CALL trc_sink( kt, Kbb, Kmm, wsbio3, sinkfer , jpsfe, rfact2 )
+      CALL trc_sink( kt, Kbb, Kmm, wsbio4, sinking2, jpgoc, rfact2 )
+      CALL trc_sink( kt, Kbb, Kmm, wsbio4, sinkfer2, jpbfe, rfact2 )
+      CALL trc_sink( kt, Kbb, Kmm, wsbio4, sinksil , jpgsi, rfact2 )
+      CALL trc_sink( kt, Kbb, Kmm, wsbio4, sinkcal , jpcal, rfact2 )
 
       IF( ln_p5z ) THEN
@@ -117,8 +116,8 @@
          !   Compute the sedimentation term using p4zsink2 for all the sinking particles
          !   -----------------------------------------------------
-         CALL trc_sink( kt, wsbio3, sinkingn , jppon, rfact2 )
-         CALL trc_sink( kt, wsbio3, sinkingp , jppop, rfact2 )
-         CALL trc_sink( kt, wsbio4, sinking2n, jpgon, rfact2 )
-         CALL trc_sink( kt, wsbio4, sinking2p, jpgop, rfact2 )
+         CALL trc_sink( kt, Kbb, Kmm, wsbio3, sinkingn , jppon, rfact2 )
+         CALL trc_sink( kt, Kbb, Kmm, wsbio3, sinkingp , jppop, rfact2 )
+         CALL trc_sink( kt, Kbb, Kmm, wsbio4, sinking2n, jpgon, rfact2 )
+         CALL trc_sink( kt, Kbb, Kmm, wsbio4, sinking2p, jpgop, rfact2 )
       ENDIF
 
@@ -142,8 +141,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="('sink')")
          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
       !
