Index: NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/PISCES/P4Z/p5zmeso.F90
===================================================================
--- NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/PISCES/P4Z/p5zmeso.F90	(revision 12260)
+++ NEMO/branches/2019/dev_r11943_MERGE_2019/src/TOP/PISCES/P4Z/p5zmeso.F90	(revision 12340)
@@ -51,4 +51,6 @@
    LOGICAL,  PUBLIC ::  bmetexc2     !: Use of excess carbon for respiration
 
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -96,253 +98,249 @@
       IF ( bmetexc2 ) zmetexcess = 1.0
 
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zcompam   = MAX( ( tr(ji,jj,jk,jpmes,Kbb) - 1.e-9 ), 0.e0 )
-               zfact     = xstep * tgfunc2(ji,jj,jk) * zcompam
-
-               !   Michaelis-Menten mortality rates of mesozooplankton
-               !   ---------------------------------------------------
-               zrespz   = resrat2 * zfact * ( tr(ji,jj,jk,jpmes,Kbb) / ( xkmort + tr(ji,jj,jk,jpmes,Kbb) )  &
-               &          + 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   = mzrat2 * 1.e6 * zfact * tr(ji,jj,jk,jpmes,Kbb) * (1. - nitrfac(ji,jj,jk))
-
-               !   Computation of the abundance of the preys
-               !   A threshold can be specified in the namelist
-               !   --------------------------------------------
-               zcompadi  = MAX( ( tr(ji,jj,jk,jpdia,Kbb) - xthresh2dia ), 0.e0 )
-               zcompaz   = MAX( ( tr(ji,jj,jk,jpzoo,Kbb) - xthresh2zoo ), 0.e0 )
-               zcompaph  = MAX( ( tr(ji,jj,jk,jpphy,Kbb) - xthresh2phy ), 0.e0 )
-               zcompapoc = MAX( ( tr(ji,jj,jk,jppoc,Kbb) - xthresh2poc ), 0.e0 )
-               zcompames = MAX( ( tr(ji,jj,jk,jpmes,Kbb) - xthresh2mes ), 0.e0 )
-
-               !   Mesozooplankton grazing
-               !   ------------------------
-               zfood     = xpref2d * zcompadi + xpref2z * zcompaz + xpref2n * zcompaph + xpref2c * zcompapoc   &
-               &           + xpref2m * zcompames 
-               zfoodlim  = MAX( 0., zfood - MIN( 0.5 * zfood, xthresh2 ) )
-               zdenom    = zfoodlim / ( xkgraz2 + zfoodlim )
-               zgraze2   = grazrat2 * xstep * tgfunc2(ji,jj,jk) * tr(ji,jj,jk,jpmes,Kbb) * (1. - nitrfac(ji,jj,jk)) 
-
-               !   An active switching parameterization is used here.
-               !   We don't use the KTW parameterization proposed by 
-               !   Vallina et al. because it tends to produce to steady biomass
-               !   composition and the variance of Chl is too low as it grazes
-               !   too strongly on winning organisms. Thus, instead of a square
-               !   a 1.5 power value is used which decreases the pressure on the
-               !   most abundant species
-               !   ------------------------------------------------------------  
-               ztmp1 = xpref2n * zcompaph**1.5
-               ztmp2 = xpref2m * zcompames**1.5
-               ztmp3 = xpref2c * zcompapoc**1.5
-               ztmp4 = xpref2d * zcompadi**1.5
-               ztmp5 = xpref2z * zcompaz**1.5
-               ztmptot = ztmp1 + ztmp2 + ztmp3 + ztmp4 + ztmp5 + rtrn
-               ztmp1 = ztmp1 / ztmptot
-               ztmp2 = ztmp2 / ztmptot
-               ztmp3 = ztmp3 / ztmptot
-               ztmp4 = ztmp4 / ztmptot
-               ztmp5 = ztmp5 / ztmptot
-
-               !   Mesozooplankton regular grazing on the different preys
-               !   ------------------------------------------------------
-               zgrazdc   = zgraze2 * ztmp4 * zdenom
-               zgrazdn   = zgrazdc * tr(ji,jj,jk,jpndi,Kbb) / ( tr(ji,jj,jk,jpdia,Kbb) + rtrn)
-               zgrazdp   = zgrazdc * tr(ji,jj,jk,jppdi,Kbb) / ( tr(ji,jj,jk,jpdia,Kbb) + rtrn)
-               zgrazdf   = zgrazdc * tr(ji,jj,jk,jpdfe,Kbb) / ( tr(ji,jj,jk,jpdia,Kbb) + rtrn)
-               zgrazz    = zgraze2 * ztmp5 * zdenom
-               zgrazm    = zgraze2 * ztmp2 * zdenom
-               zgraznc   = zgraze2 * ztmp1 * zdenom
-               zgraznn   = zgraznc * tr(ji,jj,jk,jpnph,Kbb) / ( tr(ji,jj,jk,jpphy,Kbb) + rtrn)
-               zgraznp   = zgraznc * tr(ji,jj,jk,jppph,Kbb) / ( tr(ji,jj,jk,jpphy,Kbb) + rtrn)
-               zgraznf   = zgraznc * tr(ji,jj,jk,jpnfe,Kbb) / ( tr(ji,jj,jk,jpphy,Kbb) + rtrn)
-               zgrazpoc  = zgraze2 * ztmp3 * zdenom
-               zgrazpon  = zgrazpoc * tr(ji,jj,jk,jppon,Kbb) / ( tr(ji,jj,jk,jppoc,Kbb) + rtrn)
-               zgrazpop  = zgrazpoc * tr(ji,jj,jk,jppop,Kbb) / ( tr(ji,jj,jk,jppoc,Kbb) + rtrn)
-               zgrazpof  = zgrazpoc * tr(ji,jj,jk,jpsfe,Kbb) / ( tr(ji,jj,jk,jppoc,Kbb) + rtrn)
-
-               !   Mesozooplankton flux feeding on GOC
-               !   ----------------------------------
-               zgrazffeg = grazflux  * xstep * wsbio4(ji,jj,jk)      &
-               &           * tgfunc2(ji,jj,jk) * tr(ji,jj,jk,jpgoc,Kbb) * tr(ji,jj,jk,jpmes,Kbb)  &
-               &           * (1. - nitrfac(ji,jj,jk))
-               zgrazfffg = zgrazffeg * tr(ji,jj,jk,jpbfe,Kbb) / (tr(ji,jj,jk,jpgoc,Kbb) + rtrn)
-               zgrazffng = zgrazffeg * tr(ji,jj,jk,jpgon,Kbb) / (tr(ji,jj,jk,jpgoc,Kbb) + rtrn)
-               zgrazffpg = zgrazffeg * tr(ji,jj,jk,jpgop,Kbb) / (tr(ji,jj,jk,jpgoc,Kbb) + rtrn)
-               zgrazffep = grazflux  * xstep *  wsbio3(ji,jj,jk)     &
-               &           * tgfunc2(ji,jj,jk) * tr(ji,jj,jk,jppoc,Kbb) * tr(ji,jj,jk,jpmes,Kbb)   &
-               &           * (1. - nitrfac(ji,jj,jk))
-               zgrazfffp = zgrazffep * tr(ji,jj,jk,jpsfe,Kbb) / (tr(ji,jj,jk,jppoc,Kbb) + rtrn)
-               zgrazffnp = zgrazffep * tr(ji,jj,jk,jppon,Kbb) / (tr(ji,jj,jk,jppoc,Kbb) + rtrn)
-               zgrazffpp = zgrazffep * tr(ji,jj,jk,jppop,Kbb) / (tr(ji,jj,jk,jppoc,Kbb) + rtrn)
-               !
-               zgraztotc  = zgrazdc + zgrazz + zgraznc + zgrazm + zgrazpoc + zgrazffep + zgrazffeg
-
-               !   Compute the proportion of filter feeders
-               !   ----------------------------------------  
-               zproport  = (zgrazffep + zgrazffeg)/(rtrn + zgraztotc)
-
-               !   Compute fractionation of aggregates. It is assumed that 
-               !   diatoms based aggregates are more prone to fractionation
-               !   since they are more porous (marine snow instead of fecal pellets)
-               !   ----------------------------------------------------------------
-               zratio    = tr(ji,jj,jk,jpgsi,Kbb) / ( tr(ji,jj,jk,jpgoc,Kbb) + rtrn )
-               zratio2   = zratio * zratio
-               zfracc    = zproport * grazflux  * xstep * wsbio4(ji,jj,jk)      &
-               &          * tr(ji,jj,jk,jpgoc,Kbb) * tr(ji,jj,jk,jpmes,Kbb)          &
-               &          * ( 0.2 + 3.8 * zratio2 / ( 1.**2 + zratio2 ) )
-               zfracfe   = zfracc * tr(ji,jj,jk,jpbfe,Kbb) / (tr(ji,jj,jk,jpgoc,Kbb) + rtrn)
-               zfracn    = zfracc * tr(ji,jj,jk,jpgon,Kbb) / (tr(ji,jj,jk,jpgoc,Kbb) + rtrn)
-               zfracp    = zfracc * tr(ji,jj,jk,jpgop,Kbb) / (tr(ji,jj,jk,jpgoc,Kbb) + rtrn)
-
-               zgrazffep = zproport * zgrazffep   ;   zgrazffeg = zproport * zgrazffeg
-               zgrazfffp = zproport * zgrazfffp   ;   zgrazfffg = zproport * zgrazfffg
-               zgrazffnp = zproport * zgrazffnp   ;   zgrazffng = zproport * zgrazffng
-               zgrazffpp = zproport * zgrazffpp   ;   zgrazffpg = zproport * zgrazffpg
-
-               zgraztotc  = zgrazdc + zgrazz + zgraznc + zgrazm + zgrazpoc + zgrazffep + zgrazffeg
-               zgraztotf  = zgrazdf + zgraznf + ( zgrazz + zgrazm ) * ferat3 + zgrazpof &
-               &            + zgrazfffp + zgrazfffg
-               zgraztotn  = zgrazdn + (zgrazm + zgrazz) * no3rat3 + zgraznn + zgrazpon  &
-               &            + zgrazffnp + zgrazffng
-               zgraztotp  = zgrazdp + (zgrazz + zgrazm) * po4rat3 + zgraznp + zgrazpop  &
-               &            + zgrazffpp + zgrazffpg
-
-
-               ! Total grazing ( grazing by microzoo is already computed in p5zmicro )
-               zgrazing(ji,jj,jk) = zgraztotc
-
-               !   Stoichiometruc ratios of the food ingested by zooplanton 
-               !   --------------------------------------------------------
-               zgrasratf  =  (zgraztotf + rtrn) / ( zgraztotc + rtrn )
-               zgrasratn  =  (zgraztotn + rtrn) / ( zgraztotc + rtrn )
-               zgrasratp  =  (zgraztotp + rtrn) / ( zgraztotc + rtrn )
-
-               !   Growth efficiency is made a function of the quality 
-               !   and the quantity of the preys
-               !   ---------------------------------------------------
-               zepshert  = MIN( 1., zgrasratn/ no3rat3, zgrasratp/ po4rat3, zgrasratf / ferat3)
-               zbeta     = MAX(0., (epsher2 - epsher2min) )
-               zepsherf  = epsher2min + zbeta / ( 1.0 + 0.04E6 * 12. * zfood * zbeta )
-               zepsherv  = zepsherf * zepshert
-
-               !   Respiration of mesozooplankton
-               !   Excess carbon in the food is used preferentially
-               !   ----------------  ------------------------------
-               zexcess  = zgraztotc * zepsherf * (1.0 - zepshert) * zmetexcess 
-               zbasresb = MAX(0., zrespz - zexcess)
-               zbasresi = zexcess + MIN(0., zrespz - zexcess)
-               zrespirc = srespir2 * zepsherv * zgraztotc + zbasresb
-
-               !   When excess carbon is used, the other elements in excess
-               !   are also used proportionally to their abundance
-               !   --------------------------------------------------------
-               zexcess  = ( zgrasratn/ no3rat3 - zepshert ) / ( 1.0 - zepshert + rtrn)
-               zbasresn = zbasresi * zexcess * zgrasratn
-               zexcess  = ( zgrasratp/ po4rat3 - zepshert ) / ( 1.0 - zepshert + rtrn)
-               zbasresp = zbasresi * zexcess * zgrasratp
-               zexcess  = ( zgrasratf/ ferat3 - zepshert ) / ( 1.0 - zepshert + rtrn)
-               zbasresf = zbasresi * zexcess * zgrasratf
-
-               !   Voiding of the excessive elements as organic matter
-               !   --------------------------------------------------------
-               zgradoct = (1. - unass2c - zepsherv) * zgraztotc - zbasresi
-               zgradont = (1. - unass2n) * zgraztotn - zepsherv * no3rat3 * zgraztotc - zbasresn
-               zgradopt = (1. - unass2p) * zgraztotp - zepsherv * po4rat3 * zgraztotc - zbasresp
-               zgrareft = (1. - unass2c) * zgraztotf - zepsherv * ferat3 * zgraztotc - zbasresf
-               ztmp1   = ( 1. - epsher2 - unass2c ) /( 1. - 0.8 * epsher2 ) * ztortz
-               zgradoc = (zgradoct + ztmp1) * ssigma2
-               zgradon = (zgradont + no3rat3 * ztmp1) * ssigma2
-               zgradop = (zgradopt + po4rat3 * ztmp1) * ssigma2
-               zgratmp = 0.2 * epsher2 /( 1. - 0.8 * epsher2 ) * ztortz
-
-               !  Since only semilabile DOM is represented in PISCES
-               !  part of DOM is in fact labile and is then released
-               !  as dissolved inorganic compounds (ssigma2)
-               !  --------------------------------------------------
-               zgrarem = zgratmp + ( zgradoct + ztmp1 ) * (1.0 - ssigma2)
-               zgraren = no3rat3 * zgratmp + ( zgradont + no3rat3 * ztmp1 ) * (1.0 - ssigma2)
-               zgrarep = po4rat3 * zgratmp + ( zgradopt + po4rat3 * ztmp1 ) * (1.0 - ssigma2)
-               zgraref = zgrareft + ferat3 * ( ztmp1 + zgratmp )
-
-               !   Defecation as a result of non assimilated products
-               !   --------------------------------------------------
-               zgrapoc  = zgraztotc * unass2c + unass2c / ( 1. - 0.8 * epsher2 ) * ztortz
-               zgrapon  = zgraztotn * unass2n + no3rat3 * unass2n / ( 1. - 0.8 * epsher2 ) * ztortz
-               zgrapop  = zgraztotp * unass2p + po4rat3 * unass2p / ( 1. - 0.8 * epsher2 ) * ztortz
-               zgrapof  = zgraztotf * unass2c + ferat3  * unass2c / ( 1. - 0.8 * epsher2 ) * ztortz
-
-               !  Addition of respiration to the release of inorganic nutrients
-               !  -------------------------------------------------------------
-               zgrarem = zgrarem + zbasresi + zrespirc
-               zgraren = zgraren + zbasresn + zrespirc * no3rat3
-               zgrarep = zgrarep + zbasresp + zrespirc * po4rat3
-               zgraref = zgraref + zbasresf + zrespirc * ferat3
-
-               !   Update the arrays TRA which contain the biological sources and
-               !   sinks
-               !   --------------------------------------------------------------
-               tr(ji,jj,jk,jppo4,Krhs) = tr(ji,jj,jk,jppo4,Krhs) + zgrarep 
-               tr(ji,jj,jk,jpnh4,Krhs) = tr(ji,jj,jk,jpnh4,Krhs) + zgraren
-               tr(ji,jj,jk,jpdoc,Krhs) = tr(ji,jj,jk,jpdoc,Krhs) + zgradoc
-               !
-               IF( ln_ligand ) THEN
-                  tr(ji,jj,jk,jplgw,Krhs)  = tr(ji,jj,jk,jplgw,Krhs) + zgradoc * ldocz
-                  zz2ligprod(ji,jj,jk) = zgradoc * ldocz
-               ENDIF
-               !
-               tr(ji,jj,jk,jpdon,Krhs) = tr(ji,jj,jk,jpdon,Krhs) + zgradon
-               tr(ji,jj,jk,jpdop,Krhs) = tr(ji,jj,jk,jpdop,Krhs) + zgradop
-               tr(ji,jj,jk,jpoxy,Krhs) = tr(ji,jj,jk,jpoxy,Krhs) - o2ut * zgrarem
-               tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + zgraref
-               zfezoo2(ji,jj,jk)   = zgraref
-               tr(ji,jj,jk,jpdic,Krhs) = tr(ji,jj,jk,jpdic,Krhs) + zgrarem
-               tr(ji,jj,jk,jptal,Krhs) = tr(ji,jj,jk,jptal,Krhs) + rno3 * zgraren
-               tr(ji,jj,jk,jpmes,Krhs) = tr(ji,jj,jk,jpmes,Krhs) + zepsherv * zgraztotc - zrespirc   &
-               &                     - ztortz - zgrazm
-               tr(ji,jj,jk,jpdia,Krhs) = tr(ji,jj,jk,jpdia,Krhs) - zgrazdc
-               tr(ji,jj,jk,jpndi,Krhs) = tr(ji,jj,jk,jpndi,Krhs) - zgrazdn
-               tr(ji,jj,jk,jppdi,Krhs) = tr(ji,jj,jk,jppdi,Krhs) - zgrazdp
-               tr(ji,jj,jk,jpdfe,Krhs) = tr(ji,jj,jk,jpdfe,Krhs) - zgrazdf
-               tr(ji,jj,jk,jpzoo,Krhs) = tr(ji,jj,jk,jpzoo,Krhs) - zgrazz
-               tr(ji,jj,jk,jpphy,Krhs) = tr(ji,jj,jk,jpphy,Krhs) - zgraznc
-               tr(ji,jj,jk,jpnph,Krhs) = tr(ji,jj,jk,jpnph,Krhs) - zgraznn
-               tr(ji,jj,jk,jppph,Krhs) = tr(ji,jj,jk,jppph,Krhs) - zgraznp
-               tr(ji,jj,jk,jpnfe,Krhs) = tr(ji,jj,jk,jpnfe,Krhs) - zgraznf
-               tr(ji,jj,jk,jpnch,Krhs) = tr(ji,jj,jk,jpnch,Krhs) - zgraznc * 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) - zgrazdc * 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) - zgrazdc * 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) + zgrazdc * tr(ji,jj,jk,jpdsi,Kbb) / ( tr(ji,jj,jk,jpdia,Kbb) + rtrn )
-
-               tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) - zgrazpoc - zgrazffep + zfracc
-               prodpoc(ji,jj,jk) = prodpoc(ji,jj,jk) + zfracc
-               conspoc(ji,jj,jk) = conspoc(ji,jj,jk) - zgrazpoc - zgrazffep
-               tr(ji,jj,jk,jppon,Krhs) = tr(ji,jj,jk,jppon,Krhs) - zgrazpon - zgrazffnp + zfracn
-               tr(ji,jj,jk,jppop,Krhs) = tr(ji,jj,jk,jppop,Krhs) - zgrazpop - zgrazffpp + zfracp
-               tr(ji,jj,jk,jpgoc,Krhs) = tr(ji,jj,jk,jpgoc,Krhs) - zgrazffeg + zgrapoc - zfracc
-               prodgoc(ji,jj,jk) = prodgoc(ji,jj,jk) + zgrapoc
-               consgoc(ji,jj,jk) = consgoc(ji,jj,jk) - zgrazffeg - zfracc
-               tr(ji,jj,jk,jpgon,Krhs) = tr(ji,jj,jk,jpgon,Krhs) - zgrazffng + zgrapon - zfracn
-               tr(ji,jj,jk,jpgop,Krhs) = tr(ji,jj,jk,jpgop,Krhs) - zgrazffpg + zgrapop - zfracp
-               tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) - zgrazpof - zgrazfffp + zfracfe
-               tr(ji,jj,jk,jpbfe,Krhs) = tr(ji,jj,jk,jpbfe,Krhs) - zgrazfffg + zgrapof - zfracfe
-               zfracal = tr(ji,jj,jk,jpcal,Kbb) / ( tr(ji,jj,jk,jpgoc,Kbb) + rtrn )
-               zgrazcal = zgrazffeg * (1. - part2) * zfracal
-
-               !  calcite production
-               !  ------------------
-               zprcaca = xfracal(ji,jj,jk) * zgraznc
-               prodcal(ji,jj,jk) = prodcal(ji,jj,jk) + zprcaca  ! prodcal=prodcal(nanophy)+prodcal(microzoo)+prodcal(mesozoo)
-               zprcaca = part2 * zprcaca
-               tr(ji,jj,jk,jpdic,Krhs) = tr(ji,jj,jk,jpdic,Krhs) + zgrazcal - zprcaca
-               tr(ji,jj,jk,jptal,Krhs) = tr(ji,jj,jk,jptal,Krhs) + 2. * ( zgrazcal - zprcaca )
-               tr(ji,jj,jk,jpcal,Krhs) = tr(ji,jj,jk,jpcal,Krhs) - zgrazcal + zprcaca
-            END DO
-         END DO
-      END DO
+      DO_3D_11_11( 1, jpkm1 )
+         zcompam   = MAX( ( tr(ji,jj,jk,jpmes,Kbb) - 1.e-9 ), 0.e0 )
+         zfact     = xstep * tgfunc2(ji,jj,jk) * zcompam
+
+         !   Michaelis-Menten mortality rates of mesozooplankton
+         !   ---------------------------------------------------
+         zrespz   = resrat2 * zfact * ( tr(ji,jj,jk,jpmes,Kbb) / ( xkmort + tr(ji,jj,jk,jpmes,Kbb) )  &
+         &          + 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   = mzrat2 * 1.e6 * zfact * tr(ji,jj,jk,jpmes,Kbb) * (1. - nitrfac(ji,jj,jk))
+
+         !   Computation of the abundance of the preys
+         !   A threshold can be specified in the namelist
+         !   --------------------------------------------
+         zcompadi  = MAX( ( tr(ji,jj,jk,jpdia,Kbb) - xthresh2dia ), 0.e0 )
+         zcompaz   = MAX( ( tr(ji,jj,jk,jpzoo,Kbb) - xthresh2zoo ), 0.e0 )
+         zcompaph  = MAX( ( tr(ji,jj,jk,jpphy,Kbb) - xthresh2phy ), 0.e0 )
+         zcompapoc = MAX( ( tr(ji,jj,jk,jppoc,Kbb) - xthresh2poc ), 0.e0 )
+         zcompames = MAX( ( tr(ji,jj,jk,jpmes,Kbb) - xthresh2mes ), 0.e0 )
+
+         !   Mesozooplankton grazing
+         !   ------------------------
+         zfood     = xpref2d * zcompadi + xpref2z * zcompaz + xpref2n * zcompaph + xpref2c * zcompapoc   &
+         &           + xpref2m * zcompames 
+         zfoodlim  = MAX( 0., zfood - MIN( 0.5 * zfood, xthresh2 ) )
+         zdenom    = zfoodlim / ( xkgraz2 + zfoodlim )
+         zgraze2   = grazrat2 * xstep * tgfunc2(ji,jj,jk) * tr(ji,jj,jk,jpmes,Kbb) * (1. - nitrfac(ji,jj,jk)) 
+
+         !   An active switching parameterization is used here.
+         !   We don't use the KTW parameterization proposed by 
+         !   Vallina et al. because it tends to produce to steady biomass
+         !   composition and the variance of Chl is too low as it grazes
+         !   too strongly on winning organisms. Thus, instead of a square
+         !   a 1.5 power value is used which decreases the pressure on the
+         !   most abundant species
+         !   ------------------------------------------------------------  
+         ztmp1 = xpref2n * zcompaph**1.5
+         ztmp2 = xpref2m * zcompames**1.5
+         ztmp3 = xpref2c * zcompapoc**1.5
+         ztmp4 = xpref2d * zcompadi**1.5
+         ztmp5 = xpref2z * zcompaz**1.5
+         ztmptot = ztmp1 + ztmp2 + ztmp3 + ztmp4 + ztmp5 + rtrn
+         ztmp1 = ztmp1 / ztmptot
+         ztmp2 = ztmp2 / ztmptot
+         ztmp3 = ztmp3 / ztmptot
+         ztmp4 = ztmp4 / ztmptot
+         ztmp5 = ztmp5 / ztmptot
+
+         !   Mesozooplankton regular grazing on the different preys
+         !   ------------------------------------------------------
+         zgrazdc   = zgraze2 * ztmp4 * zdenom
+         zgrazdn   = zgrazdc * tr(ji,jj,jk,jpndi,Kbb) / ( tr(ji,jj,jk,jpdia,Kbb) + rtrn)
+         zgrazdp   = zgrazdc * tr(ji,jj,jk,jppdi,Kbb) / ( tr(ji,jj,jk,jpdia,Kbb) + rtrn)
+         zgrazdf   = zgrazdc * tr(ji,jj,jk,jpdfe,Kbb) / ( tr(ji,jj,jk,jpdia,Kbb) + rtrn)
+         zgrazz    = zgraze2 * ztmp5 * zdenom
+         zgrazm    = zgraze2 * ztmp2 * zdenom
+         zgraznc   = zgraze2 * ztmp1 * zdenom
+         zgraznn   = zgraznc * tr(ji,jj,jk,jpnph,Kbb) / ( tr(ji,jj,jk,jpphy,Kbb) + rtrn)
+         zgraznp   = zgraznc * tr(ji,jj,jk,jppph,Kbb) / ( tr(ji,jj,jk,jpphy,Kbb) + rtrn)
+         zgraznf   = zgraznc * tr(ji,jj,jk,jpnfe,Kbb) / ( tr(ji,jj,jk,jpphy,Kbb) + rtrn)
+         zgrazpoc  = zgraze2 * ztmp3 * zdenom
+         zgrazpon  = zgrazpoc * tr(ji,jj,jk,jppon,Kbb) / ( tr(ji,jj,jk,jppoc,Kbb) + rtrn)
+         zgrazpop  = zgrazpoc * tr(ji,jj,jk,jppop,Kbb) / ( tr(ji,jj,jk,jppoc,Kbb) + rtrn)
+         zgrazpof  = zgrazpoc * tr(ji,jj,jk,jpsfe,Kbb) / ( tr(ji,jj,jk,jppoc,Kbb) + rtrn)
+
+         !   Mesozooplankton flux feeding on GOC
+         !   ----------------------------------
+         zgrazffeg = grazflux  * xstep * wsbio4(ji,jj,jk)      &
+         &           * tgfunc2(ji,jj,jk) * tr(ji,jj,jk,jpgoc,Kbb) * tr(ji,jj,jk,jpmes,Kbb)  &
+         &           * (1. - nitrfac(ji,jj,jk))
+         zgrazfffg = zgrazffeg * tr(ji,jj,jk,jpbfe,Kbb) / (tr(ji,jj,jk,jpgoc,Kbb) + rtrn)
+         zgrazffng = zgrazffeg * tr(ji,jj,jk,jpgon,Kbb) / (tr(ji,jj,jk,jpgoc,Kbb) + rtrn)
+         zgrazffpg = zgrazffeg * tr(ji,jj,jk,jpgop,Kbb) / (tr(ji,jj,jk,jpgoc,Kbb) + rtrn)
+         zgrazffep = grazflux  * xstep *  wsbio3(ji,jj,jk)     &
+         &           * tgfunc2(ji,jj,jk) * tr(ji,jj,jk,jppoc,Kbb) * tr(ji,jj,jk,jpmes,Kbb)   &
+         &           * (1. - nitrfac(ji,jj,jk))
+         zgrazfffp = zgrazffep * tr(ji,jj,jk,jpsfe,Kbb) / (tr(ji,jj,jk,jppoc,Kbb) + rtrn)
+         zgrazffnp = zgrazffep * tr(ji,jj,jk,jppon,Kbb) / (tr(ji,jj,jk,jppoc,Kbb) + rtrn)
+         zgrazffpp = zgrazffep * tr(ji,jj,jk,jppop,Kbb) / (tr(ji,jj,jk,jppoc,Kbb) + rtrn)
+         !
+         zgraztotc  = zgrazdc + zgrazz + zgraznc + zgrazm + zgrazpoc + zgrazffep + zgrazffeg
+
+         !   Compute the proportion of filter feeders
+         !   ----------------------------------------  
+         zproport  = (zgrazffep + zgrazffeg)/(rtrn + zgraztotc)
+
+         !   Compute fractionation of aggregates. It is assumed that 
+         !   diatoms based aggregates are more prone to fractionation
+         !   since they are more porous (marine snow instead of fecal pellets)
+         !   ----------------------------------------------------------------
+         zratio    = tr(ji,jj,jk,jpgsi,Kbb) / ( tr(ji,jj,jk,jpgoc,Kbb) + rtrn )
+         zratio2   = zratio * zratio
+         zfracc    = zproport * grazflux  * xstep * wsbio4(ji,jj,jk)      &
+         &          * tr(ji,jj,jk,jpgoc,Kbb) * tr(ji,jj,jk,jpmes,Kbb)          &
+         &          * ( 0.2 + 3.8 * zratio2 / ( 1.**2 + zratio2 ) )
+         zfracfe   = zfracc * tr(ji,jj,jk,jpbfe,Kbb) / (tr(ji,jj,jk,jpgoc,Kbb) + rtrn)
+         zfracn    = zfracc * tr(ji,jj,jk,jpgon,Kbb) / (tr(ji,jj,jk,jpgoc,Kbb) + rtrn)
+         zfracp    = zfracc * tr(ji,jj,jk,jpgop,Kbb) / (tr(ji,jj,jk,jpgoc,Kbb) + rtrn)
+
+         zgrazffep = zproport * zgrazffep   ;   zgrazffeg = zproport * zgrazffeg
+         zgrazfffp = zproport * zgrazfffp   ;   zgrazfffg = zproport * zgrazfffg
+         zgrazffnp = zproport * zgrazffnp   ;   zgrazffng = zproport * zgrazffng
+         zgrazffpp = zproport * zgrazffpp   ;   zgrazffpg = zproport * zgrazffpg
+
+         zgraztotc  = zgrazdc + zgrazz + zgraznc + zgrazm + zgrazpoc + zgrazffep + zgrazffeg
+         zgraztotf  = zgrazdf + zgraznf + ( zgrazz + zgrazm ) * ferat3 + zgrazpof &
+         &            + zgrazfffp + zgrazfffg
+         zgraztotn  = zgrazdn + (zgrazm + zgrazz) * no3rat3 + zgraznn + zgrazpon  &
+         &            + zgrazffnp + zgrazffng
+         zgraztotp  = zgrazdp + (zgrazz + zgrazm) * po4rat3 + zgraznp + zgrazpop  &
+         &            + zgrazffpp + zgrazffpg
+
+
+         ! Total grazing ( grazing by microzoo is already computed in p5zmicro )
+         zgrazing(ji,jj,jk) = zgraztotc
+
+         !   Stoichiometruc ratios of the food ingested by zooplanton 
+         !   --------------------------------------------------------
+         zgrasratf  =  (zgraztotf + rtrn) / ( zgraztotc + rtrn )
+         zgrasratn  =  (zgraztotn + rtrn) / ( zgraztotc + rtrn )
+         zgrasratp  =  (zgraztotp + rtrn) / ( zgraztotc + rtrn )
+
+         !   Growth efficiency is made a function of the quality 
+         !   and the quantity of the preys
+         !   ---------------------------------------------------
+         zepshert  = MIN( 1., zgrasratn/ no3rat3, zgrasratp/ po4rat3, zgrasratf / ferat3)
+         zbeta     = MAX(0., (epsher2 - epsher2min) )
+         zepsherf  = epsher2min + zbeta / ( 1.0 + 0.04E6 * 12. * zfood * zbeta )
+         zepsherv  = zepsherf * zepshert
+
+         !   Respiration of mesozooplankton
+         !   Excess carbon in the food is used preferentially
+         !   ----------------  ------------------------------
+         zexcess  = zgraztotc * zepsherf * (1.0 - zepshert) * zmetexcess 
+         zbasresb = MAX(0., zrespz - zexcess)
+         zbasresi = zexcess + MIN(0., zrespz - zexcess)
+         zrespirc = srespir2 * zepsherv * zgraztotc + zbasresb
+
+         !   When excess carbon is used, the other elements in excess
+         !   are also used proportionally to their abundance
+         !   --------------------------------------------------------
+         zexcess  = ( zgrasratn/ no3rat3 - zepshert ) / ( 1.0 - zepshert + rtrn)
+         zbasresn = zbasresi * zexcess * zgrasratn
+         zexcess  = ( zgrasratp/ po4rat3 - zepshert ) / ( 1.0 - zepshert + rtrn)
+         zbasresp = zbasresi * zexcess * zgrasratp
+         zexcess  = ( zgrasratf/ ferat3 - zepshert ) / ( 1.0 - zepshert + rtrn)
+         zbasresf = zbasresi * zexcess * zgrasratf
+
+         !   Voiding of the excessive elements as organic matter
+         !   --------------------------------------------------------
+         zgradoct = (1. - unass2c - zepsherv) * zgraztotc - zbasresi
+         zgradont = (1. - unass2n) * zgraztotn - zepsherv * no3rat3 * zgraztotc - zbasresn
+         zgradopt = (1. - unass2p) * zgraztotp - zepsherv * po4rat3 * zgraztotc - zbasresp
+         zgrareft = (1. - unass2c) * zgraztotf - zepsherv * ferat3 * zgraztotc - zbasresf
+         ztmp1   = ( 1. - epsher2 - unass2c ) /( 1. - 0.8 * epsher2 ) * ztortz
+         zgradoc = (zgradoct + ztmp1) * ssigma2
+         zgradon = (zgradont + no3rat3 * ztmp1) * ssigma2
+         zgradop = (zgradopt + po4rat3 * ztmp1) * ssigma2
+         zgratmp = 0.2 * epsher2 /( 1. - 0.8 * epsher2 ) * ztortz
+
+         !  Since only semilabile DOM is represented in PISCES
+         !  part of DOM is in fact labile and is then released
+         !  as dissolved inorganic compounds (ssigma2)
+         !  --------------------------------------------------
+         zgrarem = zgratmp + ( zgradoct + ztmp1 ) * (1.0 - ssigma2)
+         zgraren = no3rat3 * zgratmp + ( zgradont + no3rat3 * ztmp1 ) * (1.0 - ssigma2)
+         zgrarep = po4rat3 * zgratmp + ( zgradopt + po4rat3 * ztmp1 ) * (1.0 - ssigma2)
+         zgraref = zgrareft + ferat3 * ( ztmp1 + zgratmp )
+
+         !   Defecation as a result of non assimilated products
+         !   --------------------------------------------------
+         zgrapoc  = zgraztotc * unass2c + unass2c / ( 1. - 0.8 * epsher2 ) * ztortz
+         zgrapon  = zgraztotn * unass2n + no3rat3 * unass2n / ( 1. - 0.8 * epsher2 ) * ztortz
+         zgrapop  = zgraztotp * unass2p + po4rat3 * unass2p / ( 1. - 0.8 * epsher2 ) * ztortz
+         zgrapof  = zgraztotf * unass2c + ferat3  * unass2c / ( 1. - 0.8 * epsher2 ) * ztortz
+
+         !  Addition of respiration to the release of inorganic nutrients
+         !  -------------------------------------------------------------
+         zgrarem = zgrarem + zbasresi + zrespirc
+         zgraren = zgraren + zbasresn + zrespirc * no3rat3
+         zgrarep = zgrarep + zbasresp + zrespirc * po4rat3
+         zgraref = zgraref + zbasresf + zrespirc * ferat3
+
+         !   Update the arrays TRA which contain the biological sources and
+         !   sinks
+         !   --------------------------------------------------------------
+         tr(ji,jj,jk,jppo4,Krhs) = tr(ji,jj,jk,jppo4,Krhs) + zgrarep 
+         tr(ji,jj,jk,jpnh4,Krhs) = tr(ji,jj,jk,jpnh4,Krhs) + zgraren
+         tr(ji,jj,jk,jpdoc,Krhs) = tr(ji,jj,jk,jpdoc,Krhs) + zgradoc
+         !
+         IF( ln_ligand ) THEN
+            tr(ji,jj,jk,jplgw,Krhs)  = tr(ji,jj,jk,jplgw,Krhs) + zgradoc * ldocz
+            zz2ligprod(ji,jj,jk) = zgradoc * ldocz
+         ENDIF
+         !
+         tr(ji,jj,jk,jpdon,Krhs) = tr(ji,jj,jk,jpdon,Krhs) + zgradon
+         tr(ji,jj,jk,jpdop,Krhs) = tr(ji,jj,jk,jpdop,Krhs) + zgradop
+         tr(ji,jj,jk,jpoxy,Krhs) = tr(ji,jj,jk,jpoxy,Krhs) - o2ut * zgrarem
+         tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + zgraref
+         zfezoo2(ji,jj,jk)   = zgraref
+         tr(ji,jj,jk,jpdic,Krhs) = tr(ji,jj,jk,jpdic,Krhs) + zgrarem
+         tr(ji,jj,jk,jptal,Krhs) = tr(ji,jj,jk,jptal,Krhs) + rno3 * zgraren
+         tr(ji,jj,jk,jpmes,Krhs) = tr(ji,jj,jk,jpmes,Krhs) + zepsherv * zgraztotc - zrespirc   &
+         &                     - ztortz - zgrazm
+         tr(ji,jj,jk,jpdia,Krhs) = tr(ji,jj,jk,jpdia,Krhs) - zgrazdc
+         tr(ji,jj,jk,jpndi,Krhs) = tr(ji,jj,jk,jpndi,Krhs) - zgrazdn
+         tr(ji,jj,jk,jppdi,Krhs) = tr(ji,jj,jk,jppdi,Krhs) - zgrazdp
+         tr(ji,jj,jk,jpdfe,Krhs) = tr(ji,jj,jk,jpdfe,Krhs) - zgrazdf
+         tr(ji,jj,jk,jpzoo,Krhs) = tr(ji,jj,jk,jpzoo,Krhs) - zgrazz
+         tr(ji,jj,jk,jpphy,Krhs) = tr(ji,jj,jk,jpphy,Krhs) - zgraznc
+         tr(ji,jj,jk,jpnph,Krhs) = tr(ji,jj,jk,jpnph,Krhs) - zgraznn
+         tr(ji,jj,jk,jppph,Krhs) = tr(ji,jj,jk,jppph,Krhs) - zgraznp
+         tr(ji,jj,jk,jpnfe,Krhs) = tr(ji,jj,jk,jpnfe,Krhs) - zgraznf
+         tr(ji,jj,jk,jpnch,Krhs) = tr(ji,jj,jk,jpnch,Krhs) - zgraznc * 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) - zgrazdc * 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) - zgrazdc * 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) + zgrazdc * tr(ji,jj,jk,jpdsi,Kbb) / ( tr(ji,jj,jk,jpdia,Kbb) + rtrn )
+
+         tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) - zgrazpoc - zgrazffep + zfracc
+         prodpoc(ji,jj,jk) = prodpoc(ji,jj,jk) + zfracc
+         conspoc(ji,jj,jk) = conspoc(ji,jj,jk) - zgrazpoc - zgrazffep
+         tr(ji,jj,jk,jppon,Krhs) = tr(ji,jj,jk,jppon,Krhs) - zgrazpon - zgrazffnp + zfracn
+         tr(ji,jj,jk,jppop,Krhs) = tr(ji,jj,jk,jppop,Krhs) - zgrazpop - zgrazffpp + zfracp
+         tr(ji,jj,jk,jpgoc,Krhs) = tr(ji,jj,jk,jpgoc,Krhs) - zgrazffeg + zgrapoc - zfracc
+         prodgoc(ji,jj,jk) = prodgoc(ji,jj,jk) + zgrapoc
+         consgoc(ji,jj,jk) = consgoc(ji,jj,jk) - zgrazffeg - zfracc
+         tr(ji,jj,jk,jpgon,Krhs) = tr(ji,jj,jk,jpgon,Krhs) - zgrazffng + zgrapon - zfracn
+         tr(ji,jj,jk,jpgop,Krhs) = tr(ji,jj,jk,jpgop,Krhs) - zgrazffpg + zgrapop - zfracp
+         tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) - zgrazpof - zgrazfffp + zfracfe
+         tr(ji,jj,jk,jpbfe,Krhs) = tr(ji,jj,jk,jpbfe,Krhs) - zgrazfffg + zgrapof - zfracfe
+         zfracal = tr(ji,jj,jk,jpcal,Kbb) / ( tr(ji,jj,jk,jpgoc,Kbb) + rtrn )
+         zgrazcal = zgrazffeg * (1. - part2) * zfracal
+
+         !  calcite production
+         !  ------------------
+         zprcaca = xfracal(ji,jj,jk) * zgraznc
+         prodcal(ji,jj,jk) = prodcal(ji,jj,jk) + zprcaca  ! prodcal=prodcal(nanophy)+prodcal(microzoo)+prodcal(mesozoo)
+         zprcaca = part2 * zprcaca
+         tr(ji,jj,jk,jpdic,Krhs) = tr(ji,jj,jk,jpdic,Krhs) + zgrazcal - zprcaca
+         tr(ji,jj,jk,jptal,Krhs) = tr(ji,jj,jk,jptal,Krhs) + 2. * ( zgrazcal - zprcaca )
+         tr(ji,jj,jk,jpcal,Krhs) = tr(ji,jj,jk,jpcal,Krhs) - zgrazcal + zprcaca
+      END_3D
       !
       IF( lk_iomput .AND. knt == nrdttrc ) THEN
