Changeset 9450

Timestamp:

2018-03-30T18:08:29+02:00 (7 years ago)

Author:

aumont

Message:

debug PISCES code

Location:

branches/CNRS/dev_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES

Files:

• P4Z/p4zbio.F90 (modified) (1 diff)
• P4Z/p4zfechem.F90 (modified) (7 diffs)
• P4Z/p4zligand.F90 (modified) (4 diffs)
• P4Z/p4zlim.F90 (modified) (6 diffs)
• P4Z/p4zlys.F90 (modified) (1 diff)
• P4Z/p4zmeso.F90 (modified) (5 diffs)
• P4Z/p4zmicro.F90 (modified) (5 diffs)
• P4Z/p4zpoc.F90 (modified) (9 diffs)
• P4Z/p4zprod.F90 (modified) (18 diffs)
• P4Z/p4zrem.F90 (modified) (6 diffs)
• P4Z/p4zsbc.F90 (modified) (1 diff)
• P4Z/p5zmeso.F90 (modified) (7 diffs)
• P4Z/p5zmicro.F90 (modified) (8 diffs)
• P4Z/p5zmort.F90 (modified) (1 diff)
• P4Z/p5zprod.F90 (modified) (14 diffs)
• sms_pisces.F90 (modified) (2 diffs)
• trcini_pisces.F90 (modified) (4 diffs)

branches/CNRS/dev_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zbio.F90

@@ -20 +20 @@
    USE p4zmicro        !  Sources and sinks of microzooplankton
    USE p4zmeso         !  Sources and sinks of mesozooplankton
-   USE p5zlim          !  Co-limitations of differents nutrients
    USE p5zprod         !  Growth rate of the 2 phyto groups
    USE p5zmort         !  Mortality terms for phytoplankton

branches/CNRS/dev_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90

@@ -73 +73 @@
       REAL(wp) :: zaggliga, zaggligb
       REAL(wp) :: dissol, zligco
+      REAL(wp) :: zrfact2
       CHARACTER (len=25) :: charout
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zTL1, zFe3, ztotlig, precip
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zFeL1, zFeL2, zTL2, zFe2, zFeP
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zcoll3d, zscav3d, zlcoll3d
       REAL(wp), POINTER, DIMENSION(:,:  ) :: zstrn, zstrn2
       !!---------------------------------------------------------------------
@@ -83 +85 @@
       ! Allocate temporary workspace
       CALL wrk_alloc( jpi, jpj, jpk, zFe3, zFeL1, zTL1, ztotlig, precip )
+      CALL wrk_alloc( jpi, jpj, jpk, zcoll3d, zscav3d )
+      IF( ln_ligand ) THEN
+         CALL wrk_alloc( jpi, jpj, jpk, zlcoll3d )
+      ENDIF
       zFe3 (:,:,:) = 0.
       zFeL1(:,:,:) = 0.
@@ -301 +307 @@
                tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zscave * zdenom1 + zaggdfea
                tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zscave * zdenom2 + zaggdfeb
+               zscav3d(ji,jj,jk)   = zscave
+               zcoll3d(ji,jj,jk)   = zaggdfea + zaggdfeb
                !
             END DO
@@ -321 +329 @@
                   !
                   zlam1b   = 3.53E3 *   trb(ji,jj,jk,jpgoc) * xdiss(ji,jj,jk)
-                  zligco   = MAX( ( 0.1 * trb(ji,jj,jk,jplgw) ), ( trb(ji,jj,jk,jplgw) - fe3sol ) )
+                  IF ( ln_fecolloid) THEN
+                     zligco   = MAX( ( 0.1 * trb(ji,jj,jk,jplgw) ), ( trb(ji,jj,jk,jplgw) - fe3sol ) )
+                  ELSE
+                     zligco  = 0.5 * trn(ji,jj,jk,jplgw)
+                  ENDIF
                   zaggliga = zlam1a * xstep * zligco
                   zaggligb = zlam1b * xstep * zligco
                   tra(ji,jj,jk,jpfep) = tra(ji,jj,jk,jpfep) + precip(ji,jj,jk)
                   tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) - zaggliga - zaggligb
+                  zlcoll3d(ji,jj,jk)  = zaggliga + zaggligb
                END DO
             END DO
@@ -340 +353 @@
       IF( lk_iomput ) THEN
          IF( knt == nrdttrc ) THEN
+         zrfact2 = 1.e3 * rfact2r  ! conversion from mol/L/timestep into mol/m3/s
          IF( iom_use("Fe3")    )  CALL iom_put("Fe3"    , zFe3   (:,:,:)       * tmask(:,:,:) )   ! Fe3+
          IF( iom_use("FeL1")   )  CALL iom_put("FeL1"   , zFeL1  (:,:,:)       * tmask(:,:,:) )   ! FeL1
@@ -345 +359 @@
          IF( iom_use("Totlig") )  CALL iom_put("Totlig" , ztotlig(:,:,:)       * tmask(:,:,:) )   ! TL
          IF( iom_use("Biron")  )  CALL iom_put("Biron"  , biron  (:,:,:) * 1e9 * tmask(:,:,:) )   ! biron
+         IF( iom_use("FESCAV") )  CALL iom_put("FESCAV" , zscav3d(:,:,:)  * 1e9 * tmask(:,:,:) * zrfact2 )
+         IF( iom_use("FECOLL") )  CALL iom_put("FECOLL" , zcoll3d(:,:,:)  * 1e9 * tmask(:,:,:) * zrfact2 )
+         IF( iom_use("LGWCOLL"))  CALL iom_put("LGWCOLL", zlcoll3d(:,:,:) * 1e9 * tmask(:,:,:) * zrfact2 )
          IF( ln_fechem ) THEN
             IF( iom_use("Fe2")  ) CALL iom_put("Fe2"    , zFe2   (:,:,:)       * tmask(:,:,:) )   ! Fe2+
@@ -361 +378 @@
       !
       CALL wrk_dealloc( jpi, jpj, jpk, zFe3, zFeL1, zTL1, ztotlig, precip )
+      CALL wrk_dealloc( jpi, jpj, jpk, zscav3d, zcoll3d )
       IF( ln_fechem )  THEN
          CALL wrk_dealloc( jpi, jpj,      zstrn, zstrn2 )
          CALL wrk_dealloc( jpi, jpj, jpk, zFe2, zFeL2, zTL2, zFeP )
+      ENDIF
+      IF ( ln_ligand ) THEN
+         CALL wrk_dealloc( jpi, jpj, jpk, zlcoll3d )
       ENDIF
       !

branches/CNRS/dev_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zligand.F90

@@ -13 +13 @@
    USE sms_pisces      !  PISCES Source Minus Sink variables
    USE prtctl_trc      !  print control for debugging
+   USE iom             !  I/O manager
 
    IMPLICIT NONE
@@ -47 +48 @@
       INTEGER  ::   ji, jj, jk
       REAL(wp) ::   zlgwp, zlgwpr, zlgwr, zlablgw, zrfepa, zfepr
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zligrem, zligpr, zrligprod
       CHARACTER (len=25) :: charout
       !!---------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('p4z_ligand')
+      !
+      CALL wrk_alloc( jpi, jpj, jpk, zligrem, zligpr, zrligprod )
       !
       ! ------------------------------------------------------------------
@@ -68 +72 @@
                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
+               zrligprod(ji,jj,jk) = zlgwp
             END DO
          END DO
@@ -90 +97 @@
          END DO
       END DO
+
+      !  Output of some diagnostics variables
+      !     ---------------------------------
+      IF( lk_iomput ) THEN
+         IF( knt == nrdttrc ) THEN
+            IF( iom_use( "LIGREM" ) )   CALL iom_put("LIGREM"   , zligrem(:,:,:) * 1e9 * 1.e+3 * rfact2r * tmask(:,:,:) )
+            IF( iom_use( "LIGPR" ) )    CALL iom_put("LIGPR"    , zligpr(:,:,:) * 1e9 * 1.e+3 * rfact2r * tmask(:,:,:) )
+            IF( iom_use( "LIGPRODR" ) ) CALL iom_put("LPRODR"   , zrligprod(:,:,:) * 1e9 * 1.e+3 * rfact2r * tmask(:,:,:) )
+         ENDIF
+      ENDIF
+      !
+      CALL wrk_dealloc( jpi, jpj, jpk, zligrem, zligpr, zrligprod )
 
       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)

branches/CNRS/dev_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90

@@ -20 +20 @@
 
    PUBLIC p4z_lim    
+   PUBLIC p5z_lim
    PUBLIC p4z_lim_init    
+   PUBLIC p5z_lim_init
    PUBLIC p4z_lim_alloc
 
    !! * Shared module variables
    REAL(wp), PUBLIC ::  concnno3    !:  NO3, PO4 half saturation   
+   REAL(wp), PUBLIC ::  concpno3    !:  NO3, PO4 half saturation   
    REAL(wp), PUBLIC ::  concdno3    !:  Phosphate half saturation for diatoms  
    REAL(wp), PUBLIC ::  concnnh4    !:  NH4 half saturation for phyto  
+   REAL(wp), PUBLIC ::  concpnh4    !:  NH4 half saturation for phyto  
    REAL(wp), PUBLIC ::  concdnh4    !:  NH4 half saturation for diatoms
+   REAL(wp), PUBLIC ::  concnpo4    !:  NH4 half saturation for diatoms
+   REAL(wp), PUBLIC ::  concppo4    !:  NH4 half saturation for diatoms
+   REAL(wp), PUBLIC ::  concdpo4    !:  NH4 half saturation for diatoms
    REAL(wp), PUBLIC ::  concnfer    !:  Iron half saturation for nanophyto 
+   REAL(wp), PUBLIC ::  concpfer    !:  Iron half saturation for nanophyto 
    REAL(wp), PUBLIC ::  concdfer    !:  Iron half saturation for diatoms  
    REAL(wp), PUBLIC ::  concbno3    !:  NO3 half saturation  for bacteria 
    REAL(wp), PUBLIC ::  concbnh4    !:  NH4 half saturation for bacteria
+   REAL(wp), PUBLIC ::  concbpo4    !:  PO4 half saturation for bacteria
    REAL(wp), PUBLIC ::  xsizedia    !:  Minimum size criteria for diatoms
+   REAL(wp), PUBLIC ::  xsizepic    !:  Minimum size criteria for diatoms
    REAL(wp), PUBLIC ::  xsizephy    !:  Minimum size criteria for nanophyto
    REAL(wp), PUBLIC ::  xsizern     !:  Size ratio for nanophytoplankton
+   REAL(wp), PUBLIC ::  xsizerp     !:  Size ratio for nanophytoplankton
    REAL(wp), PUBLIC ::  xsizerd     !:  Size ratio for diatoms
    REAL(wp), PUBLIC ::  xksi1       !:  half saturation constant for Si uptake 
@@ -41 +52 @@
    REAL(wp), PUBLIC ::  concbfe     !:  Fe half saturation for bacteria 
    REAL(wp), PUBLIC ::  oxymin      !:  half saturation constant for anoxia
+   REAL(wp), PUBLIC ::  qfnopt      !:  optimal Fe quota for nanophyto
+   REAL(wp), PUBLIC ::  qfpopt      !:  optimal Fe quota for nanophyto
+   REAL(wp), PUBLIC ::  qfdopt      !:  optimal Fe quota for diatoms
+   REAL(wp), PUBLIC ::  qnnmin      !:  optimal Fe quota for diatoms
+   REAL(wp), PUBLIC ::  qnnmax      !:  optimal Fe quota for diatoms
+   REAL(wp), PUBLIC ::  qpnmin      !:  optimal Fe quota for diatoms
+   REAL(wp), PUBLIC ::  qpnmax      !:  optimal Fe quota for diatoms
+   REAL(wp), PUBLIC ::  qnpmin      !:  optimal Fe quota for diatoms
+   REAL(wp), PUBLIC ::  qnpmax      !:  optimal Fe quota for diatoms
+   REAL(wp), PUBLIC ::  qppmin      !:  optimal Fe quota for diatoms
+   REAL(wp), PUBLIC ::  qppmax      !:  optimal Fe quota for diatoms
+   REAL(wp), PUBLIC ::  qndmin      !:  optimal Fe quota for diatoms
+   REAL(wp), PUBLIC ::  qndmax      !:  optimal Fe quota for diatoms
+   REAL(wp), PUBLIC ::  qpdmin      !:  optimal Fe quota for diatoms
+   REAL(wp), PUBLIC ::  qpdmax      !:  optimal Fe quota for diatoms
+   REAL(wp), PUBLIC ::  qfnmax      !:  optimal Fe quota for diatoms
+   REAL(wp), PUBLIC ::  qfpmax      !:  optimal Fe quota for diatoms
+   REAL(wp), PUBLIC ::  qfdmax      !:  optimal Fe quota for diatoms
    REAL(wp), PUBLIC ::  qnfelim     !:  optimal Fe quota for nanophyto
    REAL(wp), PUBLIC ::  qdfelim     !:  optimal Fe quota for diatoms
    REAL(wp), PUBLIC ::  caco3r      !:  mean rainratio 
+   REAL(wp), PUBLIC ::  zpsinh4
+   REAL(wp), PUBLIC ::  zpsino3
+   REAL(wp), PUBLIC ::  zpsiuptk
 
    !!* Phytoplankton limitation terms
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xnanono3   !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xpicono3   !: ???
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xdiatno3   !: ???
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xnanonh4   !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xpiconh4   !: ???
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xdiatnh4   !: ???
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xnanopo4   !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xpicopo4   !: ???
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xdiatpo4   !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xnanodop   !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xpicodop   !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xdiatdop   !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xnanofer   !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xpicofer   !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xdiatfer   !: ???
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimphy    !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimpic    !: ???
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimdia    !: ???
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimnfe    !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimpfe    !: ???
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimdfe    !: ???
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimsi     !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimsi2    !: ???
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimbac    !: ??
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimbacl   !: ??
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   concdfe    !: ???
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   concnfe    !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   fvnuptk
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   fvpuptk
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   fvduptk
+
+   !!*  Allometric variations of the quotas
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqnnmin    !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqnnmax    !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqpnmin    !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqpnmax    !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqnpmin    !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqnpmax    !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqppmin    !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqppmax    !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqndmin    !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqndmax    !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqpdmin    !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE,   DIMENSION(:,:,:)  ::   xqpdmax    !: ?
 
    ! Coefficient for iron limitation
-   REAL(wp) ::  xcoef1   = 0.0016  / 55.85  
+   REAL(wp) ::  xcoef1   = 0.00167  / 55.85  
    REAL(wp) ::  xcoef2   = 1.21E-5 * 14. / 55.85 / 7.625 * 0.5 * 1.5
    REAL(wp) ::  xcoef3   = 1.15E-4 * 14. / 55.85 / 7.625 * 0.5 
@@ -227 +288 @@
       !
    END SUBROUTINE p4z_lim
+
+
+   SUBROUTINE p5z_lim( kt, knt )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE p5z_lim  ***
+      !!
+      !! ** Purpose :   Compute the co-limitations by the various nutrients
+      !!              for the various phytoplankton species
+      !!
+      !! ** Method  : - ???
+      !!---------------------------------------------------------------------
+      !
+      INTEGER, INTENT(in)  :: kt, knt
+      !
+      INTEGER  ::   ji, jj, jk
+      REAL(wp) ::   zlim1, zlim2, zlim3, zlim4, zno3, zferlim
+      REAL(wp) ::   z1_trndia, z1_trnpic, z1_trnphy, ztem1, ztem2, zetot1
+      REAL(wp) ::   zratio, zration, zratiof, znutlim, zfalim
+      REAL(wp) ::   zconc1d, zconc1dnh4, zconc0n, zconc0nnh4, zconc0npo4, zconc0dpo4
+      REAL(wp) ::   zconc0p, zconc0pnh4, zconc0ppo4, zconcpfe, zconcnfe, zconcdfe
+      REAL(wp) ::   fanano, fananop, fananof, fadiat, fadiatp, fadiatf
+      REAL(wp) ::   fapico, fapicop, fapicof
+      REAL(wp) ::   zrpho, zrass, zcoef, zfuptk, zratchl
+      REAL(wp) ::   zfvn, zfvp, zfvf, zsizen, zsizep, zsized, znanochl, zpicochl, zdiatchl
+      REAL(wp) ::   zqfemn, zqfemp, zqfemd, zbactno3, zbactnh4
+      !!---------------------------------------------------------------------
+      !
+      !!---------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('p5z_lim')
+      !
+      zratchl = 6.0
+      !
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ! 
+               ! Tuning of the iron concentration to a minimum level that is set to the detection limit
+               !-------------------------------------
+               zno3    = trb(ji,jj,jk,jpno3) / 40.e-6
+               zferlim = MAX( 3e-11 * zno3 * zno3, 5e-12 )
+               zferlim = MIN( zferlim, 7e-11 )
+               trb(ji,jj,jk,jpfer) = MAX( trb(ji,jj,jk,jpfer), zferlim )
+
+               ! Computation of the mean relative size of each community
+               ! -------------------------------------------------------
+               z1_trnphy   = 1. / ( trb(ji,jj,jk,jpphy) + rtrn )
+               z1_trnpic   = 1. / ( trb(ji,jj,jk,jppic) + rtrn )
+               z1_trndia   = 1. / ( trb(ji,jj,jk,jpdia) + rtrn )
+               znanochl = trb(ji,jj,jk,jpnch) * z1_trnphy
+               zpicochl = trb(ji,jj,jk,jppch) * z1_trnpic
+               zdiatchl = trb(ji,jj,jk,jpdch) * z1_trndia
+
+               ! Computation of a variable Ks for iron on diatoms taking into account
+               ! that increasing biomass is made of generally bigger cells
+               !------------------------------------------------
+               zsized            = sized(ji,jj,jk)**0.81
+               zconcdfe          = concdfer * zsized
+               zconc1d           = concdno3 * zsized
+               zconc1dnh4        = concdnh4 * zsized
+               zconc0dpo4        = concdpo4 * zsized
+
+               zsizep            = 1.
+               zconcpfe          = concpfer * zsizep
+               zconc0p           = concpno3 * zsizep
+               zconc0pnh4        = concpnh4 * zsizep
+               zconc0ppo4        = concppo4 * zsizep
+
+               zsizen            = 1.
+               zconcnfe          = concnfer * zsizen
+               zconc0n           = concnno3 * zsizen
+               zconc0nnh4        = concnnh4 * zsizen
+               zconc0npo4        = concnpo4 * zsizen
+
+               ! Allometric variations of the minimum and maximum quotas
+               ! From Talmy et al. (2014) and Maranon et al. (2013)
+               ! -------------------------------------------------------
+               xqnnmin(ji,jj,jk) = qnnmin
+               xqnnmax(ji,jj,jk) = qnnmax
+               xqndmin(ji,jj,jk) = qndmin * sized(ji,jj,jk)**(-0.27)
+               xqndmax(ji,jj,jk) = qndmax
+               xqnpmin(ji,jj,jk) = qnpmin
+               xqnpmax(ji,jj,jk) = qnpmax
+
+               ! Computation of the optimal allocation parameters
+               ! Based on the different papers by Pahlow et al., and Smith et al.
+               ! -----------------------------------------------------------------
+               znutlim = MAX( trb(ji,jj,jk,jpnh4) / zconc0nnh4,    &
+                 &         trb(ji,jj,jk,jpno3) / zconc0n)
+               fanano = MAX(0.01, MIN(0.99, 1. / ( SQRT(znutlim) + 1.) ) )
+               znutlim = trb(ji,jj,jk,jppo4) / zconc0npo4
+               fananop = MAX(0.01, MIN(0.99, 1. / ( SQRT(znutlim) + 1.) ) )
+               znutlim = biron(ji,jj,jk) / zconcnfe
+               fananof = MAX(0.01, MIN(0.99, 1. / ( SQRT(znutlim) + 1.) ) )
+               znutlim = MAX( trb(ji,jj,jk,jpnh4) / zconc0pnh4,    &
+                 &         trb(ji,jj,jk,jpno3) / zconc0p)
+               fapico = MAX(0.01, MIN(0.99, 1. / ( SQRT(znutlim) + 1.) ) )
+               znutlim = trb(ji,jj,jk,jppo4) / zconc0ppo4
+               fapicop = MAX(0.01, MIN(0.99, 1. / ( SQRT(znutlim) + 1.) ) )
+               znutlim = biron(ji,jj,jk) / zconcpfe
+               fapicof = MAX(0.01, MIN(0.99, 1. / ( SQRT(znutlim) + 1.) ) )
+               znutlim = MAX( trb(ji,jj,jk,jpnh4) / zconc1dnh4,    &
+                 &         trb(ji,jj,jk,jpno3) / zconc1d )
+               fadiat = MAX(0.01, MIN(0.99, 1. / ( SQRT(znutlim) + 1.) ) )
+               znutlim = trb(ji,jj,jk,jppo4) / zconc0dpo4
+               fadiatp = MAX(0.01, MIN(0.99, 1. / ( SQRT(znutlim) + 1.) ) )
+               znutlim = biron(ji,jj,jk) / zconcdfe
+               fadiatf = MAX(0.01, MIN(0.99, 1. / ( SQRT(znutlim) + 1.) ) )
+               !
+               ! Michaelis-Menten Limitation term for nutrients Small bacteria
+               ! -------------------------------------------------------------
+               zbactnh4 = trb(ji,jj,jk,jpnh4) / ( concbnh4 + trb(ji,jj,jk,jpnh4) )
+               zbactno3 = trb(ji,jj,jk,jpno3) / ( concbno3 + trb(ji,jj,jk,jpno3) ) * (1. - zbactnh4)
+               !
+               zlim1    = zbactno3 + zbactnh4
+               zlim2    = trb(ji,jj,jk,jppo4) / ( trb(ji,jj,jk,jppo4) + concbpo4)
+               zlim3    = biron(ji,jj,jk) / ( concbfe + biron(ji,jj,jk) )
+               zlim4    = trb(ji,jj,jk,jpdoc) / ( xkdoc   + trb(ji,jj,jk,jpdoc) )
+               xlimbacl(ji,jj,jk) = MIN( zlim1, zlim2, zlim3 )
+               xlimbac (ji,jj,jk) = xlimbacl(ji,jj,jk) * zlim4
+               !
+               ! Michaelis-Menten Limitation term for nutrients Small flagellates
+               ! -----------------------------------------------
+               zfalim = (1.-fanano) / fanano
+               xnanonh4(ji,jj,jk) = (1. - fanano) * trb(ji,jj,jk,jpnh4) / ( zfalim * zconc0nnh4 + trb(ji,jj,jk,jpnh4) )
+               xnanono3(ji,jj,jk) = (1. - fanano) * trb(ji,jj,jk,jpno3) / ( zfalim * zconc0n + trb(ji,jj,jk,jpno3) )  &
+               &                    * (1. - xnanonh4(ji,jj,jk))
+               !
+               zfalim = (1.-fananop) / fananop
+               xnanopo4(ji,jj,jk) = (1. - fananop) * trb(ji,jj,jk,jppo4) / ( trb(ji,jj,jk,jppo4) + zfalim * zconc0npo4 )
+               xnanodop(ji,jj,jk) = trb(ji,jj,jk,jpdop) / ( trb(ji,jj,jk,jpdop) + xkdoc )   &
+               &                    * ( 1.0 - xnanopo4(ji,jj,jk) )
+               xnanodop(ji,jj,jk) = 0.
+               !
+               zfalim = (1.-fananof) / fananof
+               xnanofer(ji,jj,jk) = (1. - fananof) * biron(ji,jj,jk) / ( biron(ji,jj,jk) + zfalim * zconcnfe )
+               !
+               zratiof   = trb(ji,jj,jk,jpnfe) * z1_trnphy
+               zqfemn = xcoef1 * znanochl + xcoef2 + xcoef3 * xnanono3(ji,jj,jk)
+               !
+               zration = trb(ji,jj,jk,jpnph) * z1_trnphy
+               zration = MIN(xqnnmax(ji,jj,jk), MAX( 2. * xqnnmin(ji,jj,jk), zration ))
+               fvnuptk(ji,jj,jk) = 1. / zpsiuptk * rno3 * 2. * xqnnmin(ji,jj,jk) / (zration + rtrn)  &
+               &                   * MAX(0., (1. - zratchl * znanochl / 12. ) )
+               !
+               zlim1    = max(0., (zration - 2. * xqnnmin(ji,jj,jk) )  &
+               &          / (xqnnmax(ji,jj,jk) - 2. * xqnnmin(ji,jj,jk) ) ) * xqnnmax(ji,jj,jk)  &
+               &          / (zration + rtrn)
+               zlim3    = MAX( 0.,( zratiof - zqfemn ) / qfnopt )
+               xlimnfe(ji,jj,jk) = MIN( 1., zlim3 )
+               xlimphy(ji,jj,jk) = MIN( 1., zlim1, zlim3 )
+               !
+               ! Michaelis-Menten Limitation term for nutrients picophytoplankton
+               ! ----------------------------------------------------------------
+               zfalim = (1.-fapico) / fapico
+               xpiconh4(ji,jj,jk) = (1. - fapico) * trb(ji,jj,jk,jpnh4) / ( zfalim * zconc0pnh4 + trb(ji,jj,jk,jpnh4) )
+               xpicono3(ji,jj,jk) = (1. - fapico) * trb(ji,jj,jk,jpno3) / ( zfalim * zconc0p + trb(ji,jj,jk,jpno3) )  &
+               &                    * (1. - xpiconh4(ji,jj,jk))
+               !
+               zfalim = (1.-fapicop) / fapicop
+               xpicopo4(ji,jj,jk) = (1. - fapicop) * trb(ji,jj,jk,jppo4) / ( trb(ji,jj,jk,jppo4) + zfalim * zconc0ppo4 )
+               xpicodop(ji,jj,jk) = trb(ji,jj,jk,jpdop) / ( trb(ji,jj,jk,jpdop) + xkdoc )   &
+               &                    * ( 1.0 - xpicopo4(ji,jj,jk) )
+               xpicodop(ji,jj,jk) = 0.
+               !
+               zfalim = (1.-fapicof) / fapicof
+               xpicofer(ji,jj,jk) = (1. - fapicof) * biron(ji,jj,jk) / ( biron(ji,jj,jk) + zfalim * zconcpfe )
+               !
+               zratiof   = trb(ji,jj,jk,jppfe) * z1_trnpic
+               zqfemp = xcoef1 * zpicochl + xcoef2 + xcoef3 * xpicono3(ji,jj,jk)
+               !
+               zration   = trb(ji,jj,jk,jpnpi) * z1_trnpic
+               zration = MIN(xqnpmax(ji,jj,jk), MAX( 2. * xqnpmin(ji,jj,jk), zration ))
+               fvpuptk(ji,jj,jk) = 1. / zpsiuptk * rno3 * 2. * xqnpmin(ji,jj,jk) / (zration + rtrn)  &
+               &                   * MAX(0., (1. - zratchl * zpicochl / 12. ) )
+               !
+               zlim1    = max(0., (zration - 2. * xqnpmin(ji,jj,jk) )  &
+               &          / (xqnpmax(ji,jj,jk) - 2. * xqnpmin(ji,jj,jk) ) ) * xqnpmax(ji,jj,jk)  &
+               &          / (zration + rtrn)
+               zlim3    = MAX( 0.,( zratiof - zqfemp ) / qfpopt )
+               xlimpfe(ji,jj,jk) = MIN( 1., zlim3 )
+               xlimpic(ji,jj,jk) = MIN( 1., zlim1, zlim3 )
+               !
+               !   Michaelis-Menten Limitation term for nutrients Diatoms
+               !   ------------------------------------------------------
+               zfalim = (1.-fadiat) / fadiat
+               xdiatnh4(ji,jj,jk) = (1. - fadiat) * trb(ji,jj,jk,jpnh4) / ( zfalim * zconc1dnh4 + trb(ji,jj,jk,jpnh4) )
+               xdiatno3(ji,jj,jk) = (1. - fadiat) * trb(ji,jj,jk,jpno3) / ( zfalim * zconc1d + trb(ji,jj,jk,jpno3) )  &
+               &                    * (1. - xdiatnh4(ji,jj,jk))
+               !
+               zfalim = (1.-fadiatp) / fadiatp
+               xdiatpo4(ji,jj,jk) = (1. - fadiatp) * trb(ji,jj,jk,jppo4) / ( trb(ji,jj,jk,jppo4) + zfalim * zconc0dpo4 )
+               xdiatdop(ji,jj,jk) = trb(ji,jj,jk,jpdop) / ( trb(ji,jj,jk,jpdop) + xkdoc )  &
+               &                    * ( 1.0 - xdiatpo4(ji,jj,jk) )
+               xdiatdop(ji,jj,jk) = 0.
+               !
+               zfalim = (1.-fadiatf) / fadiatf
+               xdiatfer(ji,jj,jk) = (1. - fadiatf) * biron(ji,jj,jk) / ( biron(ji,jj,jk) + zfalim * zconcdfe )
+               !
+               zratiof   = trb(ji,jj,jk,jpdfe) * z1_trndia
+               zqfemd = xcoef1 * zdiatchl + xcoef2 + xcoef3 * xdiatno3(ji,jj,jk)
+               !
+               zration   = trb(ji,jj,jk,jpndi) * z1_trndia
+               zration = MIN(xqndmax(ji,jj,jk), MAX( 2. * xqndmin(ji,jj,jk), zration ))
+               fvduptk(ji,jj,jk) = 1. / zpsiuptk * rno3 * 2. * xqndmin(ji,jj,jk) / (zration + rtrn)   &
+               &                   * MAX(0., (1. - zratchl * zdiatchl / 12. ) )
+               !
+               zlim1    = max(0., (zration - 2. * xqndmin(ji,jj,jk) )    &
+               &          / (xqndmax(ji,jj,jk) - 2. * xqndmin(ji,jj,jk) ) )   &
+               &          * xqndmax(ji,jj,jk) / (zration + rtrn)
+               zlim3    = trb(ji,jj,jk,jpsil) / ( trb(ji,jj,jk,jpsil) + xksi(ji,jj) )
+               zlim4    = MAX( 0., ( zratiof - zqfemd ) / qfdopt )
+               xlimdfe(ji,jj,jk) = MIN( 1., zlim4 )
+               xlimdia(ji,jj,jk) = MIN( 1., zlim1, zlim3, zlim4 )
+               xlimsi(ji,jj,jk)  = MIN( 1., zlim3, zlim4 )
+               xlimsi2(ji,jj,jk)  = MIN( 1., zlim1, zlim4 )
+            END DO
+         END DO
+      END DO
+      !
+      ! Compute the phosphorus quota values. It is based on Litchmann et al., 2004 and Daines et al, 2013.
+      ! The relative contribution of three fonctional pools are computed: light harvesting apparatus, 
+      ! nutrient uptake pool and assembly machinery. DNA is assumed to represent 1% of the dry mass of 
+      ! phytoplankton (see Daines et al., 2013). 
+      ! --------------------------------------------------------------------------------------------------
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ! Size estimation of nanophytoplankton
+               ! ------------------------------------
+               zfvn = zpsiuptk / rno3 * fvnuptk(ji,jj,jk)
+               sizen(ji,jj,jk) = MAX(1., MIN(xsizern, 1.0 / ( MAX(rtrn, zfvn) ) ) )
+
+               ! N/P ratio of nanophytoplankton
+               ! ------------------------------
+               zfuptk = 0.23 * zfvn
+               zrpho = 2.24 * trb(ji,jj,jk,jpnch) / ( trb(ji,jj,jk,jpnph) * rno3 * 15. + rtrn )
+               zrass = 1. - 0.2 - zrpho - zfuptk
+               xqpnmax(ji,jj,jk) = ( zfuptk + zrpho ) * 0.0128 * 16. + zrass * 1./ 7.2 * 16. + 0.13
+               xqpnmax(ji,jj,jk) = xqpnmax(ji,jj,jk) * trb(ji,jj,jk,jpnph) / ( trb(ji,jj,jk,jpphy) + rtrn ) 
+               xqpnmin(ji,jj,jk) = qpnmin
+
+               ! Size estimation of picophytoplankton
+               ! ------------------------------------
+               zfvn = zpsiuptk / rno3  * fvpuptk(ji,jj,jk)
+               sizep(ji,jj,jk) = MAX(1., MIN(xsizerp, 1.0 / ( MAX(rtrn, zfvn) ) ) )
+
+               ! N/P ratio of picophytoplankton
+               ! ------------------------------
+               zfuptk = 0.35 * zfvn
+               zrpho = 2.24 * trb(ji,jj,jk,jppch) / ( trb(ji,jj,jk,jpnpi) * rno3 * 15. + rtrn )
+               zrass = 1. - 0.4 - zrpho - zfuptk
+               xqppmax(ji,jj,jk) =  (zrpho + zfuptk) * 0.0128 * 16. + zrass * 1./ 10.9 * 16. + 0.13
+               xqppmax(ji,jj,jk) = xqppmax(ji,jj,jk) * trb(ji,jj,jk,jpnpi) / ( trb(ji,jj,jk,jppic) + rtrn ) 
+               xqppmin(ji,jj,jk) = qppmin
+
+               ! Size estimation of diatoms
+               ! --------------------------
+               zfvn = zpsiuptk / rno3 * fvduptk(ji,jj,jk)
+               sized(ji,jj,jk) = MAX(1., MIN(xsizerd, 1.0 / ( MAX(rtrn, zfvn) ) ) )
+!               zcoef = trb(ji,jj,jk,jpdia) - MIN(xsizedia, trb(ji,jj,jk,jpdia) )
+!               sized(ji,jj,jk) = 1. + xsizerd * zcoef *1E6 / ( 1. + zcoef * 1E6 )
+
+               ! N/P ratio of diatoms
+               ! --------------------
+               zfuptk = 0.2 * zfvn
+               zrpho = 2.24 * trb(ji,jj,jk,jpdch) / ( trb(ji,jj,jk,jpndi) * rno3 * 15. + rtrn )
+               zrass = 1. - 0.2 - zrpho - zfuptk
+               xqpdmax(ji,jj,jk) = ( zfuptk + zrpho ) * 0.0128 * 16. + zrass * 1./ 7.2 * 16. + 0.13
+               xqpdmax(ji,jj,jk) = xqpdmax(ji,jj,jk) * trb(ji,jj,jk,jpndi) / ( trb(ji,jj,jk,jpdia) + rtrn )
+               xqpdmin(ji,jj,jk) = qpdmin * sized(ji,jj,jk)**(-0.27)
+
+            END DO
+         END DO
+      END DO
+
+      ! Compute the fraction of nanophytoplankton that is made of calcifiers
+      ! --------------------------------------------------------------------
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zlim1 =  trb(ji,jj,jk,jpnh4) / ( trb(ji,jj,jk,jpnh4) + concnnh4 ) + trb(ji,jj,jk,jpno3)    &
+               &        / ( trb(ji,jj,jk,jpno3) + concnno3 ) * ( 1.0 - trb(ji,jj,jk,jpnh4)   &
+               &        / ( trb(ji,jj,jk,jpnh4) + concnnh4 ) )
+               zlim2  = trb(ji,jj,jk,jppo4) / ( trb(ji,jj,jk,jppo4) + concnpo4 )
+               zlim3  = trb(ji,jj,jk,jpfer) / ( trb(ji,jj,jk,jpfer) +  1.E-10 )
+               ztem1  = MAX( 0., tsn(ji,jj,jk,jp_tem) )
+               ztem2  = tsn(ji,jj,jk,jp_tem) - 10.
+               zetot1 = MAX( 0., etot(ji,jj,jk) - 1.) / ( 4. + etot(ji,jj,jk) ) * 20. / ( 20. + etot(ji,jj,jk) )
+
+!               xfracal(ji,jj,jk) = caco3r * MIN( zlim1, zlim2, zlim3 )                  &
+               xfracal(ji,jj,jk) = caco3r * xlimphy(ji,jj,jk)                &
+               &                   * ztem1 / ( 1. + ztem1 ) * MAX( 1., trb(ji,jj,jk,jpphy)*1E6 )   &
+                  &                * ( 1. + EXP(-ztem2 * ztem2 / 25. ) )         &
+                  &                * zetot1 * MIN( 1., 50. / ( hmld(ji,jj) + rtrn ) )
+               xfracal(ji,jj,jk) = MAX( 0.02, MIN( 0.8 , xfracal(ji,jj,jk) ) )
+            END DO
+         END DO
+      END DO
+      !
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ! denitrification factor computed from O2 levels
+               nitrfac(ji,jj,jk) = MAX(  0.e0, 0.4 * ( 6.e-6  - trb(ji,jj,jk,jpoxy) )    &
+                  &                                / ( oxymin + trb(ji,jj,jk,jpoxy) )  )
+               nitrfac(ji,jj,jk) = MIN( 1., nitrfac(ji,jj,jk) )
+            END DO
+         END DO
+      END DO
+      !
+      IF( lk_iomput .AND. knt == nrdttrc ) THEN        ! save output diagnostics
+        IF( iom_use( "xfracal" ) ) CALL iom_put( "xfracal", xfracal(:,:,:) * tmask(:,:,:) )  ! euphotic layer deptht
+        IF( iom_use( "LNnut"   ) ) CALL iom_put( "LNnut"  , xlimphy(:,:,:) * tmask(:,:,:) )  ! Nutrient limitation term
+        IF( iom_use( "LPnut"   ) ) CALL iom_put( "LPnut"  , xlimpic(:,:,:) * tmask(:,:,:) )  ! Nutrient limitation term
+        IF( iom_use( "LDnut"   ) ) CALL iom_put( "LDnut"  , xlimdia(:,:,:) * tmask(:,:,:) )  ! Nutrient limitation term
+        IF( iom_use( "LNFe"    ) ) CALL iom_put( "LNFe"   , xlimnfe(:,:,:) * tmask(:,:,:) )  ! Iron limitation term
+        IF( iom_use( "LPFe"    ) ) CALL iom_put( "LPFe"   , xlimpfe(:,:,:) * tmask(:,:,:) )  ! Iron limitation term
+        IF( iom_use( "LDFe"    ) ) CALL iom_put( "LDFe"   , xlimdfe(:,:,:) * tmask(:,:,:) )  ! Iron limitation term
+        IF( iom_use( "SIZEN"   ) ) CALL iom_put( "SIZEN"  , sizen(:,:,:) * tmask(:,:,:) )  ! Iron limitation term
+        IF( iom_use( "SIZEP"   ) ) CALL iom_put( "SIZEP"  , sizep(:,:,:) * tmask(:,:,:) )  ! Iron limitation term
+        IF( iom_use( "SIZED"   ) ) CALL iom_put( "SIZED"  , sized(:,:,:) * tmask(:,:,:) )  ! Iron limitation term
+      ENDIF
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('p5z_lim')
+      !
+   END SUBROUTINE p5z_lim
 
    SUBROUTINE p4z_lim_init
@@ -286 +674 @@
    END SUBROUTINE p4z_lim_init
 
+   SUBROUTINE p5z_lim_init
+
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE p5z_lim_init  ***
+      !!
+      !! ** Purpose :   Initialization of nutrient limitation parameters
+      !!
+      !! ** Method  :   Read the nampislim and nampisquota namelists and check
+      !!      the parameters called at the first timestep (nittrc000)
+      !!
+      !! ** input   :   Namelist nampislim
+      !!
+      !!----------------------------------------------------------------------
+
+      NAMELIST/namp5zlim/ concnno3, concpno3, concdno3, concnnh4, concpnh4, concdnh4,  &
+         &                concnfer, concpfer, concdfer, concbfe, concnpo4, concppo4,   &
+         &                concdpo4, concbno3, concbnh4, concbpo4, xsizedia, xsizepic,  &
+         &                xsizephy, xsizern, xsizerp, xsizerd, xksi1, xksi2, xkdoc,    &
+         &                caco3r, oxymin
+
+      NAMELIST/namp5zquota/ qnnmin, qnnmax, qpnmin, qpnmax, qnpmin, qnpmax, qppmin,      &
+         &                  qppmax, qndmin, qndmax, qpdmin, qpdmax, qfnmax, qfpmax, qfdmax,  &
+         &                  qfnopt, qfpopt, qfdopt
+
+
+      INTEGER :: ios                 ! Local integer output status for namelist read
+
+      REWIND( numnatp_ref )              ! Namelist nampislim in reference namelist : Pisces nutrient limitation parameters
+      READ  ( numnatp_ref, namp5zlim, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampislim in reference namelist', lwp )
+
+      REWIND( numnatp_cfg )              ! Namelist nampislim in configuration namelist : Pisces nutrient limitation parameters 
+      READ  ( numnatp_cfg, namp5zlim, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampislim in configuration namelist', lwp )
+      IF(lwm) WRITE ( numonp, namp5zlim )
+
+      IF(lwp) THEN                         ! control print
+         WRITE(numout,*) ' '
+         WRITE(numout,*) ' Namelist parameters for nutrient limitations, namp5zlim'
+         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
+         WRITE(numout,*) '    mean rainratio                           caco3r    = ', caco3r
+         WRITE(numout,*) '    NO3 half saturation of nanophyto         concnno3  = ', concnno3
+         WRITE(numout,*) '    NO3 half saturation of picophyto         concpno3  = ', concpno3
+         WRITE(numout,*) '    NO3 half saturation of diatoms           concdno3  = ', concdno3
+         WRITE(numout,*) '    NH4 half saturation for phyto            concnnh4  = ', concnnh4
+         WRITE(numout,*) '    NH4 half saturation for pico             concpnh4  = ', concpnh4
+         WRITE(numout,*) '    NH4 half saturation for diatoms          concdnh4  = ', concdnh4
+         WRITE(numout,*) '    PO4 half saturation for phyto            concnpo4  = ', concnpo4
+         WRITE(numout,*) '    PO4 half saturation for pico             concppo4  = ', concppo4
+         WRITE(numout,*) '    PO4 half saturation for diatoms          concdpo4  = ', concdpo4
+         WRITE(numout,*) '    half saturation constant for Si uptake   xksi1     = ', xksi1
+         WRITE(numout,*) '    half saturation constant for Si/C        xksi2     = ', xksi2
+         WRITE(numout,*) '    half-sat. of DOC remineralization        xkdoc     = ', xkdoc
+         WRITE(numout,*) '    Iron half saturation for nanophyto       concnfer  = ', concnfer
+         WRITE(numout,*) '    Iron half saturation for picophyto       concpfer  = ', concpfer
+         WRITE(numout,*) '    Iron half saturation for diatoms         concdfer  = ', concdfer
+         WRITE(numout,*) '    size ratio for nanophytoplankton         xsizern   = ', xsizern
+         WRITE(numout,*) '    size ratio for picophytoplankton         xsizerp   = ', xsizerp
+         WRITE(numout,*) '    size ratio for diatoms                   xsizerd   = ', xsizerd
+         WRITE(numout,*) '    NO3 half saturation of bacteria          concbno3  = ', concbno3
+         WRITE(numout,*) '    NH4 half saturation for bacteria         concbnh4  = ', concbnh4
+         WRITE(numout,*) '    Minimum size criteria for diatoms        xsizedia  = ', xsizedia
+         WRITE(numout,*) '    Minimum size criteria for picophyto      xsizepic  = ', xsizepic
+         WRITE(numout,*) '    Minimum size criteria for nanophyto      xsizephy  = ', xsizephy
+         WRITE(numout,*) '    Fe half saturation for bacteria          concbfe   = ', concbfe
+         WRITE(numout,*) '    halk saturation constant for anoxia       oxymin   =' , oxymin
+      ENDIF
+
+      REWIND( numnatp_ref )              ! Namelist nampislim in reference namelist : Pisces nutrient limitation parameters
+      READ  ( numnatp_ref, namp5zquota, IOSTAT = ios, ERR = 903)
+903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisquota in reference namelist', lwp )
+
+      REWIND( numnatp_cfg )              ! Namelist nampislim in configuration namelist : Pisces nutrient limitation parameters 
+      READ  ( numnatp_cfg, namp5zquota, IOSTAT = ios, ERR = 904 )
+904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisquota in configuration namelist', lwp )
+      IF(lwm) WRITE ( numonp, namp5zquota )
+
+      IF(lwp) THEN                         ! control print
+         WRITE(numout,*) ' '
+         WRITE(numout,*) ' Namelist parameters for nutrient limitations, namp5zquota'
+         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
+         WRITE(numout,*) '    optimal Fe quota for nano.               qfnopt    = ', qfnopt
+         WRITE(numout,*) '    optimal Fe quota for pico.               qfpopt    = ', qfpopt
+         WRITE(numout,*) '    Optimal Fe quota for diatoms             qfdopt    = ', qfdopt
+         WRITE(numout,*) '    Minimal N quota for nano                 qnnmin    = ', qnnmin
+         WRITE(numout,*) '    Maximal N quota for nano                 qnnmax    = ', qnnmax
+         WRITE(numout,*) '    Minimal P quota for nano                 qpnmin    = ', qpnmin
+         WRITE(numout,*) '    Maximal P quota for nano                 qpnmax    = ', qpnmax
+         WRITE(numout,*) '    Minimal N quota for pico                 qnpmin    = ', qnpmin
+         WRITE(numout,*) '    Maximal N quota for pico                 qnpmax    = ', qnpmax
+         WRITE(numout,*) '    Minimal P quota for pico                 qppmin    = ', qppmin
+         WRITE(numout,*) '    Maximal P quota for pico                 qppmax    = ', qppmax
+         WRITE(numout,*) '    Minimal N quota for diatoms              qndmin    = ', qndmin
+         WRITE(numout,*) '    Maximal N quota for diatoms              qndmax    = ', qndmax
+         WRITE(numout,*) '    Minimal P quota for diatoms              qpdmin    = ', qpdmin
+         WRITE(numout,*) '    Maximal P quota for diatoms              qpdmax    = ', qpdmax
+         WRITE(numout,*) '    Maximal Fe quota for nanophyto.          qfnmax    = ', qfnmax
+         WRITE(numout,*) '    Maximal Fe quota for picophyto.          qfpmax    = ', qfpmax
+         WRITE(numout,*) '    Maximal Fe quota for diatoms             qfdmax    = ', qfdmax
+      ENDIF
+
+      zpsino3 = 2.3 * rno3
+      zpsinh4 = 1.8 * rno3
+      zpsiuptk = 2.3 * rno3
+      !
+      nitrfac (:,:,:) = 0._wp
+      !
+   END SUBROUTINE p5z_lim_init
+
    INTEGER FUNCTION p4z_lim_alloc()
       !!----------------------------------------------------------------------
@@ -291 +788 @@
       !!----------------------------------------------------------------------
       USE lib_mpp , ONLY: ctl_warn
+      INTEGER ::   ierr(3)        ! Local variables
       !!----------------------------------------------------------------------
+      ierr(:) = 0
 
       !*  Biological arrays for phytoplankton growth
@@ -301 +800 @@
          &      xlimbac (jpi,jpj,jpk), xlimbacl(jpi,jpj,jpk),       &
          &      concnfe (jpi,jpj,jpk), concdfe (jpi,jpj,jpk),       &
-         &      xlimsi  (jpi,jpj,jpk), STAT=p4z_lim_alloc )
+         &      xlimsi  (jpi,jpj,jpk), STAT=ierr(1) )
+
+      IF (ln_p5z) THEN
+      !*  Biological arrays for phytoplankton growth
+         ALLOCATE( xpicono3(jpi,jpj,jpk), xpiconh4(jpi,jpj,jpk),       &
+         &      xpicopo4(jpi,jpj,jpk), xpicodop(jpi,jpj,jpk),       &
+         &      xnanodop(jpi,jpj,jpk), xdiatdop(jpi,jpj,jpk),       &
+         &      xnanofer(jpi,jpj,jpk), xdiatfer(jpi,jpj,jpk),       &
+         &      xpicofer(jpi,jpj,jpk), xlimpfe (jpi,jpj,jpk),       &
+         &      fvnuptk (jpi,jpj,jpk), fvduptk (jpi,jpj,jpk),       &
+         &      fvpuptk (jpi,jpj,jpk), xlimpic (jpi,jpj,jpk),       &
+         &      xlimsi2 (jpi,jpj,jpk), STAT=ierr(2) )
+
+      !*  Minimum/maximum quotas of phytoplankton
+         ALLOCATE( xqnnmin (jpi,jpj,jpk), xqnnmax(jpi,jpj,jpk),       &
+         &      xqpnmin (jpi,jpj,jpk), xqpnmax(jpi,jpj,jpk),       &
+         &      xqnpmin (jpi,jpj,jpk), xqnpmax(jpi,jpj,jpk),       &
+         &      xqppmin (jpi,jpj,jpk), xqppmax(jpi,jpj,jpk),       &
+         &      xqndmin (jpi,jpj,jpk), xqndmax(jpi,jpj,jpk),       &
+         &      xqpdmin (jpi,jpj,jpk), xqpdmax(jpi,jpj,jpk),     STAT=ierr(3) )
+      ENDIF
+      !
+
+      p4z_lim_alloc = MAXVAL( ierr ) 
       !
       IF( p4z_lim_alloc /= 0 ) CALL ctl_warn('p4z_lim_alloc : failed to allocate arrays.')

branches/CNRS/dev_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90

@@ -131 +131 @@
          IF( iom_use( "CO3sat" ) ) CALL iom_put( "CO3sat", zco3sat(:,:,:) * 1.e+3            * tmask(:,:,:) )
          IF( iom_use( "DCAL"   ) ) CALL iom_put( "DCAL"  , zcaldiss(:,:,:) * 1.e+3 * rfact2r * tmask(:,:,:) )
+         IF( iom_use( "AOU"    ) )  CALL iom_put( "AOU"  , MAX(0., (chemo2(:,:,:) - trb(:,:,:,jpoxy) ) ) * 1.E6 * tmask(:,:,:) )
       ENDIF
       !

branches/CNRS/dev_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmeso.F90

@@ -72 +72 @@
       REAL(wp) :: zgrazfffp, zgrazfffg, zgrazffep, zgrazffeg
       CHARACTER (len=25) :: charout
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zgrazing, zw3d
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zgrazing, zw3d, zfezoo2
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zz2ligprod
 
       !!---------------------------------------------------------------------
@@ -78 +79 @@
       IF( nn_timing == 1 )  CALL timing_start('p4z_meso')
       !
-      CALL wrk_alloc( jpi, jpj, jpk, zgrazing )
+      CALL wrk_alloc( jpi, jpj, jpk, zgrazing, zfezoo2 )
       zgrazing(:,:,:) = 0._wp
+      zfezoo2 (:,:,:) = 0._wp
+      !
+      IF (ln_ligand) THEN
+         CALL wrk_alloc( jpi, jpj, jpk, zz2ligprod )
+         zz2ligprod(:,:,:) = 0._wp
+      ENDIF
 
       DO jk = 1, jpkm1
@@ -176 +183 @@
                tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zgrarem2 - zgrarsig
                !
-               IF( ln_ligand ) tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) + (zgrarem2 - zgrarsig) * ldocz
+               IF( ln_ligand ) THEN
+                  tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) + (zgrarem2 - zgrarsig) * ldocz
+                  zz2ligprod(ji,jj,jk) = (zgrarem2 - zgrarsig) * ldocz
+               ENDIF
                !
                tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) - o2ut * zgrarsig
                tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + zgrafer2
+               zfezoo2(ji,jj,jk) = zgrafer2
                tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) + zgrarsig
                tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * zgrarsig              
@@ -229 +240 @@
             CALL iom_put( "PCAL", zw3d )  
          ENDIF
+         IF( iom_use( "FEZOO2" ) ) THEN
+            zw3d(:,:,:) = zfezoo2(:,:,:) * 1e9 * 1.e+3 * rfact2r * tmask(:,:,:)   !
+            CALL iom_put( "FEZOO2", zw3d )
+         ENDIF
+         IF( iom_use( "LPRODZ2" ) .AND. ln_ligand )  THEN
+            zw3d(:,:,:) = zz2ligprod(:,:,:) * 1e9 * 1.e+3 * rfact2r * tmask(:,:,:)
+            CALL iom_put( "LPRODZ2"  , zw3d )
+         ENDIF
          CALL wrk_dealloc( jpi, jpj, jpk, zw3d )
       ENDIF
@@ -238 +257 @@
       ENDIF
       !
-      CALL wrk_dealloc( jpi, jpj, jpk, zgrazing )
+      CALL wrk_dealloc( jpi, jpj, jpk, zgrazing, zfezoo2 )
+      IF ( ln_ligand ) CALL wrk_dealloc( jpi, jpj, jpk, zz2ligprod )
       !
       IF( nn_timing == 1 )  CALL timing_stop('p4z_meso')

branches/CNRS/dev_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmicro.F90

@@ -71 +71 @@
       REAL(wp) :: zgrazp, zgrazm, zgrazsd
       REAL(wp) :: zgrazmf, zgrazsf, zgrazpf
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zgrazing, zw3d
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zgrazing, zw3d, zfezoo
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zzligprod
       CHARACTER (len=25) :: charout
       !!---------------------------------------------------------------------
@@ -77 +78 @@
       IF( nn_timing == 1 )  CALL timing_start('p4z_micro')
       !
-      CALL wrk_alloc( jpi, jpj, jpk, zgrazing )
+      CALL wrk_alloc( jpi, jpj, jpk, zgrazing, zfezoo )
+      IF (ln_ligand) THEN
+         CALL wrk_alloc( jpi, jpj, jpk, zzligprod )
+         zzligprod(:,:,:) = 0._wp
+      ENDIF
       !
       DO jk = 1, jpkm1
@@ -139 +144 @@
                tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zgrarem - zgrarsig
                !
-               IF( ln_ligand ) tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) + (zgrarem - zgrarsig) * ldocz
+               IF( ln_ligand ) THEN
+                  tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) + (zgrarem - zgrarsig) * ldocz
+                  zzligprod(ji,jj,jk) = (zgrarem - zgrarsig) * ldocz
+               ENDIF
                !
                tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) - o2ut * zgrarsig
                tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + zgrafer
+               zfezoo(ji,jj,jk) = zgrafer
                tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zgrapoc
                prodpoc(ji,jj,jk) = prodpoc(ji,jj,jk) + zgrapoc
@@ -184 +193 @@
               CALL iom_put( "GRAZ1", zw3d )
            ENDIF
+           IF( iom_use( "FEZOO" ) ) THEN
+              zw3d(:,:,:) = zfezoo(:,:,:) * 1e9 * 1.e+3 * rfact2r * tmask(:,:,:)   !
+              CALL iom_put( "FEZOO", zw3d )
+           ENDIF
+           IF( iom_use( "LPRODZ" ) .AND. ln_ligand )  THEN
+              zw3d(:,:,:) = zzligprod(:,:,:) * 1e9 * 1.e+3 * rfact2r * tmask(:,:,:)
+              CALL iom_put( "LPRODZ"  , zw3d )
+           ENDIF
            CALL wrk_dealloc( jpi, jpj, jpk, zw3d )
          ENDIF
@@ -194 +211 @@
       ENDIF
       !
-      CALL wrk_dealloc( jpi, jpj, jpk, zgrazing )
+      CALL wrk_dealloc( jpi, jpj, jpk, zgrazing, zfezoo )
+      IF (ln_ligand) CALL wrk_dealloc( jpi, jpj, jpk, zzligprod )
       !
       IF( nn_timing == 1 )  CALL timing_stop('p4z_micro')

branches/CNRS/dev_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zpoc.F90

@@ -67 +67 @@
       REAL(wp), POINTER, DIMENSION(:,:,:)   :: zremipoc, zremigoc, zorem3, ztremint
       REAL(wp), POINTER, DIMENSION(:,:,:,:) :: alphag
+      REAL(wp), POINTER, DIMENSION(:,:,:)   :: zfolimi
       !!---------------------------------------------------------------------
       !
@@ -74 +75 @@
       CALL wrk_alloc( jpi, jpj,      totprod,  totthick, totcons )
       CALL wrk_alloc( jpi, jpj, jpk, zremipoc, zremigoc, zorem3, ztremint )
+      CALL wrk_alloc( jpi, jpj, jpk, zfolimi )
       ALLOCATE( alphag(jpi,jpj,jpk,jcpoc) )
 
@@ -98 +100 @@
       orem  (:,:,:)   = 0.
       ztremint(:,:,:) = 0.
+      zfolimi(:,:,:)  = 0.
 
       DO jn = 1, jcpoc
@@ -217 +220 @@
                   tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zorem2
                   tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + zofer2
+                  zfolimi(ji,jj,jk)   = zofer2
                END DO
             END DO
@@ -247 +251 @@
                   tra(ji,jj,jk,jpgop) = tra(ji,jj,jk,jpgop) - zopop2 * (1. + solgoc)
                   tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) - zofer2 * (1. + solgoc)
+                  zfolimi(ji,jj,jk)   = zofer2
                END DO
             END DO
@@ -422 +427 @@
                     tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) - zorem
                     tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) - zofer
+                    zfolimi(ji,jj,jk)   = zfolimi(ji,jj,jk) + zofer
                   ENDIF
                END DO
@@ -447 +453 @@
                 tra(ji,jj,jk,jpdop) = tra(ji,jj,jk,jpdop) + zopop 
                 tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + zofer 
+                zfolimi(ji,jj,jk)   = zfolimi(ji,jj,jk) + zofer
              END DO
            END DO
@@ -457 +464 @@
           CALL iom_put( "REMINP" , zremipoc(:,:,:) * tmask(:,:,:) )  ! Remineralisation rate
           CALL iom_put( "REMING" , zremigoc(:,:,:) * tmask(:,:,:) )  ! Remineralisation rate
+          CALL iom_put( "REMINF" , zfolimi(:,:,:) * tmask(:,:,:)  * 1.e+9 * zrfact2 )  ! Remineralisation rate
         ENDIF
      ENDIF
@@ -468 +476 @@
       CALL wrk_dealloc( jpi, jpj,      totprod,  totthick, totcons )
       CALL wrk_dealloc( jpi, jpj, jpk, zremipoc, zremigoc, zorem3, ztremint )
+      CALL wrk_dealloc( jpi, jpj, jpk, zfolimi )
       DEALLOCATE( alphag )
       !

branches/CNRS/dev_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90

@@ -41 +41 @@
    REAL(wp), PUBLIC ::  grosip          !:
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   prmax    !: optimal production = f(temperature)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   prmaxn    !: optimal production = f(temperature)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   quotan   !: proxy of N quota in Nanophyto
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   quotad   !: proxy of N quota in diatomee
@@ -83 +83 @@
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zpronewn, zpronewd
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zmxl_fac, zmxl_chl
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zpligprod1, zpligprod2
       !!---------------------------------------------------------------------
       !
@@ -92 +93 @@
       CALL wrk_alloc( jpi, jpj, jpk, zmxl_fac, zmxl_chl )
       CALL wrk_alloc( jpi, jpj, jpk, zprorcan, zprorcad, zprofed, zprofen, zpronewn, zpronewd )
+     IF( ln_ligand ) THEN
+        CALL wrk_alloc( jpi, jpj, jpk, zpligprod1, zpligprod2 )
+     ENDIF
       !
       zprorcan(:,:,:) = 0._wp ; zprorcad(:,:,:) = 0._wp ; zprofed (:,:,:) = 0._wp
@@ -100 +104 @@
 
       ! Computation of the optimal production
-      prmax(:,:,:) = 0.8_wp * r1_rday * tgfunc(:,:,:) 
+      prmaxn(:,:,:) = 0.8_wp * r1_rday * tgfunc(:,:,:) 
+
+      ! Assign the reference growth rate that is used in implicit 
+      ! biological processes
+      prmax(:,:,:)  = prmaxn(:,:,:)
 
       ! compute the day length depending on latitude and the day
@@ -132 +140 @@
       END DO
 
-      zprbio(:,:,:) = prmax(:,:,:) * zmxl_fac(:,:,:)
+      zprbio(:,:,:) = prmaxn(:,:,:) * zmxl_fac(:,:,:)
       zprdia(:,:,:) = zprbio(:,:,:)
 
@@ -173 +181 @@
                       !  Computation of production function for Chlorophyll
                       !--------------------------------------------------
-                      zpislopen = zpislopeadn(ji,jj,jk) / ( prmax(ji,jj,jk) * zmxl_chl(ji,jj,jk) * rday + rtrn )
-                      zpisloped = zpislopeadd(ji,jj,jk) / ( prmax(ji,jj,jk) * zmxl_chl(ji,jj,jk) * rday + rtrn )
-                      zprnch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislopen * enano(ji,jj,jk) ) )
-                      zprdch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediat(ji,jj,jk) ) )
+                      zpislopen = zpislopeadn(ji,jj,jk) / ( prmaxn(ji,jj,jk) * zmxl_chl(ji,jj,jk) * rday + rtrn )
+                      zpisloped = zpislopeadd(ji,jj,jk) / ( prmaxn(ji,jj,jk) * zmxl_chl(ji,jj,jk) * rday + rtrn )
+                      zprnch(ji,jj,jk) = prmaxn(ji,jj,jk) * ( 1.- EXP( -zpislopen * enano(ji,jj,jk) ) )
+                      zprdch(ji,jj,jk) = prmaxn(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediat(ji,jj,jk) ) )
                   ENDIF
                END DO
@@ -196 +204 @@
                       zpislopen = zpislopen * zmxl_fac(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
                       zpisloped = zpisloped * zmxl_fac(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
-                      zprnch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislopen * enano(ji,jj,jk) ) )
-                      zprdch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediat(ji,jj,jk) ) )
+                      zprnch(ji,jj,jk) = prmaxn(ji,jj,jk) * ( 1.- EXP( -zpislopen * enano(ji,jj,jk) ) )
+                      zprdch(ji,jj,jk) = prmaxn(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediat(ji,jj,jk) ) )
                   ENDIF
                END DO
@@ -210 +218 @@
             DO ji = 1, jpi
                 zval = MIN( xnanopo4(ji,jj,jk), ( xnanonh4(ji,jj,jk) + xnanono3(ji,jj,jk) ) )   &
-                &      * prmax(ji,jj,jk) / ( zprbio(ji,jj,jk) + rtrn )
+                &      * prmaxn(ji,jj,jk) / ( zprbio(ji,jj,jk) + rtrn )
                 quotan(ji,jj,jk) = MIN( 1., 0.2 + 0.8 * zval )
                 zval = MIN( xdiatpo4(ji,jj,jk), ( xdiatnh4(ji,jj,jk) + xdiatno3(ji,jj,jk) ) )   &
-                &      * prmax(ji,jj,jk) / ( zprdia(ji,jj,jk) + rtrn )
+                &      * prmaxn(ji,jj,jk) / ( zprdia(ji,jj,jk) + rtrn )
                 quotad(ji,jj,jk) = MIN( 1., 0.2 + 0.8 * zval )
             END DO
@@ -231 +239 @@
                    !    to mimic the very high ratios observed in the Southern Ocean (silpot2)
                   zlim  = trb(ji,jj,jk,jpsil) / ( trb(ji,jj,jk,jpsil) + xksi1 )
-                  zsilim = MIN( zprdia(ji,jj,jk) / ( prmax(ji,jj,jk) + rtrn ), xlimsi(ji,jj,jk) )
+                  zsilim = MIN( zprdia(ji,jj,jk) / ( prmaxn(ji,jj,jk) + rtrn ), xlimsi(ji,jj,jk) )
                   zsilfac = 4.4 * EXP( -4.23 * zsilim ) * MAX( 0.e0, MIN( 1., 2.2 * ( zlim - 0.5 ) )  ) + 1.e0
                   zsiborn = trb(ji,jj,jk,jpsil) * trb(ji,jj,jk,jpsil) * trb(ji,jj,jk,jpsil)
@@ -253 +261 @@
                zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
                zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
-               zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
-               zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
             END DO
          END DO
@@ -270 +276 @@
                   zratio = trb(ji,jj,jk,jpnfe) / ( trb(ji,jj,jk,jpphy) * fecnm + rtrn )
                   zmax   = MAX( 0., ( 1. - zratio ) / ABS( 1.05 - zratio ) ) 
-                  zprofen(ji,jj,jk) = fecnm * prmax(ji,jj,jk) * ( 1.0 - fr_i(ji,jj) )  &
+                  zprofen(ji,jj,jk) = fecnm * prmaxn(ji,jj,jk) * ( 1.0 - fr_i(ji,jj) )  &
                   &             * ( 4. - 4.5 * xlimnfe(ji,jj,jk) / ( xlimnfe(ji,jj,jk) + 0.5 ) )    &
                   &             * biron(ji,jj,jk) / ( biron(ji,jj,jk) + concnfe(ji,jj,jk) )  &
@@ -280 +286 @@
                   zratio = trb(ji,jj,jk,jpdfe) / ( trb(ji,jj,jk,jpdia) * fecdm + rtrn )
                   zmax   = MAX( 0., ( 1. - zratio ) / ABS( 1.05 - zratio ) ) 
-                  zprofed(ji,jj,jk) = fecdm * prmax(ji,jj,jk) * ( 1.0 - fr_i(ji,jj) )  &
+                  zprofed(ji,jj,jk) = fecdm * prmaxn(ji,jj,jk) * ( 1.0 - fr_i(ji,jj) )  &
                   &             * ( 4. - 4.5 * xlimdfe(ji,jj,jk) / ( xlimdfe(ji,jj,jk) + 0.5 ) )    &
                   &             * biron(ji,jj,jk) / ( biron(ji,jj,jk) + concdfe(ji,jj,jk) )  &
@@ -355 +361 @@
                     zfeup    = texcretn * zprofen(ji,jj,jk) + texcretd * zprofed(ji,jj,jk)
                     tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) + zdocprod * ldocp - zfeup * plig(ji,jj,jk) * lthet
+                    zpligprod1(ji,jj,jk) = zdocprod * ldocp
+                    zpligprod2(ji,jj,jk) = zfeup * plig(ji,jj,jk) * lthet
                  ENDIF
               END DO
@@ -398 +406 @@
           ENDIF
           IF( iom_use( "Mumax" ) )  THEN
-              zw3d(:,:,:) = prmax(:,:,:) * tmask(:,:,:)   ! Maximum growth rate
+              zw3d(:,:,:) = prmaxn(:,:,:) * tmask(:,:,:)   ! Maximum growth rate
               CALL iom_put( "Mumax"  , zw3d )
           ENDIF
@@ -409 +417 @@
           ENDIF
           IF( iom_use( "LNlight" ) .OR. iom_use( "LDlight" ) )  THEN
-              zw3d(:,:,:) = zprbio (:,:,:) / (prmax(:,:,:) + rtrn) * tmask(:,:,:) ! light limitation term
+              zw3d(:,:,:) = zprbio (:,:,:) / (prmaxn(:,:,:) + rtrn) * tmask(:,:,:) ! light limitation term
               CALL iom_put( "LNlight"  , zw3d )
               !
-              zw3d(:,:,:) =  zprdia (:,:,:) / (prmax(:,:,:) + rtrn) * tmask(:,:,:)  ! light limitation term
+              zw3d(:,:,:) =  zprdia (:,:,:) / (prmaxn(:,:,:) + rtrn) * tmask(:,:,:)  ! light limitation term
               CALL iom_put( "LDlight"  , zw3d )
           ENDIF
@@ -426 +434 @@
               zw3d(:,:,:) = ( zprofen(:,:,:) + zprofed(:,:,:) ) * zfact * tmask(:,:,:)  ! total biogenic iron production
               CALL iom_put( "TPBFE"  , zw3d )
+          ENDIF
+          IF( iom_use( "LPRODP" ) )  THEN
+              zw3d(:,:,:) = zpligprod1(:,:,:) * 1e9 * zfact * tmask(:,:,:)
+              CALL iom_put( "LPRODP"  , zw3d )
+          ENDIF
+          IF( iom_use( "LDETP" ) )  THEN
+              zw3d(:,:,:) = zpligprod2(:,:,:) * 1e9 * zfact * tmask(:,:,:)
+              CALL iom_put( "LDETP"  , zw3d )
           ENDIF
           IF( iom_use( "INTPPPHYN" ) .OR. iom_use( "INTPPPHYD" ) ) THEN  
@@ -485 +501 @@
      CALL wrk_dealloc( jpi, jpj, jpk, zmxl_fac, zmxl_chl )
      CALL wrk_dealloc( jpi, jpj, jpk, zprorcan, zprorcad, zprofed, zprofen, zpronewn, zpronewd )
+     IF( ln_ligand ) THEN
+        CALL wrk_dealloc( jpi, jpj, jpk, zpligprod1, zpligprod2 )
+     ENDIF
      !
      IF( nn_timing == 1 )  CALL timing_stop('p4z_prod')
@@ -550 +569 @@
       !!                     ***  ROUTINE p4z_prod_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( prmax(jpi,jpj,jpk), quotan(jpi,jpj,jpk), quotad(jpi,jpj,jpk), STAT = p4z_prod_alloc )
+      ALLOCATE( prmaxn(jpi,jpj,jpk), quotan(jpi,jpj,jpk), quotad(jpi,jpj,jpk), STAT = p4z_prod_alloc )
       !
       IF( p4z_prod_alloc /= 0 ) CALL ctl_warn('p4z_prod_alloc : failed to allocate arrays.')

branches/CNRS/dev_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90

@@ -16 +16 @@
    USE sms_pisces      !  PISCES Source Minus Sink variables
    USE p4zche          !  chemical model
-   USE p4zprod         !  Growth rate of the 2 phyto groups
    USE p4zlim
    USE prtctl_trc      !  print control for debugging
@@ -70 +69 @@
       REAL(wp), POINTER, DIMENSION(:,:  ) :: ztempbac
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zdepbac, zolimi, zdepprod, zfacsi, zw3d, zfacsib
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zfebact, zdepeff
       !!---------------------------------------------------------------------
       !
@@ -76 +76 @@
       ! Allocate temporary workspace
       CALL wrk_alloc( jpi, jpj,      ztempbac                  )
-      CALL wrk_alloc( jpi, jpj, jpk, zdepbac, zdepprod, zolimi, zfacsi, zfacsib )
+      CALL wrk_alloc( jpi, jpj, jpk, zdepbac, zdepprod, zolimi, zfacsi, zfacsib, zfebact, zdepeff )
 
       ! Initialisation of temprary arrys
       zdepprod(:,:,:) = 1._wp
+      zdepeff (:,:,:) = 0.3_wp
       ztempbac(:,:)   = 0._wp
+      zfebact (:,:,:) = 0._wp
       zfacsib(:,:,:)  = xsilab / ( 1.0 - xsilab )
       zfacsi(:,:,:)   = xsilab
@@ -100 +102 @@
                   zdepbac (ji,jj,jk) = zdepmin**0.683 * ztempbac(ji,jj)
                   zdepprod(ji,jj,jk) = zdepmin**0.273
+                  zdepeff (ji,jj,jk) = 0.3 * zdepmin**0.3
                ENDIF
             END DO
@@ -220 +223 @@
                zbactfer = feratb *  rfact2 * prmax(ji,jj,jk) * xlimbacl(ji,jj,jk)             &
                   &              * trb(ji,jj,jk,jpfer) / ( xkferb + trb(ji,jj,jk,jpfer) )    &
-                  &              * zdepprod(ji,jj,jk) * zdepbac(ji,jj,jk)
-               tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zbactfer*0.16
-               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zbactfer*0.12
-               tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zbactfer*0.04
+                  &              * zdepprod(ji,jj,jk) * zdepeff(ji,jj,jk) * zdepbac(ji,jj,jk)
+               tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zbactfer*0.33
+               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zbactfer*0.25
+               tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zbactfer*0.08
+               zfebact(ji,jj,jk)   = zbactfer * 0.33
             END DO
          END DO
@@ -283 +287 @@
           ENDIF
           !
+          IF( iom_use( "BACT" ) )  THEN
+               zw3d(:,:,:) = zdepbac(:,:,:) * 1.E6 * tmask(:,:,:)  ! Bacterial biomass
+               CALL iom_put( "BACT", zw3d )
+          ENDIF
+          IF( iom_use( "FEBACT" ) )  THEN
+               zw3d(:,:,:) = zfebact(:,:,:) * 1E9 * tmask(:,:,:) * zrfact2   ! Bacterial iron consumption
+               CALL iom_put( "FEBACT" , zw3d )  
+          ENDIF
           CALL wrk_dealloc( jpi, jpj, jpk, zw3d )
        ENDIF
       !
       CALL wrk_dealloc( jpi, jpj,      ztempbac                  )
-      CALL wrk_dealloc( jpi, jpj, jpk, zdepbac, zdepprod, zolimi, zfacsi, zfacsib )
+      CALL wrk_dealloc( jpi, jpj, jpk, zdepbac, zdepprod, zolimi, zfacsi, zfacsib, zfebact, zdepeff )
       !
       IF( nn_timing == 1 )  CALL timing_stop('p4z_rem')

branches/CNRS/dev_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90

@@ -170 +170 @@
                      rivdop(ji,jj) = ( sf_river(jr_dop)%fnow(ji,jj,1) ) &
                         &              * 1.E3 / po4r / ( 31. * zcoef + rtrn ) * tmask(ji,jj,1)
+                     rivdsi(ji,jj) =   sf_river(jr_dsi)%fnow(ji,jj,1)                                    &
+                        &              * 1.E3        / ( 28.1 * zcoef + rtrn )
                   END DO
                END DO

branches/CNRS/dev_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p5zmeso.F90

@@ -85 +85 @@
       CHARACTER (len=25) :: charout
       REAL(wp) :: zrfact2, zmetexcess
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zgrazing, zw3d
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zgrazing, zw3d, zfezoo2
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zz2ligprod
 
       !!---------------------------------------------------------------------
@@ -91 +92 @@
       IF( nn_timing == 1 )  CALL timing_start('p5z_meso')
       !
-      CALL wrk_alloc( jpi, jpj, jpk, zgrazing )
+      CALL wrk_alloc( jpi, jpj, jpk, zgrazing, zfezoo2 )
       zgrazing(:,:,:) = 0._wp
+      zfezoo2 (:,:,:) = 0._wp
+      !
+      IF (ln_ligand) THEN
+         CALL wrk_alloc( jpi, jpj, jpk, zz2ligprod )
+         zz2ligprod(:,:,:) = 0._wp
+      ENDIF
 
       zmetexcess = 0.0
@@ -138 +145 @@
                !   most abundant species
                !   ------------------------------------------------------------  
-               ztmp1 = xpref2p * zcompaph**1.5
-               ztmp2 = xpref2m * zcompames**1.5
-               ztmp3 = xpref2c * zcompapoc**1.5
-               ztmp4 = xpref2d * zcompadi**1.5
-               ztmp5 = xpref2z * zcompaz**1.5
+               ztmp1 = xpref2p * zcompaph * (0.75 * zcompaph + 0.25 * zcompadi)
+               ztmp2 = xpref2m * zcompames**2
+               ztmp3 = xpref2c * zcompapoc**2
+               ztmp4 = xpref2d * zcompadi * (0.75 * zcompadi + 0.25 * zcompaph)
+               ztmp5 = xpref2z * zcompaz**2
                ztmptot = ztmp1 + ztmp2 + ztmp3 + ztmp4 + ztmp5 + rtrn
                ztmp1 = ztmp1 / ztmptot
@@ -290 +297 @@
                tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zgradoc
                !
-               IF( ln_ligand ) tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) + zgradoc * ldocz
+               IF( ln_ligand ) THEN
+                  tra(ji,jj,jk,jplgw)  = tra(ji,jj,jk,jplgw) + zgradoc * ldocz
+                  zz2ligprod(ji,jj,jk) = zgradoc * ldocz
+               ENDIF
                !
                tra(ji,jj,jk,jpdon) = tra(ji,jj,jk,jpdon) + zgradon
@@ -296 +306 @@
                tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) - o2ut * zgrarem
                tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + zgraref
+               zfezoo2(ji,jj,jk)   = zgraref
                tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) + zgrarem
                tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * zgraren
@@ -351 +362 @@
             CALL iom_put( "PCAL", zw3d )
          ENDIF
+         IF( iom_use( "FEZOO2" ) ) THEN
+            zw3d(:,:,:) = zfezoo2(:,:,:) * 1e9 * 1.e+3 * rfact2r * tmask(:,:,:)   !
+            CALL iom_put( "FEZOO2", zw3d )
+         ENDIF
+         IF( iom_use( "LPRODZ2" ) .AND. ln_ligand )  THEN
+            zw3d(:,:,:) = zz2ligprod(:,:,:) * 1e9 * 1.e+3 * rfact2r * tmask(:,:,:)
+            CALL iom_put( "LPRODZ2"  , zw3d )
+         ENDIF
          CALL wrk_dealloc( jpi, jpj, jpk, zw3d )
       ENDIF
@@ -360 +379 @@
       ENDIF
       !
-      CALL wrk_dealloc( jpi, jpj, jpk, zgrazing )
+      CALL wrk_dealloc( jpi, jpj, jpk, zgrazing, zfezoo2 )
+      IF ( ln_ligand ) CALL wrk_dealloc( jpi, jpj, jpk, zz2ligprod )
       !
       IF( nn_timing == 1 )  CALL timing_stop('p5z_meso')

branches/CNRS/dev_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p5zmicro.F90

@@ -15 +15 @@
    USE trc             !  passive tracers common variables 
    USE sms_pisces      !  PISCES Source Minus Sink variables
-   USE p5zlim          !  Phytoplankton limitation terms
+   USE p4zlim          !  Phytoplankton limitation terms
    USE iom             !  I/O manager
    USE prtctl_trc      !  print control for debugging
@@ -82 +82 @@
       REAL(wp) :: zgrazdc, zgrazdn, zgrazdp, zgrazdf, zgraznf, zgrazz
       REAL(wp) :: zgrazpc, zgrazpn, zgrazpp, zgrazpf, zbeta, zrfact2, zmetexcess
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zgrazing, zw3d
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zgrazing, zw3d, zfezoo
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zzligprod
       CHARACTER (len=25) :: charout
       !!---------------------------------------------------------------------
@@ -88 +89 @@
       IF( nn_timing == 1 )  CALL timing_start('p5z_micro')
       !
-      CALL wrk_alloc( jpi, jpj, jpk, zgrazing )
+      CALL wrk_alloc( jpi, jpj, jpk, zgrazing, zfezoo )
+      !
+      IF (ln_ligand) THEN
+         CALL wrk_alloc( jpi, jpj, jpk, zzligprod )
+         zzligprod(:,:,:) = 0._wp
+      ENDIF
       !
       zmetexcess = 0.0
@@ -134 +140 @@
                !   most abundant species
                !   ------------------------------------------------------------  
-               ztmp1 = xprefn * zcompaph**1.5
-               ztmp2 = xprefp * zcompapi**1.5
-               ztmp3 = xprefc * zcompapoc**1.5
-               ztmp4 = xprefd * zcompadi**1.5
-               ztmp5 = xprefz * zcompaz**1.5
+               ztmp1 = xprefn * zcompaph * ( 0.75 * zcompaph + 0.25 * zcompadi )
+               ztmp2 = xprefp * zcompapi * zcompapi
+               ztmp3 = xprefc * zcompapoc * zcompapoc
+               ztmp4 = xprefd * zcompadi * ( 0.75 * zcompadi + 0.25 * zcompaph )
+               ztmp5 = xprefz * zcompaz * zcompaz
                ztmptot = ztmp1 + ztmp2 + ztmp3 + ztmp4 + ztmp5 + rtrn
                ztmp1 = ztmp1 / ztmptot
@@ -245 +251 @@
                tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zgradoc
                !
-               IF( ln_ligand ) tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) + zgradoc * ldocz
+               IF( ln_ligand ) THEN
+                  tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) + zgradoc * ldocz
+                  zzligprod(ji,jj,jk) = zgradoc * ldocz
+               ENDIF
                !
                tra(ji,jj,jk,jpdon) = tra(ji,jj,jk,jpdon) + zgradon
@@ -251 +260 @@
                tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) - o2ut * zgrarem 
                tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + zgraref
+               zfezoo(ji,jj,jk) = zgraref
                tra(ji,jj,jk,jpzoo) = tra(ji,jj,jk,jpzoo) + zepsherv * zgraztotc - zrespirc - ztortz - zgrazz
                tra(ji,jj,jk,jpphy) = tra(ji,jj,jk,jpphy) - zgraznc
@@ -295 +305 @@
             CALL iom_put( "GRAZ1", zw3d )
          ENDIF
+         IF( iom_use( "FEZOO" ) ) THEN
+            zw3d(:,:,:) = zfezoo(:,:,:) * 1e9 * 1.e+3 * rfact2r * tmask(:,:,:)   !
+            CALL iom_put( "FEZOO", zw3d )
+         ENDIF
+         IF( iom_use( "LPRODZ" ) .AND. ln_ligand )  THEN
+            zw3d(:,:,:) = zzligprod(:,:,:) * 1e9 * 1.e+3 * rfact2r * tmask(:,:,:)
+            CALL iom_put( "LPRODZ"  , zw3d )
+         ENDIF
          CALL wrk_dealloc( jpi, jpj, jpk, zw3d )
       ENDIF
@@ -304 +322 @@
       ENDIF
       !
-      CALL wrk_dealloc( jpi, jpj, jpk, zgrazing )
+      CALL wrk_dealloc( jpi, jpj, jpk, zgrazing, zfezoo )
+      IF (ln_ligand) CALL wrk_dealloc( jpi, jpj, jpk, zzligprod )
       !
       IF( nn_timing == 1 )  CALL timing_stop('p5z_micro')

branches/CNRS/dev_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p5zmort.F90

@@ -14 +14 @@
    USE trc             !  passive tracers common variables 
    USE sms_pisces      !  PISCES Source Minus Sink variables
-   USE p5zlim          !  Phytoplankton limitation terms
+   USE p4zlim          !  Phytoplankton limitation terms
    USE prtctl_trc      !  print control for debugging
 

branches/CNRS/dev_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p5zprod.F90

@@ -16 +16 @@
    USE trc             !  passive tracers common variables 
    USE sms_pisces      !  PISCES Source Minus Sink variables
-   USE p5zlim          !  Co-limitations of differents nutrients
+   USE p4zlim          !  Co-limitations of differents nutrients
    USE prtctl_trc      !  print control for debugging
    USE iom             !  I/O manager
@@ -78 +78 @@
       REAL(wp) ::   zprnutmax, zdocprod, zprochln, zprochld, zprochlp
       REAL(wp) ::   zpislopen, zpislopep, zpisloped, thetannm_n, thetandm_n, thetanpm_n
+      REAL(wp) ::   zratiosi, zquotasi, grosip2
       REAL(wp) ::   zrum, zcodel, zargu, zval, zfeup
       REAL(wp) ::   zfact, zrfact2
@@ -94 +95 @@
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zcroissn, zcroissp, zcroissd
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zmxl_fac, zmxl_chl
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zpligprod1, zpligprod2
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d
       REAL(wp), POINTER, DIMENSION(:,:  ) :: zw2d
@@ -110 +112 @@
       CALL wrk_alloc( jpi, jpj, jpk, zprchln, zprchlp, zprchld, zprodopn, zprodopp, zprodopd )
       CALL wrk_alloc( jpi, jpj, jpk, zcroissp, zcroissn, zcroissd )
+      IF( ln_ligand ) THEN
+         CALL wrk_alloc( jpi, jpj, jpk, zpligprod1, zpligprod2 )
+      ENDIF
       !
       zprorcan(:,:,:) = 0._wp ; zprorcap(:,:,:) = 0._wp ; zprorcad(:,:,:) = 0._wp
@@ -121 +126 @@
 
       ! Computation of the optimal production
-      prmaxn(:,:,:) = ( 0.65_wp * (1. + zpsino3 * qnpmax ) ) * r1_rday * tgfunc(:,:,:)
-      prmaxp(:,:,:) = 0.5 / 0.65 * prmaxn(:,:,:) 
+      prmaxn(:,:,:) = ( 0.8_wp * (1. + zpsino3 * qnpmax ) ) * r1_rday * tgfunc(:,:,:)
+      prmaxp(:,:,:) = 0.5 / 0.8 * prmaxn(:,:,:) 
       prmaxd(:,:,:) = prmaxn(:,:,:) 
-      zprnut(:,:,:) = 0.65_wp * r1_rday * tgfunc(:,:,:)
+      zprnut(:,:,:) = 0.8_wp * r1_rday * tgfunc(:,:,:)
+
+      ! Assign the reference growth rate used in implicit formulations
+      prmax(:,:,:)  = zprnut(:,:,:)
 
       ! compute the day length depending on latitude and the day
@@ -213 +221 @@
                   !    Si/C is arbitrariliy increased for very high Si concentrations
                   !    to mimic the very high ratios observed in the Southern Ocean (silpot2)
+                  IF (gphit(ji,jj) < -30 ) THEN
+                    zsilfac2 = 1. + zsiborn / ( zsiborn + xksi2**3 )
+                  ELSE
+                    zsilfac2 = 1. 
+                  ENDIF
+                  zsilim = MIN( zprdia(ji,jj,jk) / ( prmaxd(ji,jj,jk) + rtrn ), xlimsi2(ji,jj,jk) )
+                  zysopt(ji,jj,jk) = grosip * zlim * zsilfac * zsilfac2
                   zlim  = trb(ji,jj,jk,jpsil) / ( trb(ji,jj,jk,jpsil) + xksi1 )
-                  zsilim = MIN( zprdia(ji,jj,jk) / ( prmaxd(ji,jj,jk) + rtrn ), xlimsi(ji,jj,jk) )
-                  zsilfac = 3.4 * EXP( -4.23 * zsilim ) * MAX( 0.e0, MIN( 1., 2.2 * ( zlim - 0.5 ) )  ) + 1.e0
-                  zsiborn = trb(ji,jj,jk,jpsil) * trb(ji,jj,jk,jpsil) * trb(ji,jj,jk,jpsil)
-                  IF (gphit(ji,jj) < -30 ) THEN
-                    zsilfac2 = 1. + 2. * zsiborn / ( zsiborn + xksi2**3 )
-                  ELSE
-                    zsilfac2 = 1. +      zsiborn / ( zsiborn + xksi2**3 )
-                  ENDIF
-                  zysopt(ji,jj,jk) = grosip * zlim * zsilfac * zsilfac2
+                  zratiosi = trb(ji,jj,jk,jpndi)
+                  grosip2 = zsilfac2 * grosip
+                  zquotasi = MIN( grosip2, trn(ji,jj,jk,jpdsi) / (zratiosi + rtrn ) )
+                  zysopt(ji,jj,jk) = 1./3. * zratiosi * zlim * MAX(0.0, grosip2 - zquotasi)
+                  zysopt(ji,jj,jk) = zysopt(ji,jj,jk) * 1.005 * zsilim**3 / (0.005 + zsilim**3 )
               ENDIF
             END DO
@@ -332 +343 @@
                   zmax = MAX(0., MIN(1., (1. - zrat)/ (1.05 - zrat) * 1.05)) 
                   zpronmax = zprnutmax * zmax * MAX(0., MIN(1., ( zratiop - xqpdmin(ji,jj,jk) )   &
-                  &          / ( xqpdmax(ji,jj,jk) - xqpdmin(ji,jj,jk) + rtrn ), xlimdfe(ji,jj,jk) ) )
+                  &          / ( xqpdmax(ji,jj,jk) - xqpdmin(ji,jj,jk) + rtrn ), xlimsi(ji,jj,jk) ) )
                   zpronewd(ji,jj,jk) = zpronmax * zdaylen(ji,jj) * xdiatno3(ji,jj,jk)
                   zproregd(ji,jj,jk) = zpronmax * xdiatnh4(ji,jj,jk)
@@ -338 +349 @@
                   zrat = MIN( 1., zratiop / (xqpdmax(ji,jj,jk) + rtrn) )
                   zmax = MAX(0., MIN(1., (1. - zrat)/ (1.05 - zrat) * 1.05)) 
-                  zpropmax = zprnutmax * zmax * xlimdfe(ji,jj,jk)
+                  zpropmax = zprnutmax * zmax * xlimsi(ji,jj,jk)
                   zpropo4d(ji,jj,jk) = zpropmax * xdiatpo4(ji,jj,jk)
                   zprodopd(ji,jj,jk) = zpropmax * xdiatdop(ji,jj,jk)
@@ -362 +373 @@
                   thetannm_n   = MIN ( thetannm, ( thetannm / (1. - 1.14 / 43.4 *tsn(ji,jj,jk,jp_tem)))   &
                   &               * (1. - 1.14 / 43.4 * 20.))
+                  thetannm_n = thetannm
                   zprochln = thetannm_n * zprod / ( zpislopeadn(ji,jj,jk) * znanotot + rtrn )
                   zprochln = MAX(zprochln, chlcmin * 12. * zprorcan (ji,jj,jk) )
@@ -369 +381 @@
                   thetanpm_n   = MIN ( thetanpm, ( thetanpm / (1. - 1.14 / 43.4 *tsn(ji,jj,jk,jp_tem)))   &
                   &               * (1. - 1.14 / 43.4 * 20.))
+                  thetanpm_n = thetanpm
                   zprochlp = thetanpm_n * zprod / ( zpislopeadp(ji,jj,jk) * zpicotot + rtrn )
                   zprochlp = MAX(zprochlp, chlcmin * 12. * zprorcap(ji,jj,jk) )
@@ -376 +389 @@
                   thetandm_n   = MIN ( thetandm, ( thetandm / (1. - 1.14 / 43.4 *tsn(ji,jj,jk,jp_tem)))   &
                   &               * (1. - 1.14 / 43.4 * 20.))
+                  thetandm_n = thetandm
                   zprochld = thetandm_n * zprod / ( zpislopeadd(ji,jj,jk) * zdiattot + rtrn )
                   zprochld = MAX(zprochld, chlcmin * 12. * zprorcad(ji,jj,jk) )
@@ -460 +474 @@
                  zfeup    = texcretn * zprofen(ji,jj,jk) + texcretd * zprofed(ji,jj,jk) + texcretp * zprofep(ji,jj,jk)
                  tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) + zdocprod * ldocp - zfeup * plig(ji,jj,jk) * lthet
+                 zpligprod1(ji,jj,jk) = zdocprod * ldocp
+                 zpligprod2(ji,jj,jk) = zfeup * plig(ji,jj,jk) * lthet
               END DO
            END DO
@@ -547 +563 @@
               !
           ENDIF
-
+          IF( iom_use( "LPRODP" ) )  THEN
+              zw3d(:,:,:) = zpligprod1(:,:,:) * 1e9 * zfact * tmask(:,:,:)
+              CALL iom_put( "LPRODP"  , zw3d )
+          ENDIF
+          IF( iom_use( "LDETP" ) )  THEN
+              zw3d(:,:,:) = zpligprod2(:,:,:) * 1e9 * zfact * tmask(:,:,:)
+              CALL iom_put( "LDETP"  , zw3d )
+          ENDIF
+          !
           IF( iom_use( "tintpp" ) )  CALL iom_put( "tintpp" , tpp * zfact )  !  global total integrated primary production molC/s
           !
@@ -570 +594 @@
       CALL wrk_dealloc( jpi, jpj, jpk, zprchln, zprchlp, zprchld, zprodopn, zprodopp, zprodopd )
       CALL wrk_dealloc( jpi, jpj, jpk, zcroissp, zcroissn, zcroissd )
+      IF( ln_ligand ) THEN
+         CALL wrk_dealloc( jpi, jpj, jpk, zpligprod1, zpligprod2 )
+      ENDIF
       !
       IF( nn_timing == 1 )  CALL timing_stop('p5z_prod')

branches/CNRS/dev_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES/sms_pisces.F90

@@ -100 +100 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   orem       !: ??
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xdiss      !: ??
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   prmax      !: Reference maximum growth rate
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   prodcal    !: Calcite production
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   prodpoc    !: Calcite production
@@ -164 +165 @@
             &      prodcal(jpi,jpj,jpk) , xdiss   (jpi,jpj,jpk),    &
             &      prodpoc(jpi,jpj,jpk) , conspoc(jpi,jpj,jpk) ,    &
-            &      prodgoc(jpi,jpj,jpk) , consgoc(jpi,jpj,jpk) ,  STAT=ierr(4) )
+            &      prodgoc(jpi,jpj,jpk) , consgoc(jpi,jpj,jpk) ,    &
+            &      prmax  (jpi,jpj,jpk) , STAT=ierr(4) )
 
          !* Variable for chemistry of the CO2 cycle

branches/CNRS/dev_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES/trcini_pisces.F90

@@ -71 +71 @@
       USE p4zpoc          !  Remineralization of organic particles
       USE p4zligand       !  Remineralization of organic ligands
-      USE p5zlim          !  Co-limitations of differents nutrients
+      USE p4zlim          !  Co-limitations of differents nutrients
       USE p5zprod         !  Growth rate of the 2 phyto groups
       USE p5zmicro        !  Sources and sinks of microzooplankton
@@ -112 +112 @@
       ierr = ierr +  p4z_sed_alloc()
       ierr = ierr +  p4z_rem_alloc()
+      ierr = ierr +  p4z_lim_alloc()
       IF( ln_p4z ) THEN
-         ierr = ierr +  p4z_lim_alloc()
          ierr = ierr +  p4z_prod_alloc()
       ELSE
-         ierr = ierr +  p5z_lim_alloc()
          ierr = ierr +  p5z_prod_alloc()
       ENDIF
@@ -167 +166 @@
         IF( cltra == 'PICN'     )   jpnpi = jn      !: Picophytoplankton N biomass
         IF( cltra == 'PICP'     )   jpppi = jn      !: Picophytoplankton P biomass
+        IF( cltra == 'PCHL'     )   jppch = jn      !: Diatoms Chlorophyll Concentration
         IF( cltra == 'PFe'      )   jppfe = jn      !: Picophytoplankton Fe biomass
         IF( cltra == 'LGW'      )   jplgw = jn      !: Weak ligands
@@ -240 +240 @@
          xksi(:,:)    = 2.e-6
          xksimax(:,:) = xksi(:,:)
+         !
+         IF( ln_p5z ) THEN
+            sizep(:,:,:) = 1.
+            sizen(:,:,:) = 1.
+            sized(:,:,:) = 1.
+         ENDIF
       END IF