Index: NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/PISCES/P4Z/p4zmort.F90
===================================================================
--- NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/PISCES/P4Z/p4zmort.F90	(revision 12236)
+++ NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/PISCES/P4Z/p4zmort.F90	(revision 12340)
@@ -29,4 +29,6 @@
    REAL(wp), PUBLIC ::   mprat2   !:
 
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -75,48 +77,44 @@
       !
       prodcal(:,:,:) = 0._wp   ! calcite production variable set to zero
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zcompaph = MAX( ( tr(ji,jj,jk,jpphy,Kbb) - 1e-8 ), 0.e0 )
-               !     When highly limited by macronutrients, very small cells 
-               !     dominate the community. As a consequence, aggregation
-               !     due to turbulence is negligible. Mortality is also set
-               !     to 0
-               zsizerat = MIN(1., MAX( 0., (quotan(ji,jj,jk) - 0.2) / 0.3) ) * tr(ji,jj,jk,jpphy,Kbb)
-               !     Squared mortality of Phyto similar to a sedimentation term during
-               !     blooms (Doney et al. 1996)
-               zrespp = wchl * 1.e6 * xstep * xdiss(ji,jj,jk) * zcompaph * zsizerat 
-
-               !     Phytoplankton mortality. This mortality loss is slightly
-               !     increased when nutrients are limiting phytoplankton growth
-               !     as observed for instance in case of iron limitation.
-               ztortp = mprat * xstep * zcompaph / ( xkmort + tr(ji,jj,jk,jpphy,Kbb) ) * zsizerat
-
-               zmortp = zrespp + ztortp
-
-               !   Update the arrays TRA which contains the biological sources and sinks
-
-               zfactfe = tr(ji,jj,jk,jpnfe,Kbb)/(tr(ji,jj,jk,jpphy,Kbb)+rtrn)
-               zfactch = tr(ji,jj,jk,jpnch,Kbb)/(tr(ji,jj,jk,jpphy,Kbb)+rtrn)
-               tr(ji,jj,jk,jpphy,Krhs) = tr(ji,jj,jk,jpphy,Krhs) - zmortp
-               tr(ji,jj,jk,jpnch,Krhs) = tr(ji,jj,jk,jpnch,Krhs) - zmortp * zfactch
-               tr(ji,jj,jk,jpnfe,Krhs) = tr(ji,jj,jk,jpnfe,Krhs) - zmortp * zfactfe
-               zprcaca = xfracal(ji,jj,jk) * zmortp
-               !
-               prodcal(ji,jj,jk) = prodcal(ji,jj,jk) + zprcaca  ! prodcal=prodcal(nanophy)+prodcal(microzoo)+prodcal(mesozoo)
-               !
-               zfracal = 0.5 * xfracal(ji,jj,jk)
-               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
-               tr(ji,jj,jk,jpgoc,Krhs) = tr(ji,jj,jk,jpgoc,Krhs) + zfracal * zmortp
-               tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) + ( 1. - zfracal ) * zmortp
-               prodpoc(ji,jj,jk) = prodpoc(ji,jj,jk) + ( 1. - zfracal ) * zmortp
-               prodgoc(ji,jj,jk) = prodgoc(ji,jj,jk) + zfracal * zmortp
-               tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) + ( 1. - zfracal ) * zmortp * zfactfe
-               tr(ji,jj,jk,jpbfe,Krhs) = tr(ji,jj,jk,jpbfe,Krhs) + zfracal * zmortp * zfactfe
-            END DO
-         END DO
-      END DO
+      DO_3D_11_11( 1, jpkm1 )
+         zcompaph = MAX( ( tr(ji,jj,jk,jpphy,Kbb) - 1e-8 ), 0.e0 )
+         !     When highly limited by macronutrients, very small cells 
+         !     dominate the community. As a consequence, aggregation
+         !     due to turbulence is negligible. Mortality is also set
+         !     to 0
+         zsizerat = MIN(1., MAX( 0., (quotan(ji,jj,jk) - 0.2) / 0.3) ) * tr(ji,jj,jk,jpphy,Kbb)
+         !     Squared mortality of Phyto similar to a sedimentation term during
+         !     blooms (Doney et al. 1996)
+         zrespp = wchl * 1.e6 * xstep * xdiss(ji,jj,jk) * zcompaph * zsizerat 
+
+         !     Phytoplankton mortality. This mortality loss is slightly
+         !     increased when nutrients are limiting phytoplankton growth
+         !     as observed for instance in case of iron limitation.
+         ztortp = mprat * xstep * zcompaph / ( xkmort + tr(ji,jj,jk,jpphy,Kbb) ) * zsizerat
+
+         zmortp = zrespp + ztortp
+
+         !   Update the arrays TRA which contains the biological sources and sinks
+
+         zfactfe = tr(ji,jj,jk,jpnfe,Kbb)/(tr(ji,jj,jk,jpphy,Kbb)+rtrn)
+         zfactch = tr(ji,jj,jk,jpnch,Kbb)/(tr(ji,jj,jk,jpphy,Kbb)+rtrn)
+         tr(ji,jj,jk,jpphy,Krhs) = tr(ji,jj,jk,jpphy,Krhs) - zmortp
+         tr(ji,jj,jk,jpnch,Krhs) = tr(ji,jj,jk,jpnch,Krhs) - zmortp * zfactch
+         tr(ji,jj,jk,jpnfe,Krhs) = tr(ji,jj,jk,jpnfe,Krhs) - zmortp * zfactfe
+         zprcaca = xfracal(ji,jj,jk) * zmortp
+         !
+         prodcal(ji,jj,jk) = prodcal(ji,jj,jk) + zprcaca  ! prodcal=prodcal(nanophy)+prodcal(microzoo)+prodcal(mesozoo)
+         !
+         zfracal = 0.5 * xfracal(ji,jj,jk)
+         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
+         tr(ji,jj,jk,jpgoc,Krhs) = tr(ji,jj,jk,jpgoc,Krhs) + zfracal * zmortp
+         tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) + ( 1. - zfracal ) * zmortp
+         prodpoc(ji,jj,jk) = prodpoc(ji,jj,jk) + ( 1. - zfracal ) * zmortp
+         prodgoc(ji,jj,jk) = prodgoc(ji,jj,jk) + zfracal * zmortp
+         tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) + ( 1. - zfracal ) * zmortp * zfactfe
+         tr(ji,jj,jk,jpbfe,Krhs) = tr(ji,jj,jk,jpbfe,Krhs) + zfracal * zmortp * zfactfe
+      END_3D
       !
        IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
@@ -154,45 +152,41 @@
       !     ------------------------------------------------------------
 
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-
-               zcompadi = MAX( ( tr(ji,jj,jk,jpdia,Kbb) - 1e-9), 0. )
-
-               !    Aggregation term for diatoms is increased in case of nutrient
-               !    stress as observed in reality. The stressed cells become more
-               !    sticky and coagulate to sink quickly out of the euphotic zone
-               !     ------------------------------------------------------------
-               !  Phytoplankton respiration 
-               !     ------------------------
-               zlim2   = xlimdia(ji,jj,jk) * xlimdia(ji,jj,jk)
-               zlim1   = 0.25 * ( 1. - zlim2 ) / ( 0.25 + zlim2 ) 
-               zrespp2 = 1.e6 * xstep * (  wchld + wchldm * zlim1 ) * xdiss(ji,jj,jk) * zcompadi * tr(ji,jj,jk,jpdia,Kbb)
-
-               !     Phytoplankton mortality. 
-               !     ------------------------
-               ztortp2 = mprat2 * xstep * tr(ji,jj,jk,jpdia,Kbb)  / ( xkmort + tr(ji,jj,jk,jpdia,Kbb) ) * zcompadi 
-
-               zmortp2 = zrespp2 + ztortp2
-
-               !   Update the arrays tr(:,:,:,:,Krhs) which contains the biological sources and sinks
-               !   ---------------------------------------------------------------------
-               zfactch = tr(ji,jj,jk,jpdch,Kbb) / ( tr(ji,jj,jk,jpdia,Kbb) + rtrn )
-               zfactfe = tr(ji,jj,jk,jpdfe,Kbb) / ( tr(ji,jj,jk,jpdia,Kbb) + rtrn )
-               zfactsi = tr(ji,jj,jk,jpdsi,Kbb) / ( tr(ji,jj,jk,jpdia,Kbb) + rtrn )
-               tr(ji,jj,jk,jpdia,Krhs) = tr(ji,jj,jk,jpdia,Krhs) - zmortp2 
-               tr(ji,jj,jk,jpdch,Krhs) = tr(ji,jj,jk,jpdch,Krhs) - zmortp2 * zfactch
-               tr(ji,jj,jk,jpdfe,Krhs) = tr(ji,jj,jk,jpdfe,Krhs) - zmortp2 * zfactfe
-               tr(ji,jj,jk,jpdsi,Krhs) = tr(ji,jj,jk,jpdsi,Krhs) - zmortp2 * zfactsi
-               tr(ji,jj,jk,jpgsi,Krhs) = tr(ji,jj,jk,jpgsi,Krhs) + zmortp2 * zfactsi
-               tr(ji,jj,jk,jpgoc,Krhs) = tr(ji,jj,jk,jpgoc,Krhs) + zrespp2 + 0.5 * ztortp2
-               tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) + 0.5 * ztortp2
-               prodpoc(ji,jj,jk) = prodpoc(ji,jj,jk) + 0.5 * ztortp2
-               prodgoc(ji,jj,jk) = prodgoc(ji,jj,jk) + zrespp2 + 0.5 * ztortp2
-               tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) + 0.5 * ztortp2 * zfactfe
-               tr(ji,jj,jk,jpbfe,Krhs) = tr(ji,jj,jk,jpbfe,Krhs) + ( zrespp2 + 0.5 * ztortp2 ) * zfactfe
-            END DO
-         END DO
-      END DO
+      DO_3D_11_11( 1, jpkm1 )
+
+         zcompadi = MAX( ( tr(ji,jj,jk,jpdia,Kbb) - 1e-9), 0. )
+
+         !    Aggregation term for diatoms is increased in case of nutrient
+         !    stress as observed in reality. The stressed cells become more
+         !    sticky and coagulate to sink quickly out of the euphotic zone
+         !     ------------------------------------------------------------
+         !  Phytoplankton respiration 
+         !     ------------------------
+         zlim2   = xlimdia(ji,jj,jk) * xlimdia(ji,jj,jk)
+         zlim1   = 0.25 * ( 1. - zlim2 ) / ( 0.25 + zlim2 ) 
+         zrespp2 = 1.e6 * xstep * (  wchld + wchldm * zlim1 ) * xdiss(ji,jj,jk) * zcompadi * tr(ji,jj,jk,jpdia,Kbb)
+
+         !     Phytoplankton mortality. 
+         !     ------------------------
+         ztortp2 = mprat2 * xstep * tr(ji,jj,jk,jpdia,Kbb)  / ( xkmort + tr(ji,jj,jk,jpdia,Kbb) ) * zcompadi 
+
+         zmortp2 = zrespp2 + ztortp2
+
+         !   Update the arrays tr(:,:,:,:,Krhs) which contains the biological sources and sinks
+         !   ---------------------------------------------------------------------
+         zfactch = tr(ji,jj,jk,jpdch,Kbb) / ( tr(ji,jj,jk,jpdia,Kbb) + rtrn )
+         zfactfe = tr(ji,jj,jk,jpdfe,Kbb) / ( tr(ji,jj,jk,jpdia,Kbb) + rtrn )
+         zfactsi = tr(ji,jj,jk,jpdsi,Kbb) / ( tr(ji,jj,jk,jpdia,Kbb) + rtrn )
+         tr(ji,jj,jk,jpdia,Krhs) = tr(ji,jj,jk,jpdia,Krhs) - zmortp2 
+         tr(ji,jj,jk,jpdch,Krhs) = tr(ji,jj,jk,jpdch,Krhs) - zmortp2 * zfactch
+         tr(ji,jj,jk,jpdfe,Krhs) = tr(ji,jj,jk,jpdfe,Krhs) - zmortp2 * zfactfe
+         tr(ji,jj,jk,jpdsi,Krhs) = tr(ji,jj,jk,jpdsi,Krhs) - zmortp2 * zfactsi
+         tr(ji,jj,jk,jpgsi,Krhs) = tr(ji,jj,jk,jpgsi,Krhs) + zmortp2 * zfactsi
+         tr(ji,jj,jk,jpgoc,Krhs) = tr(ji,jj,jk,jpgoc,Krhs) + zrespp2 + 0.5 * ztortp2
+         tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) + 0.5 * ztortp2
+         prodpoc(ji,jj,jk) = prodpoc(ji,jj,jk) + 0.5 * ztortp2
+         prodgoc(ji,jj,jk) = prodgoc(ji,jj,jk) + zrespp2 + 0.5 * ztortp2
+         tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) + 0.5 * ztortp2 * zfactfe
+         tr(ji,jj,jk,jpbfe,Krhs) = tr(ji,jj,jk,jpbfe,Krhs) + ( zrespp2 + 0.5 * ztortp2 ) * zfactfe
+      END_3D
       !
       IF(sn_cfctl%l_prttrc) THEN      ! print mean trends (used for debugging)
