Index: NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zmeso.F90
===================================================================
--- NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zmeso.F90	(revision 12524)
+++ NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zmeso.F90	(revision 12537)
@@ -2,5 +2,5 @@
    !!======================================================================
    !!                         ***  MODULE p5zmeso  ***
-   !! TOP :   PISCES Compute the sources/sinks for mesozooplankton
+   !! TOP :   PISCES-QUOTA Compute the sources/sinks for mesozooplankton
    !!======================================================================
    !! History :   1.0  !  2002     (O. Aumont) Original code
@@ -9,6 +9,7 @@
    !!             3.6  !  2015-05  (O. Aumont) PISCES quota
    !!----------------------------------------------------------------------
-   !!   p5z_meso       :   Compute the sources/sinks for mesozooplankton
-   !!   p5z_meso_init  :   Initialization of the parameters for mesozooplankton
+   !!   p5z_meso       : Compute the sources/sinks for mesozooplankton
+   !!   p5z_meso_init  : Initialization of the parameters for mesozooplankton
+   !!   p5z_meso_alloc : Allocate variables for mesozooplankton 
    !!----------------------------------------------------------------------
    USE oce_trc         !  shared variables between ocean and passive tracers
@@ -23,5 +24,5 @@
    PUBLIC   p5z_meso              ! called in p5zbio.F90
    PUBLIC   p5z_meso_init         ! called in trcsms_pisces.F90
-   PUBLIC   p5z_meso_alloc
+   PUBLIC   p5z_meso_alloc        ! called in trcini_pisces.F90
 
    !! * Shared module variables
@@ -53,7 +54,6 @@
    LOGICAL,  PUBLIC ::  bmetexc2     !: Use of excess carbon for respiration
    LOGICAL , PUBLIC ::  ln_dvm_meso  !: Boolean to activate DVM of mesozooplankton
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: depmig
-   INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:) :: kmig
-
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: depmig  !: DVM of mesozooplankton : migration depth
+   INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:) :: kmig    !: Vertical indice of the the migration depth
 
    !!----------------------------------------------------------------------
@@ -70,8 +70,14 @@
       !!
       !! ** Purpose :   Compute the sources/sinks for mesozooplankton
-      !!
-      !! ** Method  : - ???
+      !!                This includes ingestion and assimilation, flux feeding
+      !!                and mortality. We use an active prey switching  
+      !!                parameterization Morozov and Petrovskii (2013). 
+      !!                All living compartments and mesozooplankton
+      !!                are potential preys of mesozooplankton as well as small
+      !!                sinking particles 
+      !!
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt, knt ! ocean time step
+      !
       INTEGER  :: ji, jj, jk, jkt
       REAL(wp) :: zcompadi, zcompaph, zcompapoc, zcompaz, zcompam, zcompames
@@ -115,7 +121,9 @@
       !
       ! Diurnal vertical migration of mesozooplankton
+      ! Computation of the migration depth
       ! ---------------------------------------------
       IF (ln_dvm_meso) CALL p5z_meso_depmig
 
+      ! Use of excess carbon for metabolism
       zmetexcess = 0.0
       IF ( bmetexc2 ) zmetexcess = 1.0
@@ -127,12 +135,15 @@
                zfact     = xstep * tgfunc2(ji,jj,jk) * zcompam
 
-               !   Michaelis-Menten mortality rates of mesozooplankton
-               !   ---------------------------------------------------
+               !  linear mortality of mesozooplankton
+               !  A michaelis menten modulation term is used to avoid extinction of 
+               !  mesozooplankton at very low food concentrations
+               !  -----------------------------------------------------------------
                zrespz   = resrat2 * zfact * ( trb(ji,jj,jk,jpmes) / ( xkmort + trb(ji,jj,jk,jpmes) )  &
                &          + 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
-               !   ---------------------------------------------------------------
+               !  Zooplankton quadratic mortality. A square function has been selected with
+               !  to mimic predation and disease (density dependent mortality). It also tends
+               !  to stabilise the model
+               !  -------------------------------------------------------------------------
                ztortz   = mzrat2 * 1.e6 * zfact * trb(ji,jj,jk,jpmes) * (1. - nitrfac(ji,jj,jk))
 
@@ -147,5 +158,10 @@
 
                !   Mesozooplankton grazing
-               !   ------------------------
+               ! The total amount of food is the sum of all preys accessible to mesozooplankton 
+               ! multiplied by their food preference
+               ! A threshold can be specified in the namelist (xthresh2). However, when food 
+               ! concentration is close to this threshold, it is decreased to avoid the 
+               ! accumulation of food in the mesozoopelagic domain
+               ! -------------------------------------------------------------------------------
                zfood     = xpref2d * zcompadi + xpref2z * zcompaz + xpref2n * zcompaph + xpref2c * zcompapoc   &
                &           + xpref2m * zcompames 
@@ -154,14 +170,21 @@
                zgraze2   = grazrat2 * xstep * tgfunc2(ji,jj,jk) * trb(ji,jj,jk,jpmes) * (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 too steady biomass
-               !   composition and the variance of Chl is too low as it grazes
-               !   too strongly on winning organisms. We use a generalized
-               !   switching parameterization proposed by Morozov and 
-               !   Petrovskii (2013)
-               !   ------------------------------------------------------------  
+               ! An active switching parameterization is used here.
+               ! We don't use the KTW parameterization proposed by 
+               ! Vallina et al. because it tends to produce too steady biomass
+               ! composition and the variance of Chl is too low as it grazes
+               ! too strongly on winning organisms. We use a generalized
+               ! switching parameterization proposed by Morozov and 
+               ! Petrovskii (2013)
+               ! ------------------------------------------------------------  
+               ! The width of the selection window is increased when preys
+               ! have low abundance, .i.e. zooplankton become less specific 
+               ! to avoid starvation.
+               ! ----------------------------------------------------------
                zsigma = 1.0 - zdenom**2/(0.05**2+zdenom**2)
                zsigma = 0.5 + 1.0 * zsigma
+               ! Nanophytoplankton and diatoms are the only preys considered
+               ! to be close enough to have potential interference
+               ! -----------------------------------------------------------
                zdiffdn = exp( -ABS(log(3.0 * sizen(ji,jj,jk) / (5.0 * sized(ji,jj,jk) + rtrn )) )**2 / zsigma**2 )
                ztmp1 = xpref2n * zcompaph * ( zcompaph + zdiffdn * zcompadi )
@@ -194,6 +217,7 @@
                zgrazpof  = zgrazpoc * trb(ji,jj,jk,jpsfe) / ( trb(ji,jj,jk,jppoc) + rtrn)
 
-               !   Mesozooplankton flux feeding on GOC
-               !   ----------------------------------
+               !  Mesozooplankton flux feeding on GOC and POC. The feeding pressure
+               ! is proportional to the flux
+               !  ------------------------------------------------------------------
                zgrazffeg = grazflux  * xstep * wsbio4(ji,jj,jk)      &
                &           * tgfunc2(ji,jj,jk) * trb(ji,jj,jk,jpgoc) * trb(ji,jj,jk,jpmes)  &
@@ -211,6 +235,6 @@
                zgraztotc  = zgrazdc + zgrazz + zgraznc + zgrazm + zgrazpoc + zgrazffep + zgrazffeg
 
-               !   Compute the proportion of filter feeders
-               !   ----------------------------------------  
+               ! Compute the proportion of filter feeders. It is assumed steady state.
+               ! ---------------------------------------------------------------------  
                zproport  = (zgrazffep + zgrazffeg)/(rtrn + zgraztotc)
 
@@ -228,4 +252,5 @@
                zfracp    = zfracc * trb(ji,jj,jk,jpgop) / (trb(ji,jj,jk,jpgoc) + rtrn)
 
+               ! Flux feeding is multiplied by the fractional biomass of flux feeders
                zgrazffep = zproport * zgrazffep   ;   zgrazffeg = zproport * zgrazffeg
                zgrazfffp = zproport * zgrazfffp   ;   zgrazfffg = zproport * zgrazfffg
@@ -241,5 +266,4 @@
                &            + zgrazffpp + zgrazffpg
 
-
                ! Total grazing ( grazing by microzoo is already computed in p5zmicro )
                zgrazing(ji,jj,jk) = zgraztotc
@@ -251,7 +275,12 @@
                zgrasratp  =  (zgraztotp + rtrn) / ( zgraztotc + rtrn )
 
-               !   Growth efficiency is made a function of the quality 
-               !   and the quantity of the preys
-               !   ---------------------------------------------------
+               ! Mesozooplankton efficiency. 
+               ! We adopt a formulation proposed by Mitra et al. (2007)
+               ! The gross growth efficiency is controled by the most limiting nutrient.
+               ! Growth is also further decreased when the food quality is poor. This is currently
+               ! hard coded : it can be decreased by up to 50% (zepsherq)
+               ! GGE can also be decreased when food quantity is high, zepsherf (Montagnes and 
+               ! Fulton, 2012)
+               ! -----------------------------------------------------------------------------------
                zepshert  = MIN( 1., zgrasratn/ no3rat3, zgrasratp/ po4rat3, zgrasratf / ferat3)
                zbeta     = MAX(0., (epsher2 - epsher2min) )
@@ -261,5 +290,6 @@
                !   Respiration of mesozooplankton
                !   Excess carbon in the food is used preferentially
-               !   ----------------  ------------------------------
+               !   when bmetexc2 is set to .true.
+               !   -----------------------------------------------
                zexcess  = zgraztotc * zepsherf * (1.0 - zepshert) * zmetexcess 
                zbasresb = MAX(0., zrespz - zexcess)
@@ -344,6 +374,8 @@
                zfracal = trb(ji,jj,jk,jpcal) / ( trb(ji,jj,jk,jpgoc) + rtrn )
                zgrazcal = zgrazffeg * (1. - part2) * zfracal
-               !  calcite production
-               !  ------------------
+               ! Calcite production
+               ! Calcite remineralization due to zooplankton activity
+               ! part2 of the ingested calcite is dissolving in the acidic gut
+               ! -------------------------------------------------------------
                zprcaca = xfracal(ji,jj,jk) * zgraznc
                prodcal(ji,jj,jk) = prodcal(ji,jj,jk) + zprcaca  ! prodcal=prodcal(nanophy)+prodcal(microzoo)+prodcal(mesozoo)
@@ -356,9 +388,12 @@
       END DO
 
+      ! Computation of the effect of DVM by mesozooplankton
+      ! This part is only activated if ln_dvm_meso is set to true
+      ! The parameterization has been published in Gorgues et al. (2019).
+      ! -----------------------------------------------------------------
       IF (ln_dvm_meso) THEN
          ALLOCATE( zgramigrem(jpi,jpj), zgramigref(jpi,jpj), zgramigpoc(jpi,jpj), zgramigpof(jpi,jpj) )
          ALLOCATE( zgramigrep(jpi,jpj), zgramigren(jpi,jpj), zgramigpop(jpi,jpj), zgramigpon(jpi,jpj) )
          ALLOCATE( zgramigdoc(jpi,jpj), zgramigdon(jpi,jpj), zgramigdop(jpi,jpj) )
-
          ALLOCATE( zstrn(jpi,jpj) )
          zgramigrem(:,:)  = 0.0   ;   zgramigref(:,:) = 0.0
@@ -385,9 +420,11 @@
 
 
+        ! Compute the amount of materials that will go into vertical migration
+        ! This fraction is sumed over the euphotic zone and is removed from 
+        ! the fluxes driven by mesozooplankton in the euphotic zone.
+        ! --------------------------------------------------------------------
          DO jk = 1, jpk
             DO jj = 1, jpj
                DO ji = 1, jpi
-
-                  !   Compute the amount of materials that will go into vertical migration
                   zmigreltime = (1. - zstrn(ji,jj))
                   IF ( gdept_n(ji,jj,jk) <= heup(ji,jj) ) THEN
@@ -431,4 +468,7 @@
          END DO
 
+         ! The inorganic and organic fluxes induced by migrating organisms are added at the 
+         ! the migration depth (corresponding indice is set by kmig)
+         ! --------------------------------------------------------------------------------
          DO jj = 1, jpj
             DO ji = 1, jpi
@@ -456,11 +496,14 @@
          DEALLOCATE( zgramigdoc, zgramigdon, zgramigdop )
          DEALLOCATE( zstrn )
-
+      ! End of the ln_dvm_meso part
       ENDIF
 
+      !   Update the arrays TRA which contain the biological sources and sinks
+      !   This only concerns the variables which are affected by DVM (inorganic 
+      !   nutrients, DOC agands, and particulate organic carbon). 
+      !   ---------------------------------------------------------------------
       DO jk = 1, jpk
          DO jj = 1, jpj
             DO ji = 1, jpi
-               !   Update the arrays TRA which contain the biological sources and sinks
                tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) + zgrarep(ji,jj,jk) 
                tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zgraren(ji,jj,jk)
@@ -526,11 +569,11 @@
       !! ** Purpose :   Initialization of mesozooplankton parameters
       !!
-      !! ** Method  :   Read the nampismes namelist and check the parameters
+      !! ** Method  :   Read the namp5zmes namelist and check the parameters
       !!      called at the first timestep (nittrc000)
       !!
-      !! ** input   :   Namelist nampismes
+      !! ** input   :   Namelist namp5zmes
       !!
       !!----------------------------------------------------------------------
-      INTEGER :: ios                 ! Local integer output status for namelist read
+      INTEGER :: ios    ! Local integer output status for namelist read
       !!
       NAMELIST/namp5zmes/part2, bmetexc2, grazrat2, resrat2, mzrat2, xpref2c, xpref2n, xpref2z, &
@@ -540,11 +583,11 @@
       !!----------------------------------------------------------------------
       !
-      REWIND( numnatp_ref )              ! Namelist nampismes in reference namelist : Pisces mesozooplankton
+      REWIND( numnatp_ref )              ! Namelist namp5zmes in reference namelist : Pisces mesozooplankton
       READ  ( numnatp_ref, namp5zmes, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampismes in reference namelist' )
-      !
-      REWIND( numnatp_cfg )              ! Namelist nampismes in configuration namelist : Pisces mesozooplankton
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namp5zmes in reference namelist' )
+      !
+      REWIND( numnatp_cfg )              ! Namelist namp5zmes in configuration namelist : Pisces mesozooplankton
       READ  ( numnatp_cfg, namp5zmes, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nampismes in configuration namelist' )
+902   IF( ios >  0 ) CALL ctl_nam ( ios , 'namp5zmes in configuration namelist' )
       IF(lwm) WRITE ( numonp, namp5zmes )
       !
@@ -628,4 +671,7 @@
       END DO
 
+      ! Compute the difference between surface values and the mean values in the mesopelagic
+      ! domain
+      ! ------------------------------------------------------------------------------------
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -634,7 +680,8 @@
          END DO
       END DO
-      !
+
       ! Computation of the migration depth based on the parameterization of 
-      !
+      ! Bianchi et al. (2013)
+      ! -------------------------------------------------------------------
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -648,5 +695,4 @@
       ! Computation of the corresponding jk indice 
       ! ------------------------------------------
-      ! 
       DO jk = 1, jpk-1
          DO jj = 1, jpj
@@ -665,5 +711,4 @@
       ! that it falls above the OMZ
       ! -----------------------------------------------------------------------
-      !
       DO ji =1, jpi
          DO jj = 1, jpj
@@ -685,5 +730,5 @@
    INTEGER FUNCTION p5z_meso_alloc()
       !!----------------------------------------------------------------------
-      !!                     ***  ROUTINE p4z_meso_alloc  ***
+      !!                     ***  ROUTINE p5z_meso_alloc  ***
       !!----------------------------------------------------------------------
       !
