Changeset 14276 for NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z

Timestamp:

2021-01-07T23:09:56+01:00 (3 years ago)

Author:

aumont

Message:

numerous updates to PISCES, PISCES-QUOTA and the sediment module

Location:

NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z

Files:

• p4zagg.F90 (modified) (5 diffs)
• p4zbio.F90 (modified) (1 diff)
• p4zfechem.F90 (modified) (12 diffs)
• p4zflx.F90 (modified) (1 diff)
• p4zint.F90 (modified) (3 diffs)
• p4zlim.F90 (modified) (10 diffs)
• p4zlys.F90 (modified) (3 diffs)
• p4zmeso.F90 (modified) (8 diffs)
• p4zmicro.F90 (modified) (3 diffs)
• p4zmort.F90 (modified) (7 diffs)
• p4zopt.F90 (modified) (5 diffs)
• p4zprod.F90 (modified) (13 diffs)
• p4zrem.F90 (modified) (12 diffs)
• p4zsbc.F90 (modified) (1 diff)
• p4zsed.F90 (modified) (18 diffs)
• p4zsms.F90 (modified) (13 diffs)
• p5zlim.F90 (modified) (23 diffs)
• p5zmeso.F90 (modified) (14 diffs)
• p5zmicro.F90 (modified) (8 diffs)
• p5zmort.F90 (modified) (8 diffs)
• p5zprod.F90 (modified) (14 diffs)

NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zagg.F90

@@ -75 +75 @@
                   ! Part I : Coagulation dependent on turbulence
                   ! The stickiness has been assumed to be 0.1
-                  zagg1 = 25.9  * zfact * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jppoc)
-                  zagg2 = 4452. * zfact * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jpgoc)
+                  zagg1 = 12.5  * zfact * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jppoc)
+                  zagg2 = 169.7 * zfact * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jpgoc)
 
                   ! Part II : Differential settling
                   ! Aggregation of small into large particles
                   ! The stickiness has been assumed to be 0.1
-                  zagg3 =  47.1 * xstep * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jpgoc)
-                  zagg4 =  3.3  * xstep * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jppoc)
+                  zagg3 =  8.63  * xstep * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jppoc)
+                  zagg4 =  132.8 * xstep * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jpgoc)
 
                   zagg   = zagg1 + zagg2 + zagg3 + zagg4
@@ -92 +92 @@
                   ! 3rd term is differential settling of DOC-POC
                   ! 1/3 of DOC is supposed to experience aggregation (HMW)
-                  zaggdoc  = ( ( 0.369 * 0.3 * trb(ji,jj,jk,jpdoc) + 102.4 * trb(ji,jj,jk,jppoc) ) * zfact       &
-                  &            + 2.4 * xstep * trb(ji,jj,jk,jppoc) ) * 0.3 * trb(ji,jj,jk,jpdoc)
+                  zaggdoc  = ( ( 12.0 * 0.3 * trb(ji,jj,jk,jpdoc) + 9.05 * trb(ji,jj,jk,jppoc) ) * zfact       &
+                  &            + 2.49 * xstep * trb(ji,jj,jk,jppoc) ) * 0.3 * trb(ji,jj,jk,jpdoc)
                   ! transfer of DOC to GOC : 
                   ! 1st term is shear aggregation
-                  ! 2nd term is differential settling 
-                  ! 1/3 of DOC is supposed to experience aggregation (HMW)
-                  zaggdoc2 = ( 3.53E3 * zfact + 0.1 * xstep ) * trb(ji,jj,jk,jpgoc) * 0.3 * trb(ji,jj,jk,jpdoc)
+                  ! 1/3 of DOC is supposed to experience aggregation (HMW)
+                  zaggdoc2 = ( 1.94 * zfact + 1.37 * xstep ) * trb(ji,jj,jk,jpgoc) * 0.3 * trb(ji,jj,jk,jpdoc)
                   ! tranfer of DOC to POC due to brownian motion
                   ! The temperature dependency has been omitted.
-                  zaggdoc3 =  114. * 0.3 * trb(ji,jj,jk,jpdoc) *xstep * 0.3 * trb(ji,jj,jk,jpdoc)
+                  zaggdoc3 =  ( 127.8 * 0.3 * trb(ji,jj,jk,jpdoc) + 725.7 * trb(ji,jj,jk,jppoc) ) * xstep * 0.3 * trb(ji,jj,jk,jpdoc)
 
                   !  Update the trends
@@ -150 +149 @@
                   ! 3rd term is differential settling of DOC-POC
                   ! 1/3 of DOC is supposed to experience aggregation (HMW)
-                  zaggtmp = ( ( 0.369 * 0.3 * trb(ji,jj,jk,jpdoc) + 102.4 * trb(ji,jj,jk,jppoc) ) * zfact       &
-                  &            + 2.4 * xstep * trb(ji,jj,jk,jppoc) )
+                  zaggtmp = ( ( 0.37 * 0.3 * trb(ji,jj,jk,jpdoc) + 20.5 * trb(ji,jj,jk,jppoc) ) * zfact       &
+                  &            + 0.15 * xstep * trb(ji,jj,jk,jppoc) )
                   zaggdoc  = zaggtmp * 0.3 * trb(ji,jj,jk,jpdoc)
                   zaggdon  = zaggtmp * 0.3 * trb(ji,jj,jk,jpdon)
@@ -160 +159 @@
                   ! 2nd term is differential settling 
                   ! 1/3 of DOC is supposed to experience aggregation (HMW)
-                  zaggtmp = ( 3.53E3 * zfact + 0.1 * xstep ) * trb(ji,jj,jk,jpgoc)
+                  zaggtmp = 655.4 * zfact * trb(ji,jj,jk,jpgoc)
                   zaggdoc2 = zaggtmp * 0.3 * trb(ji,jj,jk,jpdoc)
                   zaggdon2 = zaggtmp * 0.3 * trb(ji,jj,jk,jpdon)
@@ -167 +166 @@
                   ! tranfer of DOC to POC due to brownian motion
                   ! 1/3 of DOC is supposed to experience aggregation (HMW)
-                  zaggtmp = ( 114. * 0.3 * trb(ji,jj,jk,jpdoc) ) * xstep
+                  zaggtmp = ( 260.2 * 0.3 * trb(ji,jj,jk,jpdoc) +  418.5 * trb(ji,jj,jk,jppoc) ) * xstep
                   zaggdoc3 =  zaggtmp * 0.3 * trb(ji,jj,jk,jpdoc)
                   zaggdon3 =  zaggtmp * 0.3 * trb(ji,jj,jk,jpdon)

NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zbio.F90

@@ -110 +110 @@
 
       ! Update of the size of the different phytoplankton groups
-      sized(:,:,:) = sizeda(:,:,:)
-      sizen(:,:,:) = sizena(:,:,:)
+      sized(:,:,:) = MAX(1.0, sizeda(:,:,:) )
+      sizen(:,:,:) = MAX(1.0, sizena(:,:,:) )
       IF (ln_p5z) THEN
-         sizep(:,:,:) = sizepa(:,:,:)
+         sizep(:,:,:) = MAX(1.0, sizepa(:,:,:) )
       ENDIF
       !                                                             !

NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zfechem.F90

@@ -30 +30 @@
    REAL(wp), PUBLIC ::   ligand       !: ligand concentration in the ocean 
    REAL(wp), PUBLIC ::   kfep         !: rate constant for nanoparticle formation
+   REAL(wp), PUBLIC ::   scaveff      !: Fraction of scavenged iron that is considered as being subject to solubilization
 
    !!----------------------------------------------------------------------
@@ -50 +51 @@
       !
       INTEGER  ::   ji, jj, jk, jic, jn
-      REAL(wp) ::   zdep, zlam1a, zlam1b, zlamfac
-      REAL(wp) ::   zkeq, zfeequi, zfesatur, zfecoll, fe3sol, zligco
-      REAL(wp) ::   zdenom1, zscave, zaggdfea, zaggdfeb, zcoag, ztrc, zdust
-      REAL(wp) ::   zdenom2, ztfe, zhplus, zxlam, zaggliga, zaggligb
+      REAL(wp) ::   zlam1a, zlam1b
+      REAL(wp) ::   zkeq, zfesatur, zfecoll, fe3sol, zligco
+      REAL(wp) ::   zscave, zaggdfea, zaggdfeb, ztrc, zdust, zklight
+      REAL(wp) ::   ztfe, zhplus, zxlam, zaggliga, zaggligb
       REAL(wp) ::   zrfact2
       CHARACTER (len=25) :: charout
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zTL1, zFe3, ztotlig, precip, zFeL1
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zcoll3d, zscav3d, zlcoll3d
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zTL1, zFe3, ztotlig, precip, precipno3, zFeL1
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zcoll3d, zscav3d, zlcoll3d, zprecip3d
       !!---------------------------------------------------------------------
       !
@@ -66 +67 @@
 
       ! Total ligand concentration : Ligands can be chosen to be constant or variable
-      ! Parameterization from Tagliabue and Voelker (2011)
+      ! Parameterization from Pham and Ito (2018)
       ! -------------------------------------------------
       IF( ln_ligvar ) THEN
-         ztotlig(:,:,:) =  0.09 * trb(:,:,:,jpdoc) * 1E6 + ligand * 1E9
+         ztotlig(:,:,:) =  0.09 * 0.667 * trb(:,:,:,jpdoc) * 1E6 + ligand * 1E9 + MAX(0., chemo2(:,:,:) - trb(:,:,:,jpoxy) ) / 400.E-6
          ztotlig(:,:,:) =  MIN( ztotlig(:,:,:), 10. )
       ELSE
         IF( ln_ligand ) THEN  ;   ztotlig(:,:,:) = trb(:,:,:,jplgw) * 1E9
-        ELSE                  ;   ztotlig(:,:,:) = ligand * 1E9
+        ELSE                  ;   ztotlig(:,:,:) = ligand * 1E9 
         ENDIF
       ENDIF
@@ -79 +80 @@
       ! ------------------------------------------------------------
       !  from Aumont and Bopp (2006)
-      ! This model is based on one ligand and Fe' 
+      ! This model is based on one ligand, Fe2+ and Fe3+ 
       ! Chemistry is supposed to be fast enough to be at equilibrium
       ! ------------------------------------------------------------
@@ -87 +88 @@
                zTL1(ji,jj,jk)  = ztotlig(ji,jj,jk)
                zkeq            = fekeq(ji,jj,jk)
+               zklight         = 4.77E-7 * etot(ji,jj,jk) * 0.5 / 10**-6.3
                zfesatur        = zTL1(ji,jj,jk) * 1E-9
-               ztfe            = trb(ji,jj,jk,jpfer) 
+               ztfe            = (1.0 + zklight) * trb(ji,jj,jk,jpfer) 
                ! Fe' is the root of a 2nd order polynom
-               zFe3 (ji,jj,jk) = ( -( 1. + zfesatur * zkeq - zkeq * ztfe )               &
-                  &              + SQRT( ( 1. + zfesatur * zkeq - zkeq * ztfe )**2       &
+               zFe3 (ji,jj,jk) = ( -( 1. + zfesatur * zkeq + zklight + consfe3(ji,jj,jk)/10**-6.3 - zkeq * trb(ji,jj,jk,jpfer) )               &
+                  &              + SQRT( ( 1. + zfesatur * zkeq + zklight + consfe3(ji,jj,jk)/10**-6.3 - zkeq * trb(ji,jj,jk,jpfer) )**2       &
                   &              + 4. * ztfe * zkeq) ) / ( 2. * zkeq )
-               zFe3 (ji,jj,jk) = zFe3(ji,jj,jk) * 1E9
-               zFeL1(ji,jj,jk) = MAX( 0., trb(ji,jj,jk,jpfer) * 1E9 - zFe3(ji,jj,jk) )
+               zFeL1(ji,jj,jk) = MAX( 0., trb(ji,jj,jk,jpfer) - zFe3(ji,jj,jk) )
            END DO
          END DO
       END DO
-         !
-
+      !
+      plig(:,:,:) =  MAX( 0., ( zFeL1(:,:,:) / ( trb(:,:,:,jpfer) + rtrn ) ) )
+      !
       zdust = 0.         ! if no dust available
       DO jk = 1, jpkm1
@@ -113 +115 @@
                &         + fesol(ji,jj,jk,5) / zhplus )
                !
-               zfeequi = zFe3(ji,jj,jk) * 1E-9
-               zfecoll = 0.5 * zFeL1(ji,jj,jk) * 1E-9
+               zfecoll = 0.5 * zFeL1(ji,jj,jk)
                ! precipitation of Fe3+, creation of nanoparticles
-               precip(ji,jj,jk) = MAX( 0., ( zFe3(ji,jj,jk) * 1E-9 - fe3sol ) ) * kfep * xstep
+               precip(ji,jj,jk) = MAX( 0., ( zFe3(ji,jj,jk) - fe3sol ) ) * kfep * xstep * ( 1.0 - nitrfac(ji,jj,jk) ) 
+               ! Precipitation of Fe2+ due to oxidation by NO3 (Croot et al., 2019)
+               ! This occurs in anoxic waters only
+               precipno3(ji,jj,jk) = 2.0 * 130.0 * trb(ji,jj,jk,jpno3) * nitrfac(ji,jj,jk) * xstep * zFe3(ji,jj,jk)
                !
                ztrc   = ( trb(ji,jj,jk,jppoc) + trb(ji,jj,jk,jpgoc) + trb(ji,jj,jk,jpcal) + trb(ji,jj,jk,jpgsi) ) * 1.e6 
+               ztrc = MAX( rtrn, ztrc )
                IF( ln_dust )  zdust  = dust(ji,jj) / ( wdust / rday ) * tmask(ji,jj,jk)
-               IF (ln_ligand) THEN
-                  zxlam  = xlam1 * MAX( 1.E-3, EXP(-2 * etot(ji,jj,jk) / 10. ) * (1. - EXP(-2 * trb(ji,jj,jk,jpoxy) / 100.E-6 ) ))
-               ELSE
-                  zxlam  = xlam1 * 1.0
-               ENDIF
-               zlam1b = 3.e-5 + xlamdust * zdust + zxlam * ztrc
-               zscave = zfeequi * zlam1b * xstep
-
-               ! Compute the different ratios for scavenging of iron
-               ! to later allocate scavenged iron to the different organic pools
-               ! ---------------------------------------------------------
-               zdenom1 = zxlam * trb(ji,jj,jk,jppoc) / zlam1b
-               zdenom2 = zxlam * trb(ji,jj,jk,jpgoc) / zlam1b
-
-               ! Increased scavenging for very high iron concentrations found near the coasts 
-               ! due to increased lithogenic particles and let say it is unknown processes (precipitation, ...)
-               !  -----------------------------------------------------------
-               zlamfac = MAX( 0.e0, ( gphit(ji,jj) + 55.) / 30. )
-               zlamfac = MIN( 1.  , zlamfac )
-               zdep    = MIN( 1., 1000. / gdept_n(ji,jj,jk) )
-               zcoag   = 1E-4 * ( 1. - zlamfac ) * zdep * xstep * trb(ji,jj,jk,jpfer)
+               zxlam  = MAX( 1.E-3, (1. - EXP(-2 * trb(ji,jj,jk,jpoxy) / 100.E-6 ) ))
+               zlam1b = 3.e-5 + ( xlamdust * zdust + xlam1 * ztrc ) * zxlam
+               zscave = zFe3(ji,jj,jk) * zlam1b * xstep
 
                ! Compute the coagulation of colloidal iron. This parameterization 
@@ -146 +133 @@
                ! It requires certainly some more work as it is very poorly constrained.
                ! ----------------------------------------------------------------
-               zlam1a   = ( 0.369  * 0.3 * trb(ji,jj,jk,jpdoc) + 102.4  * trb(ji,jj,jk,jppoc) ) * xdiss(ji,jj,jk)    &
-                   &    + ( 114.   * 0.3 * trb(ji,jj,jk,jpdoc) )
+               zlam1a   = ( 12.0  * 0.3 * trb(ji,jj,jk,jpdoc) + 9.05  * trb(ji,jj,jk,jppoc) ) * xdiss(ji,jj,jk)    &
+                   &    + ( 2.49  * trb(ji,jj,jk,jppoc) )     &
+                   &    + ( 127.8 * 0.3 * trb(ji,jj,jk,jpdoc) + 725.7 * trb(ji,jj,jk,jppoc) )
                zaggdfea = zlam1a * xstep * zfecoll
                !
-               zlam1b   = 3.53E3 * trb(ji,jj,jk,jpgoc) * xdiss(ji,jj,jk)
+               zlam1b   = ( 1.94 * xdiss(ji,jj,jk) + 1.37 ) * trb(ji,jj,jk,jpgoc)
                zaggdfeb = zlam1b * xstep * zfecoll
                !
                tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zscave - zaggdfea - zaggdfeb &
-               &                     - zcoag - precip(ji,jj,jk)
-               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
+               &                     - precip(ji,jj,jk) - precipno3(ji,jj,jk)
+               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zscave * scaveff * trb(ji,jj,jk,jppoc) / ztrc
+               tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zscave * scaveff * trb(ji,jj,jk,jpgoc) / ztrc
+
+               ! Precipitated iron is supposed to be permanently lost.
+               ! Scavenged iron is supposed to be released back to seawater
+               ! when POM is solubilized. This is highly uncertain as probably
+               ! a significant part of it may be rescavenged back onto 
+               ! the particles. An efficiency factor is applied that is read
+               ! in the namelist. 
+               ! See for instance Tagliabue et al. (2019).
+               ! Aggregated FeL is considered as biogenic Fe as it 
+               ! probably remains  complexed when the particle is solubilized.
+               ! -------------------------------------------------------------
+               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zaggdfea
+               tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zaggdfeb
+               zscav3d(ji,jj,jk)   = zscave 
                zcoll3d(ji,jj,jk)   = zaggdfea + zaggdfeb
+               zprecip3d(ji,jj,jk) = precip(ji,jj,jk) + precipno3(ji,jj,jk)
                !
             END DO
@@ -175 +177 @@
 
                   ! Coagulation of ligands due to various processes (Brownian, shear, diff. sedimentation
-                  ! Coefficients are taken from the p4zagg
+                  ! Coefficients are taken from p4zagg
                   ! -------------------------------------------------------------------------------------
-                  zlam1a   = ( 0.369  * 0.3 * trb(ji,jj,jk,jpdoc) + 102.4  * trb(ji,jj,jk,jppoc) ) * xdiss(ji,jj,jk)    &
-                      &    + ( 114.   * 0.3 * trb(ji,jj,jk,jpdoc) )
+                  zlam1a   = ( 12.0  * 0.3 * trb(ji,jj,jk,jpdoc) + 9.05  * trb(ji,jj,jk,jppoc) ) * xdiss(ji,jj,jk)    &
+                      &    + ( 2.49  * trb(ji,jj,jk,jppoc) )     &
+                      &    + ( 127.8 * 0.3 * trb(ji,jj,jk,jpdoc) + 725.7 * trb(ji,jj,jk,jppoc) )
                   !
-                  zlam1b   = 3.53E3 *   trb(ji,jj,jk,jpgoc) * xdiss(ji,jj,jk)
-
+                  zlam1b   = ( 1.94 * xdiss(ji,jj,jk) + 1.37 ) * trb(ji,jj,jk,jpgoc)
                   ! 50% of the ligands are supposed to be in the colloidal size fraction
-                  zligco   = 0.5 * trn(ji,jj,jk,jplgw)
-                  zaggliga = zlam1a * xstep * zligco
+                  ! as for FeL
+                  zligco   = 0.5 * trb(ji,jj,jk,jplgw)
+                  zaggliga = zlam1a * xstep * zligco 
                   zaggligb = zlam1b * xstep * zligco
                   tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) - zaggliga - zaggligb
@@ -191 +194 @@
             END DO
          END DO
-         !
-         plig(:,:,:) =  MAX( 0., ( ( zFeL1(:,:,:) * 1E-9 ) / ( trb(:,:,:,jpfer) +rtrn ) ) )
          !
       ENDIF
@@ -206 +207 @@
          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("FEPREC") )  CALL iom_put("FEPREC" , zprecip3d(:,:,:)  * 1e9 * tmask(:,:,:) * zrfact2 )
          IF( iom_use("LGWCOLL"))  CALL iom_put("LGWCOLL", zlcoll3d(:,:,:) * 1e9 * tmask(:,:,:) * zrfact2 )
       ENDIF
@@ -234 +236 @@
       INTEGER ::   ios   ! Local integer 
       !!
-      NAMELIST/nampisfer/ ln_ligvar, xlam1, xlamdust, ligand, kfep 
+      NAMELIST/nampisfer/ ln_ligvar, xlam1, xlamdust, ligand, kfep, scaveff 
       !!----------------------------------------------------------------------
       !
@@ -258 +260 @@
          WRITE(numout,*) '      ligand concentration in the ocean         ligand       =', ligand
          WRITE(numout,*) '      rate constant for nanoparticle formation  kfep         =', kfep
+         WRITE(numout,*) '      Scavenged iron that is added to POFe      scaveff      =', scaveff
       ENDIF
       ! 

NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zflx.F90

@@ -359 +359 @@
       IF( ln_presatm ) THEN
          CALL fld_read( kt, 1, sf_patm )               !* input Patm provided at kt + 1/2
-         patm(:,:) = sf_patm(1)%fnow(:,:,1)                        ! atmospheric pressure
+         patm(:,:) = sf_patm(1)%fnow(:,:,1)/101325.0     ! atmospheric pressure
       ENDIF
       !

NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zint.F90

@@ -36 +36 @@
       !
       INTEGER  :: ji, jj                 ! dummy loop indices
-      REAL(wp) :: zvar                   ! local variable
+      REAL(wp) :: zrum, zcodel, zargu, zvar
       !!---------------------------------------------------------------------
       !
@@ -44 +44 @@
       ! -------------------------------------------
       ! Generic temperature dependence (Eppley, 1972)
-      tgfunc (:,:,:) = EXP( 0.063913 * tsn(:,:,:,jp_tem) )
+      tgfunc (:,:,:) = EXP( 0.0631 * tsn(:,:,:,jp_tem) )
       ! Temperature dependence of mesozooplankton (Buitenhuis et al. (2005))
-      tgfunc2(:,:,:) = EXP( 0.07608  * tsn(:,:,:,jp_tem) )
-      ! Temperature dependence of picophytoplankton (Stawiarsky et al., 2016)
-      tgfunc3(:,:,:) = EXP( 0.0825   * tsn(:,:,:,jp_tem) )
+      tgfunc2(:,:,:) = EXP( 0.0761  * tsn(:,:,:,jp_tem) )
 
       ! Computation of the silicon dependant half saturation  constant for silica uptake
@@ -69 +67 @@
       ENDIF
       !
+      ! compute the day length depending on latitude and the day
+      ! Astronomical parameterization taken from HAMOCC3
+      zrum = REAL( nday_year - 80, wp ) / REAL( nyear_len(1), wp )
+      zcodel = ASIN(  SIN( zrum * rpi * 2._wp ) * SIN( rad * 23.5_wp )  )
+
+      ! day length in hours
+      strn(:,:) = 0.
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zargu = TAN( zcodel ) * TAN( gphit(ji,jj) * rad )
+            zargu = MAX( -1., MIN(  1., zargu ) )
+            strn(ji,jj) = MAX( 0.0, 24. - 2. * ACOS( zargu ) / rad / 15. )
+         END DO
+      END DO
+      !
       IF( ln_timing )   CALL timing_stop('p4z_int')
       !

NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zlim.F90

@@ -63 +63 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xnanofer   !: Limitation of Fe uptake by nanophyto
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xdiatfer   !: Limitation of Fe uptake by diatoms
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xqfuncfecd, xqfuncfecn
 
    ! Coefficient for iron limitation following Flynn and Hipkin (1999)
@@ -90 +91 @@
       !
       INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   zlim1, zlim2, zlim3, zlim4, zno3, zferlim, zcoef
+      REAL(wp) ::   zlim1, zlim2, zlim3, zlim4, zcoef
       REAL(wp) ::   z1_trbdia, z1_trbphy, ztem1, ztem2, zetot1, zetot2
-      REAL(wp) ::   zdenom, zratio, zironmin
+      REAL(wp) ::   zdenom, zratio, zironmin, zbactno3, zbactnh4
       REAL(wp) ::   zconc1d, zconc1dnh4, zconc0n, zconc0nnh4   
       REAL(wp) ::   fananof, fadiatf, znutlim, zfalim
+      REAL(wp) ::   znutlimtot, zlimno3, zlimnh4, zbiron
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('p4z_lim')
       !
-      sizena(:,:,:) = 0.0  ;  sizeda(:,:,:) = 0.0
+      sizena(:,:,:) = 1.0  ;  sizeda(:,:,:) = 1.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 a variable Ks of diatoms taking into account
                ! that increasing biomass is made of generally bigger cells
@@ -134 +127 @@
 
                ! Nanophytoplankton
-               znutlim = biron(ji,jj,jk) / concnfe(ji,jj,jk)
+               zbiron = ( 75.0 * ( 1.0 - plig(ji,jj,jk) ) + plig(ji,jj,jk) ) * biron(ji,jj,jk)
+               znutlim = zbiron / concnfe(ji,jj,jk)
                fananof = MAX(0.01, MIN(0.99, 1. / ( SQRT(znutlim) + 1.) ) )
 
                ! Diatoms
-               znutlim = biron(ji,jj,jk) / concdfe(ji,jj,jk)
+               znutlim = zbiron / concdfe(ji,jj,jk)
                fadiatf = MAX(0.01, MIN(0.99, 1. / ( SQRT(znutlim) + 1.) ) )
 
@@ -144 +138 @@
                ! heterotrophic bacteria
                ! -------------------------------------------------
-               zdenom = 1. /  ( concbno3 * concbnh4 + concbnh4 * trb(ji,jj,jk,jpno3) + concbno3 * trb(ji,jj,jk,jpnh4) )
-               xnanono3(ji,jj,jk) = trb(ji,jj,jk,jpno3) * concbnh4 * zdenom
-               xnanonh4(ji,jj,jk) = trb(ji,jj,jk,jpnh4) * concbno3 * zdenom
+               zlimnh4 = trb(ji,jj,jk,jpnh4) / ( concbno3 + trb(ji,jj,jk,jpnh4) )
+               zlimno3 = trb(ji,jj,jk,jpno3) / ( concbno3 + trb(ji,jj,jk,jpno3) )
+               znutlimtot = ( trb(ji,jj,jk,jpnh4) + trb(ji,jj,jk,jpno3) ) / ( concbno3 + trb(ji,jj,jk,jpnh4) + trb(ji,jj,jk,jpno3) )
+               zbactnh4 = znutlimtot * 5.0 * zlimnh4 / ( zlimno3 + 5.0 * zlimnh4 + rtrn )
+               zbactno3 = znutlimtot * zlimno3 / ( zlimno3 + 5.0 * zlimnh4 + rtrn )
                !
-               zlim1    = xnanono3(ji,jj,jk) + xnanonh4(ji,jj,jk)
+               zlim1    = zbactno3 + zbactnh4
                zlim2    = trb(ji,jj,jk,jppo4) / ( trb(ji,jj,jk,jppo4) + concbnh4 )
                zlim3    = trb(ji,jj,jk,jpfer) / ( concbfe + trb(ji,jj,jk,jpfer) )
@@ -166 +162 @@
                ! Limitation of Fe uptake (Quota formalism)
                zfalim = (1.-fananof) / fananof
-               xnanofer(ji,jj,jk) = (1. - fananof) * biron(ji,jj,jk) / ( biron(ji,jj,jk) + zfalim * concnfe(ji,jj,jk) )
+               xnanofer(ji,jj,jk) = (1. - fananof) * zbiron / ( zbiron + zfalim * concnfe(ji,jj,jk) )
 
                ! Limitation of nanophytoplankton growth
-               zdenom = 1. /  ( zconc0n * zconc0nnh4 + zconc0nnh4 * trb(ji,jj,jk,jpno3) + zconc0n * trb(ji,jj,jk,jpnh4) )
-               xnanono3(ji,jj,jk) = trb(ji,jj,jk,jpno3) * zconc0nnh4 * zdenom
-               xnanonh4(ji,jj,jk) = trb(ji,jj,jk,jpnh4) * zconc0n    * zdenom
+               zlimnh4 = trb(ji,jj,jk,jpnh4) / ( zconc0n + trb(ji,jj,jk,jpnh4) )
+               zlimno3 = trb(ji,jj,jk,jpno3) / ( zconc0n + trb(ji,jj,jk,jpno3) )
+               znutlimtot = ( trb(ji,jj,jk,jpnh4) + trb(ji,jj,jk,jpno3) ) / ( zconc0n + trb(ji,jj,jk,jpnh4) + trb(ji,jj,jk,jpno3) )
+               xnanonh4(ji,jj,jk) = znutlimtot * 5.0 * zlimnh4 / ( zlimno3 + 5.0 * zlimnh4 + rtrn )
+               xnanono3(ji,jj,jk) = znutlimtot * zlimno3 / ( zlimno3 + 5.0 * zlimnh4 + rtrn )
                !
                zlim1    = xnanono3(ji,jj,jk) + xnanonh4(ji,jj,jk)
@@ -181 +179 @@
                ! proposed by Flynn and Hipkin (1999)
                zironmin = xcoef1 * trb(ji,jj,jk,jpnch) * z1_trbphy + xcoef2 * zlim1 + xcoef3 * xnanono3(ji,jj,jk)
+               xqfuncfecn(ji,jj,jk) = zironmin + qnfelim
                zlim3    = MAX( 0.,( zratio - zironmin ) / qnfelim )
                xnanopo4(ji,jj,jk) = zlim2
@@ -190 +189 @@
                ! Limitation of Fe uptake (Quota formalism)
                zfalim = (1.-fadiatf) / fadiatf
-               xdiatfer(ji,jj,jk) = (1. - fadiatf) * biron(ji,jj,jk) / ( biron(ji,jj,jk) + zfalim * concdfe(ji,jj,jk) )
+               xdiatfer(ji,jj,jk) = (1. - fadiatf) * zbiron / ( zbiron + zfalim * concdfe(ji,jj,jk) )
 
                ! Limitation of diatoms growth
-               zdenom   = 1. / ( zconc1d * zconc1dnh4 + zconc1dnh4 * trb(ji,jj,jk,jpno3) + zconc1d * trb(ji,jj,jk,jpnh4) )
-               xdiatno3(ji,jj,jk) = trb(ji,jj,jk,jpno3) * zconc1dnh4 * zdenom
-               xdiatnh4(ji,jj,jk) = trb(ji,jj,jk,jpnh4) * zconc1d    * zdenom
+               zlimnh4 = trb(ji,jj,jk,jpnh4) / ( zconc1d + trb(ji,jj,jk,jpnh4) )
+               zlimno3 = trb(ji,jj,jk,jpno3) / ( zconc1d + trb(ji,jj,jk,jpno3) )
+               znutlimtot = ( trb(ji,jj,jk,jpnh4) + trb(ji,jj,jk,jpno3) ) / ( zconc1d + trb(ji,jj,jk,jpnh4) + trb(ji,jj,jk,jpno3) )
+               xdiatnh4(ji,jj,jk) = znutlimtot * 5.0 * zlimnh4 / ( zlimno3 + 5.0 * zlimnh4 + rtrn ) 
+               xdiatno3(ji,jj,jk) = znutlimtot * zlimno3 / ( zlimno3 + 5.0 * zlimnh4 + rtrn )
                !
                zlim1    = xdiatno3(ji,jj,jk) + xdiatnh4(ji,jj,jk)
@@ -206 +207 @@
                ! proposed by Flynn and Hipkin (1999)
                zironmin = xcoef1 * trb(ji,jj,jk,jpdch) * z1_trbdia + xcoef2 * zlim1 + xcoef3 * xdiatno3(ji,jj,jk)
+               xqfuncfecd(ji,jj,jk) = zironmin + qdfelim
                zlim4    = MAX( 0., ( zratio - zironmin ) / qdfelim )
                xdiatpo4(ji,jj,jk) = zlim2
@@ -238 +240 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zlim1 =  ( trb(ji,jj,jk,jpno3) * concnnh4 + trb(ji,jj,jk,jpnh4) * concnno3 )    &
-                  &   / ( concnno3 * concnnh4 + concnnh4 * trb(ji,jj,jk,jpno3) + concnno3 * trb(ji,jj,jk,jpnh4) ) 
+               zlim1  = xnanonh4(ji,jj,jk) + xnanono3(ji,jj,jk) 
                zlim2  = trb(ji,jj,jk,jppo4) / ( trb(ji,jj,jk,jppo4) + concnnh4 )
-               zlim3  = trb(ji,jj,jk,jpfer) / ( trb(ji,jj,jk,jpfer) +  5.E-11   )
-               ztem1  = MAX( 0., tsn(ji,jj,jk,jp_tem) )
+               zlim3  = trb(ji,jj,jk,jpfer) / ( trb(ji,jj,jk,jpfer) +  6.E-11   )
+               ztem1  = MAX( 0., tsn(ji,jj,jk,jp_tem) + 1.8)
                ztem2  = tsn(ji,jj,jk,jp_tem) - 10.
                zetot1 = MAX( 0., etot_ndcy(ji,jj,jk) - 1.) / ( 4. + etot_ndcy(ji,jj,jk) ) 
-               zetot2 = 30. / ( 30. + etot_ndcy(ji,jj,jk) ) 
+               zetot2 = 30. / ( 30.0 + etot_ndcy(ji,jj,jk) )
 
                xfracal(ji,jj,jk) = caco3r * MIN( zlim1, zlim2, zlim3 )                  &
@@ -375 +376 @@
          &      xlimbac (jpi,jpj,jpk), xlimbacl(jpi,jpj,jpk),       &
          &      concnfe (jpi,jpj,jpk), concdfe (jpi,jpj,jpk),       &
+         &      xqfuncfecn(jpi,jpj,jpk), xqfuncfecd(jpi,jpj,jpk),   &
          &      xlimsi  (jpi,jpj,jpk), STAT=p4z_lim_alloc )
       !

NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zlys.F90

@@ -57 +57 @@
       !
       INTEGER  ::   ji, jj, jk, jn
-      REAL(wp) ::   zdispot, zfact, zcalcon
+      REAL(wp) ::   zdispot, zfact, zcalcon, zkd
       REAL(wp) ::   zomegaca, zexcess, zexcess0
       CHARACTER (len=25) ::   charout
@@ -65 +65 @@
       IF( ln_timing )  CALL timing_start('p4z_lys')
       !
-      zco3    (:,:,:) = 0.
-      zcaldiss(:,:,:) = 0.
+      zco3    (:,:,:) = 0._wp
+      zcaldiss(:,:,:) = 0._wp     ;     zco3sat(:,:,:) = 0._wp
       zhinit  (:,:,:) = hi(:,:,:) * 1000. / ( rhop(:,:,:) + rtrn )
       !
@@ -106 +106 @@
                excess(ji,jj,jk) = 1._wp - zomegaca
                zexcess0 = MAX( 0., excess(ji,jj,jk) )
-               zexcess  = zexcess0**nca
+               IF (zomegaca < 0.8) THEN
+                  zexcess  = zexcess0**nca
+                 ! AMOUNT CACO3 THAT RE-ENTERS SOLUTION
+                  zdispot = kdca * zexcess * trb(ji,jj,jk,jpcal)
+               ELSE
+                  zkd = kdca * 0.1**(nca - 0.4)
+                  zexcess  = zexcess0**0.4
+                  zdispot = zkd * zexcess * trb(ji,jj,jk,jpcal)
+               ENDIF
 
-               ! AMOUNT CACO3 THAT RE-ENTERS SOLUTION
-               zdispot = kdca * zexcess * trb(ji,jj,jk,jpcal)
                ! CHANGE OF [CO3--] , [ALK], PARTICULATE [CACO3],
                ! AND [SUM(CO2)] DUE TO CACO3 DISSOLUTION/PRECIPITATION

NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zmeso.F90

@@ -89 +89 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zgrarem, zgraref, zgrapoc, zgrapof, zgrabsi
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zgramigrem, zgramigref, zgramigpoc, zgramigpof
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zstrn, zgramigbsi
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zgramigbsi
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zw3d, zz2ligprod
       !!---------------------------------------------------------------------
@@ -236 +236 @@
                zgraztotn = zgrazdc * quotad(ji,jj,jk) + zgrazz + zgraznc * quotan(ji,jj,jk)   &
                &   + zgrazpoc + zgrazffep + zgrazffeg
-               zgraztotf = zgrazdf + zgraznf + zgrazz * ferat3 + zgrazpof + zgrazfffp + zgrazfffg
+               zgraztotf = zgrazdf + zgraznf + zgrazz * feratz + zgrazpof + zgrazfffp + zgrazfffg
 
                ! Total grazing ( grazing by microzoo is already computed in p4zmicro )
@@ -251 +251 @@
                zgrasratf =  ( zgraztotf + rtrn )/ ( zgraztotc + rtrn )
                zgrasratn =  ( zgraztotn + rtrn )/ ( zgraztotc + rtrn )
-               zepshert  = MIN( 1., zgrasratn, zgrasratf / ferat3)
+               zepshert  = MIN( 1., zgrasratn, zgrasratf / feratm)
                zbeta     = MAX(0., (epsher2 - epsher2min) )
                ! Food quantity deprivation of GGE
@@ -298 +298 @@
                zprcaca = part2 * zprcaca
                tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) + zgrazcal - zprcaca
-               tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) - 2. * ( zgrazcal + zprcaca )
+               tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + 2. * ( zgrazcal - zprcaca )
                tra(ji,jj,jk,jpcal) = tra(ji,jj,jk,jpcal) - zgrazcal + zprcaca
  
@@ -305 +305 @@
                zgrarem(ji,jj,jk) = zgraztotc * ( 1. - zepsherv - unass2 ) &
                &         + ( 1. - epsher2 - unass2 ) / ( 1. - epsher2 ) * ztortz
-               zgraref(ji,jj,jk) = zgraztotc * MAX( 0. , ( 1. - unass2 ) * zgrasratf - ferat3 * zepsherv )    &
-               &         + ferat3 * ( ( 1. - epsher2 - unass2 ) /( 1. - epsher2 ) * ztortz )
+               zgraref(ji,jj,jk) = zgraztotc * MAX( 0. , ( 1. - unass2 ) * zgrasratf - feratm * zepsherv )    &
+               &         + feratm * ( ( 1. - epsher2 - unass2 ) /( 1. - epsher2 ) * ztortz )
                zgrapoc(ji,jj,jk) = zgraztotc * unass2 + zmortzgoc
-               zgrapof(ji,jj,jk) = zgraztotf * unass2 + ferat3 * zmortzgoc
+               zgrapof(ji,jj,jk) = zgraztotf * unass2 + feratm * zmortzgoc
             END DO 
          END DO
@@ -320 +320 @@
          ALLOCATE( zgramigrem(jpi,jpj), zgramigref(jpi,jpj), zgramigpoc(jpi,jpj), zgramigpof(jpi,jpj) )
          ALLOCATE( zgramigbsi(jpi,jpj) )
-         ALLOCATE( zstrn(jpi,jpj) )
          zgramigrem(:,:) = 0.0    ;   zgramigref(:,:) = 0.0
          zgramigpoc(:,:) = 0.0    ;   zgramigpof(:,:) = 0.0
          zgramigbsi(:,:) = 0.0
-
-         ! compute the day length depending on latitude and the day
-         zrum = REAL( nday_year - 80, wp ) / REAL( nyear_len(1), wp )
-         zcodel = ASIN(  SIN( zrum * rpi * 2._wp ) * SIN( rad * 23.5_wp )  )
-
-         ! day length in hours
-         zstrn(:,:) = 0.
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zargu = TAN( zcodel ) * TAN( gphit(ji,jj) * rad )
-               zargu = MAX( -1., MIN(  1., zargu ) )
-               zstrn(ji,jj) = MAX( 0.0, 24. - 2. * ACOS( zargu ) / rad / 15. )
-               zstrn(ji,jj) = MIN(0.75, MAX( 0.25, zstrn(ji,jj) / 24.) )
-            END DO
-         END DO
 
         ! Compute the amount of materials that will go into vertical migration
@@ -347 +331 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  zmigreltime = (1. - zstrn(ji,jj))
+                  zmigreltime = (1. - strn(ji,jj))
                   IF ( gdept_n(ji,jj,jk) <= heup(ji,jj) ) THEN
                      zgramigrem(ji,jj) = zgramigrem(ji,jj) + xfracmig * zgrarem(ji,jj,jk) * (1. - zmigreltime )    &
@@ -389 +373 @@
          ! ------------------------------
          DEALLOCATE( zgramigrem, zgramigref, zgramigpoc, zgramigpof, zgramigbsi )
-         DEALLOCATE( zstrn )
-
       ! End of the ln_dvm_meso part
       ENDIF

NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zmicro.F90

@@ -187 +187 @@
                zgrasratf = ( zgraztotf + rtrn ) / ( zgraztotc + rtrn )
                zgrasratn = ( zgraztotn + rtrn ) / ( zgraztotc + rtrn )
-               zepshert  =  MIN( 1., zgrasratn, zgrasratf / ferat3)
+               zepshert  =  MIN( 1., zgrasratn, zgrasratf / feratz)
                zbeta     = MAX(0., (epsher - epshermin) )
                ! Food quantity deprivation of the GGE
@@ -196 +196 @@
                zepsherv  = zepsherf * zepshert * zepsherq
                ! Excretion of Fe
-               zgrafer   = zgraztotc * MAX( 0. , ( 1. - unass ) * zgrasratf - ferat3 * zepsherv ) 
+               zgrafer   = zgraztotc * MAX( 0. , ( 1. - unass ) * zgrasratf - feratz * zepsherv ) 
                ! Excretion of C, N, P
                zgrarem   = zgraztotc * ( 1. - zepsherv - unass )
@@ -236 +236 @@
                prodpoc(ji,jj,jk) = prodpoc(ji,jj,jk) + zmortz
                conspoc(ji,jj,jk) = conspoc(ji,jj,jk) - zgrazpoc
-               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + ferat3 * zmortz - zgrazpof
+               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + feratz * zmortz - zgrazpof
                !
                ! Calcite remineralization due to zooplankton activity

NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zmort.F90

@@ -25 +25 @@
    REAL(wp), PUBLIC ::   wchln    !: Quadratic mortality rate of nanophytoplankton
    REAL(wp), PUBLIC ::   wchld    !: Quadratic mortality rate of diatoms
-   REAL(wp), PUBLIC ::   wchldm   !: Maximum quadratic mortality rate of diatoms
    REAL(wp), PUBLIC ::   mpratn   !: Linear mortality rate of nanophytoplankton
    REAL(wp), PUBLIC ::   mpratd   !: Linear mortality rate of diatoms
@@ -63 +62 @@
       !!---------------------------------------------------------------------
       INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   zsizerat, zcompaph
+      REAL(wp) ::   zcompaph
       REAL(wp) ::   zfactfe, zfactch, zprcaca, zfracal
-      REAL(wp) ::   ztortp , zrespp , zmortp 
+      REAL(wp) ::   ztortp , zrespp , zmortp, zlim1, zlim2 
       CHARACTER (len=25) ::   charout
       !!---------------------------------------------------------------------
@@ -76 +75 @@
             DO ji = 1, jpi
                zcompaph = MAX( ( trb(ji,jj,jk,jpphy) - 1e-8 ), 0.e0 )
-               !     When highly limited by macronutrients, very small cells 
-               !     dominate the community. As a consequence, aggregation
-               !     due to turbulence is negligible. Mortality is also set
-               !     to 0
-               zsizerat = MIN(1., MAX( 0., (quotan(ji,jj,jk) - 0.2) / 0.3) ) * trb(ji,jj,jk,jpphy)
-
                ! Quadratic mortality of nano due to aggregation during
                ! blooms (Doney et al. 1996)
                ! -----------------------------------------------------
-               zrespp = wchln * 1.e6 * xstep * xdiss(ji,jj,jk) * zcompaph * zsizerat 
+               zlim2   = xlimphy(ji,jj,jk) * xlimphy(ji,jj,jk)
+               zlim1   = 0.25 * ( 1. - zlim2 ) / ( 0.25 + zlim2 ) * trb(ji,jj,jk,jpphy)
+               zrespp  = wchln * 1.e6 * xstep * zlim1 * xdiss(ji,jj,jk) * zcompaph
 
                ! Phytoplankton linear mortality
@@ -91 +86 @@
                ! extinction of nanophyto in highly limited areas
                ! ----------------------------------------------------
-               ztortp = mpratn * xstep * zcompaph / ( xkmort + trb(ji,jj,jk,jpphy) ) * zsizerat
+               ztortp = mpratn * xstep * zcompaph / ( xkmort + trb(ji,jj,jk,jpphy) ) * trb(ji,jj,jk,jpphy)
 
                zmortp = zrespp + ztortp
@@ -170 +165 @@
                zlim2   = xlimdia(ji,jj,jk) * xlimdia(ji,jj,jk)
                zlim1   = 0.25 * ( 1. - zlim2 ) / ( 0.25 + zlim2 ) 
-               zrespp2 = 1.e6 * xstep * (  wchld + wchldm * zlim1 ) * xdiss(ji,jj,jk) * zcompadi * trb(ji,jj,jk,jpdia)
+               zrespp2 = 1.e6 * xstep * wchld * zlim1 * xdiss(ji,jj,jk) * zcompadi * trb(ji,jj,jk,jpdia)
 
                ! Phytoplankton linear mortality
@@ -228 +223 @@
       INTEGER ::   ios   ! Local integer
       !
-      NAMELIST/namp4zmort/ wchln, wchld, wchldm, mpratn, mpratd
+      NAMELIST/namp4zmort/ wchln, wchld, mpratn, mpratd
       !!----------------------------------------------------------------------
       !
@@ -249 +244 @@
          WRITE(numout,*) '      quadratic mortality of phytoplankton        wchln  =', wchln
          WRITE(numout,*) '      maximum quadratic mortality of diatoms      wchld  =', wchld
-         WRITE(numout,*) '      maximum quadratic mortality of diatoms      wchldm =', wchldm
          WRITE(numout,*) '      phytoplankton mortality rate                mpratn =', mpratn
          WRITE(numout,*) '      Diatoms mortality rate                      mpratd =', mpratd

NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zopt.F90

@@ -67 +67 @@
       REAL(wp), DIMENSION(jpi,jpj    ) :: zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4
       REAL(wp), DIMENSION(jpi,jpj    ) :: zqsr100, zqsr_corr
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpar, ze0, ze1, ze2, ze3, zchl3d
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze0, ze1, ze2, ze3, zchl3d
       !!---------------------------------------------------------------------
       !
@@ -206 +206 @@
                  heup(ji,jj) = gdepw_n(ji,jj,jk+1)     ! Euphotic layer depth
               ENDIF
-              IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >= 0.50 )  THEN
+              IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >= 0.10 )  THEN
                  heup_01(ji,jj) = gdepw_n(ji,jj,jk+1)  ! Euphotic layer depth (light level definition)
               ENDIF
@@ -235 +235 @@
       END DO
       !
-      emoy(:,:,:) = etot(:,:,:)       ! PAR
-      zpar(:,:,:) = etot_ndcy(:,:,:)  ! diagnostic : PAR with no diurnal cycle 
+      emoy(:,:,:)  = etot(:,:,:)       ! PAR
+      etotm(:,:,:) = etot_ndcy(:,:,:)  ! diagnostic : PAR with no diurnal cycle 
       !
       DO jk = 1, nksrp
@@ -244 +244 @@
                   z1_dep = 1. / ( zdepmoy(ji,jj) + rtrn )
                   emoy (ji,jj,jk) = zetmp1(ji,jj) * z1_dep
-                  zpar (ji,jj,jk) = zetmp2(ji,jj) * z1_dep
+                  etotm(ji,jj,jk) = zetmp2(ji,jj) * z1_dep
                ENDIF
             END DO
@@ -313 +313 @@
         IF( knt == nrdttrc ) THEN
            IF( iom_use( "Heup"  ) ) CALL iom_put( "Heup" , heup(:,:  ) * tmask(:,:,1) )  ! euphotic layer deptht
-           IF( iom_use( "PARDM" ) ) CALL iom_put( "PARDM", zpar(:,:,:) * tmask(:,:,:) )  ! Photosynthetically Available Radiation
+           IF( iom_use( "PARDM" ) ) CALL iom_put( "PARDM", etotm(:,:,:) * tmask(:,:,:) )  ! Photosynthetically Available Radiation
            IF( iom_use( "PAR"   ) ) CALL iom_put( "PAR"  , emoy(:,:,:) * tmask(:,:,:) )  ! Photosynthetically Available Radiation
         ENDIF

NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zprod.F90

@@ -67 +67 @@
       INTEGER  ::   ji, jj, jk
       REAL(wp) ::   zsilfac, znanotot, zdiattot, zconctemp, zconctemp2
-      REAL(wp) ::   zratio, zmax, zsilim, ztn, zadap, zlim, zsilfac2, zsiborn
+      REAL(wp) ::   zratio, zmax, zsilim, ztn, zadap, zlim, zsiborn
       REAL(wp) ::   zprod, zproreg, zproreg2, zprochln, zprochld
       REAL(wp) ::   zdocprod, zpislopen, zpisloped, zfact
-      REAL(wp) ::   zratiosi, zmaxsi, zlimfac, zsizetmp
-      REAL(wp) ::   zrum, zcodel, zargu, zval, zfeup, chlcnm_n, chlcdm_n
+      REAL(wp) ::   zratiosi, zmaxsi, zlimfac, zsizetmp, zfecnm, zfecdm
+      REAL(wp) ::   zrum, zcodel, zargu, zval, zfeup
       CHARACTER (len=25) :: charout
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zw2d
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d
-      REAL(wp), DIMENSION(jpi,jpj    ) :: zstrn
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprmaxn,zprmaxd
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpislopeadn, zpislopeadd, zysopt  
@@ -94 +93 @@
       zprbio  (:,:,:) = 0._wp ; zprchld (:,:,:) = 0._wp ; zprchln (:,:,:) = 0._wp 
       zmxl_fac(:,:,:) = 0._wp ; zmxl_chl(:,:,:) = 0._wp 
+      consfe3 (:,:,:) = 0._wp
 
       ! Computation of the maximimum production. Based on a Q10 description
@@ -100 +100 @@
       zprmaxn(:,:,:) = 0.8_wp * r1_rday * tgfunc(:,:,:)
       zprmaxd(:,:,:) = zprmaxn(:,:,:)
-
-      ! compute the day length depending on latitude and the day
-      ! Astronomical parameterization taken from HAMOCC3
-      zrum = REAL( nday_year - 80, wp ) / REAL( nyear_len(1), wp )
-      zcodel = ASIN(  SIN( zrum * rpi * 2._wp ) * SIN( rad * 23.5_wp )  )
-
-      ! day length in hours
-      zstrn(:,:) = 0.
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zargu = TAN( zcodel ) * TAN( gphit(ji,jj) * rad )
-            zargu = MAX( -1., MIN(  1., zargu ) )
-            zstrn(ji,jj) = MAX( 0.0, 24. - 2. * ACOS( zargu ) / rad / 15. )
-         END DO
-      END DO
 
       ! Impact of the day duration and light intermittency on phytoplankton growth
@@ -127 +112 @@
             DO ji = 1, jpi
                IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
-                  zval = MAX( 1., zstrn(ji,jj) )
+                  zval = MAX( 1., strn(ji,jj) )
                   IF( gdepw_n(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
                      zval = zval * MIN(1., heup_01(ji,jj) / ( hmld(ji,jj) + rtrn ))
@@ -225 +210 @@
                    ! Si/C of diatoms
                    ! ------------------------
-                   ! Si/C increases with iron stress and silicate availability (zsilfac)
+                   ! Si/C increases with iron stress and silicate availability
                    ! Si/C is arbitrariliy increased for very high Si concentrations
-                   ! to mimic the very high ratios observed in the Southern Ocean (zsilfac2)
+                   ! to mimic the very high ratios observed in the Southern Ocean (zsilfac)
                    ! A parameterization derived from Flynn (2003) is used for the control
                    ! when Si is not limiting which is similar to the parameterisation
@@ -236 +221 @@
                  zsiborn = trb(ji,jj,1,jpsil)**3
                  IF (gphit(ji,jj) < -30.0 ) THEN
-                   zsilfac2 = 1. + 2. * zsiborn / ( zsiborn + xksi2**3 )
+                   zsilfac = 1. + 2. * zsiborn / ( zsiborn + xksi2**3 )
                  ELSE
-                   zsilfac2 = 1. +      zsiborn / ( zsiborn + xksi2**3 )
+                   zsilfac = 1. +      zsiborn / ( zsiborn + xksi2**3 )
                  ENDIF
-                 zratiosi = 1.0 - trb(ji,jj,jk,jpdsi) / ( trb(ji,jj,jk,jpdia) + rtrn ) / ( zsilfac2 * grosip * 3.0 + rtrn )
+                 zratiosi = 1.0 - trb(ji,jj,jk,jpdsi) / ( trb(ji,jj,jk,jpdia) + rtrn ) / ( zsilfac * grosip * 3.0 + rtrn )
                  zratiosi = MAX(0., MIN(1.0, zratiosi) )
                  zmaxsi  = (1.0 + 0.1**4) * zratiosi**4 / ( zratiosi**4 + 0.1**4 )
                  IF ( xlimsi(ji,jj,jk) /= xlimdia(ji,jj,jk) ) THEN
-                    zysopt(ji,jj,jk) = zlim * zsilfac2 * grosip * 1.0 * zmaxsi
+                    zysopt(ji,jj,jk) = zlim * zsilfac * grosip * 1.0 * zmaxsi
                  ELSE
-                    zysopt(ji,jj,jk) = zlim * zsilfac2 * grosip * 1.0 * zsilim**0.7 * zmaxsi
+                    zysopt(ji,jj,jk) = zlim * zsilfac * grosip * 1.0 * zsilim**0.7 * zmaxsi
                  ENDIF
               ENDIF
@@ -291 +276 @@
                   ! do not suggest it for accimated cells. Uptake is
                   ! downregulated when the quota is close to the maximum quota
-                  zratio = 1.0 - MIN(1.0,trb(ji,jj,jk,jpnfe) / ( trb(ji,jj,jk,jpphy) * fecnm + rtrn ) )
+                  zfecnm = xqfuncfecn(ji,jj,jk) + ( fecnm - xqfuncfecn(ji,jj,jk) ) * ( xnanono3(ji,jj,jk) + xnanonh4(ji,jj,jk) )
+                  zratio = 1.0 - MIN(1.0,trb(ji,jj,jk,jpnfe) / ( trb(ji,jj,jk,jpphy) * zfecnm + rtrn ) )
                   zmax   = MAX( 0., MIN( 1.0, zratio**2/ (0.05**2+zratio**2) ) ) 
-                  zprofen(ji,jj,jk) = fecnm * zprmaxn(ji,jj,jk) * ( 1.0 - fr_i(ji,jj) )  &
+                  zprofen(ji,jj,jk) = zfecnm * zprmaxn(ji,jj,jk) * ( 1.0 - fr_i(ji,jj) )  &
                   &          * (1. + 0.8 * xnanono3(ji,jj,jk) / ( rtrn + xnanono3(ji,jj,jk)  &
                   &          + xnanonh4(ji,jj,jk) ) * (1. - xnanofer(ji,jj,jk) ) )   &
@@ -317 +303 @@
                   ! do not suggest it for accimated cells. Uptake is
                   ! downregulated when the quota is close to the maximum quota
-                  zratio = 1.0 - MIN(1.0, trb(ji,jj,jk,jpdfe) / ( trb(ji,jj,jk,jpdia) * fecdm + rtrn ) )
+                  zfecdm = xqfuncfecd(ji,jj,jk) + ( fecdm - xqfuncfecd(ji,jj,jk) ) * ( xdiatno3(ji,jj,jk) + xdiatnh4(ji,jj,jk) )
+                  zratio = 1.0 - MIN(1.0, trb(ji,jj,jk,jpdfe) / ( trb(ji,jj,jk,jpdia) * zfecdm + rtrn ) )
                   zmax   = MAX( 0., MIN( 1.0, zratio**2/ (0.05**2+zratio**2) ) ) 
-                  zprofed(ji,jj,jk) = fecdm * zprmaxd(ji,jj,jk) * (1.0 - fr_i(ji,jj) )  &
+                  zprofed(ji,jj,jk) = zfecdm * zprmaxd(ji,jj,jk) * (1.0 - fr_i(ji,jj) )  &
                   &          * (1. + 0.8 * xdiatno3(ji,jj,jk) / ( rtrn + xdiatno3(ji,jj,jk)  &
                   &          + xdiatnh4(ji,jj,jk) ) * (1. - xdiatfer(ji,jj,jk) ) )   &
@@ -339 +326 @@
                   zprod    = rday * zprorcan(ji,jj,jk) * zprchln(ji,jj,jk) * xlimphy(ji,jj,jk)
                   zprochln = chlcmin * 12. * zprorcan (ji,jj,jk)
-
-                  ! The maximum reachable Chl quota is modulated by temperature
-                  ! following Geider (1987)
-                  chlcnm_n   = MIN ( chlcnm, ( chlcnm / (1. - 1.14 / 43.4 *tsn(ji,jj,jk,jp_tem))) * (1. - 1.14 / 43.4 * 20.))
-                  zprochln = zprochln + (chlcnm_n-chlcmin) * 12. * zprod / &
+                  zprochln = zprochln + (chlcnm - chlcmin) * 12. * zprod / &
                                         & (  zpislopeadn(ji,jj,jk) * znanotot +rtrn)
 
@@ -350 +333 @@
                   zprod    = rday * zprorcad(ji,jj,jk) * zprchld(ji,jj,jk) * xlimdia(ji,jj,jk)
                   zprochld = chlcmin * 12. * zprorcad(ji,jj,jk)
-
-                  ! The maximum reachable Chl quota is modulated by temperature
-                  ! following Geider (1987)
-                  chlcdm_n   = MIN ( chlcdm, ( chlcdm / (1. - 1.14 / 43.4 * tsn(ji,jj,jk,jp_tem))) * (1. - 1.14 / 43.4 * 20.))
-                  zprochld = zprochld + (chlcdm_n-chlcmin) * 12. * zprod / &
+                  zprochld = zprochld + (chlcdm - chlcmin) * 12. * zprod / &
                                         & ( zpislopeadd(ji,jj,jk) * zdiattot +rtrn )
 
@@ -386 +365 @@
                  &                   + ( o2ut + o2nit ) * ( zpronewn(ji,jj,jk) + zpronewd(ji,jj,jk) )
                  !
-                 zfeup = texcretn * zprofen(ji,jj,jk) + texcretd * zprofed(ji,jj,jk)
-                 tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zfeup
+                 tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - ( texcretn * zprofen(ji,jj,jk)    &
+                 &                   + texcretd * zprofed(ji,jj,jk) )
+                 consfe3(ji,jj,jk) = ( texcretn * zprofen(ji,jj,jk) + texcretd * zprofed(ji,jj,jk) )     &
+                 &                    * 75.0 / ( rtrn + ( plig(ji,jj,jk) + 75.0 * (1.0 - plig(ji,jj,jk) ) ) * trb(ji,jj,jk,jpfer) ) / rfact2
                  tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) - zprmaxd(ji,jj,jk) * zysopt(ji,jj,jk)   &
                  &                   * rfact2 * trb(ji,jj,jk,jpdia)
@@ -401 +382 @@
      ! when ln_ligand is set to .true. in the namelist. Ligand uptake is small 
      ! and based on the FeL model by Morel et al. (2008) and on the study of
-     ! Shaked and Lis (2012)
+     ! Shaked et al. (2020)
      ! -------------------------------------------------------------------------
      IF( ln_ligand ) THEN
@@ -412 +393 @@
                     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) / ( rtrn + plig(ji,jj,jk) + 2.E3 * (1.0 - plig(ji,jj,jk) ) ) * lthet
+                    &       - zfeup * plig(ji,jj,jk) / ( rtrn + plig(ji,jj,jk) + 75.0 * (1.0 - plig(ji,jj,jk) ) ) * lthet
                     zpligprod1(ji,jj,jk) = zdocprod * ldocp
                     zpligprod2(ji,jj,jk) = zfeup * plig(ji,jj,jk) / ( rtrn + plig(ji,jj,jk) &
-                    &                      + 2.E3 * (1.0 - plig(ji,jj,jk) ) ) * lthet
+                    &                      + 75.0 * (1.0 - plig(ji,jj,jk) ) ) * lthet
                  ENDIF
               END DO

NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zrem.F90

@@ -35 +35 @@
    REAL(wp), PUBLIC ::   xremikn    !: remineralisation rate of DON (p5z) 
    REAL(wp), PUBLIC ::   xremikp    !: remineralisation rate of DOP (p5z) 
-   REAL(wp), PUBLIC ::   xremik     !: remineralisation rate of DOC (p4z) 
    REAL(wp), PUBLIC ::   nitrif     !: NH4 nitrification rate 
    REAL(wp), PUBLIC ::   xsirem     !: remineralisation rate of biogenic silica
@@ -71 +70 @@
       REAL(wp) ::   zremik, zremikc, zremikn, zremikp, zsiremin, zfact 
       REAL(wp) ::   zsatur, zsatur2, znusil, znusil2, zdep, zdepmin, zfactdep
-      REAL(wp) ::   zbactfer, zolimit, zonitr, zrfact2
-      REAL(wp) ::   zammonic, zoxyremc, zoxyremn, zoxyremp
-      REAL(wp) ::   zosil, ztem, zdenitnh4, zolimic, zolimin, zolimip, zdenitrn, zdenitrp
+      REAL(wp) ::   zbactfer, zonitr
+      REAL(wp) ::   zammonic, zoxyremc, zosil, ztem, zdenitnh4, zolimic
       CHARACTER (len=25) :: charout
-      REAL(wp), DIMENSION(jpi,jpj    ) :: ztempbac
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdepbac, zolimi, zdepprod, zfacsi, zfacsib, zdepeff, zfebact
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdepbac, zolimi, zfacsi, zfacsib, zdepeff, zfebact
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d
       !!---------------------------------------------------------------------
@@ -83 +80 @@
       !
       ! Initialisation of arrays
-      zdepprod(:,:,:) = 1._wp
       zdepeff (:,:,:) = 0.3_wp
-      ztempbac(:,:)   = 0._wp
       zfacsib(:,:,:)  = xsilab / ( 1.0 - xsilab )
       zfebact(:,:,:)  = 0._wp
@@ -99 +94 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zdep = MAX( hmld(ji,jj), heup(ji,jj) )
-               IF( gdept_n(ji,jj,jk) < zdep ) THEN
-                  zdepbac(ji,jj,jk) = MIN( 0.7 * ( trb(ji,jj,jk,jpzoo) + 2.* trb(ji,jj,jk,jpmes) ), 4.e-6 )
-                  ztempbac(ji,jj)   = zdepbac(ji,jj,jk)
-               ELSE
+               zdep = MAX( hmld(ji,jj), heup_01(ji,jj) )
+               zdepbac(ji,jj,jk) = 0.6 * ( MAX(0.0, trb(ji,jj,jk,jpzoo) + trb(ji,jj,jk,jpmes) ) * 1.0E6 )**0.6 * 1.E-6
+               IF( gdept_n(ji,jj,jk) >= zdep ) THEN
                   zdepmin = MIN( 1., zdep / gdept_n(ji,jj,jk) )
-                  zdepbac (ji,jj,jk) = zdepmin**0.683 * ztempbac(ji,jj)
-                  zdepprod(ji,jj,jk) = zdepmin**0.273
                   zdepeff (ji,jj,jk) = zdepeff(ji,jj,jk) * zdepmin**0.3
                ENDIF
@@ -113 +104 @@
       END DO
 
-      IF( ln_p4z ) THEN ! Standard PISCES code
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ! DOC ammonification. Depends on a limitation term of the bacterial activity
-                  ! and on the implicit bacteria concentration
-                  ! --------------------------------------------------------------------------
-                  zremik = xremik * xstep / 1.e-6 * xlimbac(ji,jj,jk) * zdepbac(ji,jj,jk) 
-                  zremik = MAX( zremik, 2.74e-4 * xstep )
-
-                  ! Ammonification in oxic waters with oxygen consumption
-                  ! -----------------------------------------------------
-                  zolimit = zremik * ( 1.- nitrfac(ji,jj,jk) ) * trb(ji,jj,jk,jpdoc) 
-                  zolimi(ji,jj,jk) = MIN( ( trb(ji,jj,jk,jpoxy) - rtrn ) / o2ut, zolimit ) 
-
-                  ! Ammonification in suboxic waters with denitrification
-                  ! -----------------------------------------------------
-                  zammonic = zremik * nitrfac(ji,jj,jk) * trb(ji,jj,jk,jpdoc)
-                  denitr(ji,jj,jk)  = zammonic * ( 1. - nitrfac2(ji,jj,jk) )
-                  denitr(ji,jj,jk)  = MIN( ( trb(ji,jj,jk,jpno3) - rtrn ) / rdenit, denitr(ji,jj,jk) )
-
-                  ! Ammonification in waters depleted in O2 and NO3 based on 
-                  ! other redox processes
-                  ! --------------------------------------------------------
-                  zoxyremc          = zammonic - denitr(ji,jj,jk)
-                  !
-                  zolimi (ji,jj,jk) = MAX( 0.e0, zolimi (ji,jj,jk) )
-                  denitr (ji,jj,jk) = MAX( 0.e0, denitr (ji,jj,jk) )
-                  zoxyremc          = MAX( 0.e0, zoxyremc )
-                  ! Update of the TRA arrays
-                  tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) + zolimi (ji,jj,jk) + denitr(ji,jj,jk) + zoxyremc
-                  tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zolimi (ji,jj,jk) + denitr(ji,jj,jk) + zoxyremc
-                  tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) - denitr (ji,jj,jk) * rdenit
-                  tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - zolimi (ji,jj,jk) - denitr(ji,jj,jk) - zoxyremc
-                  tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) - zolimi (ji,jj,jk) * o2ut
-                  tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) + zolimi (ji,jj,jk) + denitr(ji,jj,jk) + zoxyremc
-                  tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * ( zolimi(ji,jj,jk) + zoxyremc    &
-                  &                     + ( rdenit + 1.) * denitr(ji,jj,jk) )
-               END DO
-            END DO
-         END DO
-      ELSE ! PISCES-QUOTA
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ! DOC ammonification. Depends on a limitation term of the bacterial activity
-                  ! and on the implicit bacteria concentration
-                  ! ---------------------------------------------------------------
-                  zremik = xstep / 1.e-6 * MAX(0.01, xlimbac(ji,jj,jk)) * zdepbac(ji,jj,jk) 
-                  zremik = MAX( zremik, 2.74e-4 * xstep / xremikc )
-
-                  zremikc = xremikc * zremik
-                  zremikn = xremikn / xremikc
-                  zremikp = xremikp / xremikc
-
-                  ! Ammonification in oxic waters with oxygen consumption
-                  ! -----------------------------------------------------
-                  zolimit = zremikc * ( 1.- nitrfac(ji,jj,jk) ) * trb(ji,jj,jk,jpdoc) 
-                  zolimic = MAX( 0.e0, MIN( ( trb(ji,jj,jk,jpoxy) - rtrn ) / o2ut, zolimit ) ) 
-                  zolimi(ji,jj,jk) = zolimic
-                  zolimin = zremikn * zolimic * trb(ji,jj,jk,jpdon) / ( trb(ji,jj,jk,jpdoc) + rtrn )
-                  zolimip = zremikp * zolimic * trb(ji,jj,jk,jpdop) / ( trb(ji,jj,jk,jpdoc) + rtrn ) 
-
-                  ! Ammonification in suboxic waters with denitrification
-                  ! -------------------------------------------------------
-                  zammonic = zremikc * nitrfac(ji,jj,jk) * trb(ji,jj,jk,jpdoc)
-                  denitr(ji,jj,jk)  = zammonic * ( 1. - nitrfac2(ji,jj,jk) )
-                  denitr(ji,jj,jk)  = MAX(0., MIN(  ( trb(ji,jj,jk,jpno3) - rtrn ) / rdenit, denitr(ji,jj,jk) ) )
-                  zoxyremc          = MAX(0., zammonic - denitr(ji,jj,jk))
-                  zdenitrn  = zremikn * denitr(ji,jj,jk) * trb(ji,jj,jk,jpdon) / ( trb(ji,jj,jk,jpdoc) + rtrn )
-                  zdenitrp  = zremikp * denitr(ji,jj,jk) * trb(ji,jj,jk,jpdop) / ( trb(ji,jj,jk,jpdoc) + rtrn )
-                  zoxyremn  = zremikn * zoxyremc * trb(ji,jj,jk,jpdon) / ( trb(ji,jj,jk,jpdoc) + rtrn )
-                  zoxyremp  = zremikp * zoxyremc * trb(ji,jj,jk,jpdop) / ( trb(ji,jj,jk,jpdoc) + rtrn )
-                  ! Update of the TRA arrays
-                  tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) + zolimip + zdenitrp + zoxyremp
-                  tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zolimin + zdenitrn + zoxyremn
-                  tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) - denitr(ji,jj,jk) * rdenit
-                  tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - zolimic - denitr(ji,jj,jk) - zoxyremc
-                  tra(ji,jj,jk,jpdon) = tra(ji,jj,jk,jpdon) - zolimin - zdenitrn - zoxyremn
-                  tra(ji,jj,jk,jpdop) = tra(ji,jj,jk,jpdop) - zolimip - zdenitrp - zoxyremp
-                  tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) - zolimic * o2ut
-                  tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) + zolimic + denitr(ji,jj,jk) + zoxyremc
-                  tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * ( zolimin + zoxyremn + ( rdenit + 1.) * zdenitrn )
-               END DO
-            END DO
-         END DO
-         !
-      ENDIF
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ! DOC ammonification. Depends on a limitation term of the bacterial activity
+               ! and on the implicit bacteria concentration
+               ! --------------------------------------------------------------------------
+               zremik = xstep / 1.e-6 * xlimbac(ji,jj,jk) * zdepbac(ji,jj,jk) 
+               zremik = MAX( zremik, 2.74e-4 * xstep / xremikc )
+               zremikc = xremikc * zremik
+
+               ! Ammonification in oxic waters with oxygen consumption
+               ! -----------------------------------------------------
+               zolimic = zremikc * ( 1.- nitrfac(ji,jj,jk) ) * trb(ji,jj,jk,jpdoc)
+               zolimic = MAX( 0.e0, MIN( ( trb(ji,jj,jk,jpoxy) - rtrn ) / o2ut, zolimic ) )
+               zolimi(ji,jj,jk) = zolimic
+
+               ! Ammonification in suboxic waters with denitrification
+               ! -----------------------------------------------------
+               zammonic = zremikc * nitrfac(ji,jj,jk) * trb(ji,jj,jk,jpdoc)
+               denitr(ji,jj,jk)  = zammonic * ( 1. - nitrfac2(ji,jj,jk) )
+               denitr(ji,jj,jk)  = MAX(0., MIN(  ( trb(ji,jj,jk,jpno3) - rtrn ) / rdenit, denitr(ji,jj,jk) ) )
+
+               ! Ammonification in waters depleted in O2 and NO3 based on 
+               ! other redox processes
+               ! --------------------------------------------------------
+               zoxyremc          = MAX(0., zammonic - denitr(ji,jj,jk) )
+
+               ! Update of the TRA arrays
+               tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) - denitr(ji,jj,jk) * rdenit
+               tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - ( zolimic + denitr(ji,jj,jk) + zoxyremc )
+               tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) - zolimic * o2ut
+               tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) + zolimic + denitr(ji,jj,jk) + zoxyremc
+               IF ( ln_p4z ) THEN ! PISCES-std
+                  tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) + zolimic + denitr(ji,jj,jk) + zoxyremc
+                  tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zolimic + denitr(ji,jj,jk) + zoxyremc
+                  tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * ( zolimic + zoxyremc + ( rdenit + 1.) * denitr(ji,jj,jk) )
+               ELSE  ! PISCES-QUOTA (p5z)
+                  zremikn = xremikn / xremikc * trb(ji,jj,jk,jpdon) / ( trb(ji,jj,jk,jpdoc) + rtrn )
+                  zremikp = xremikp / xremikc * trb(ji,jj,jk,jpdop) / ( trb(ji,jj,jk,jpdoc) + rtrn )
+                  tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) + zremikp * ( zolimic + denitr(ji,jj,jk) + zoxyremc )
+                  tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zremikn * ( zolimic + denitr(ji,jj,jk) + zoxyremc )
+                  tra(ji,jj,jk,jpdon) = tra(ji,jj,jk,jpdon) - zremikn * ( zolimic + denitr(ji,jj,jk) + zoxyremc )
+                  tra(ji,jj,jk,jpdop) = tra(ji,jj,jk,jpdop) - zremikp * ( zolimic + denitr(ji,jj,jk) + zoxyremc )
+                  tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * zremikn * ( zolimic + zoxyremc + ( rdenit + 1.) * denitr(ji,jj,jk) )
+               ENDIF
+            END DO
+         END DO
+      END DO
 
       DO jk = 1, jpkm1
@@ -211 +162 @@
                &         / ( 1.+ emoy(ji,jj,jk) ) * ( 1. + fr_i(ji,jj) * emoy(ji,jj,jk) ) 
                zdenitnh4 = nitrif * xstep * trb(ji,jj,jk,jpnh4) * nitrfac(ji,jj,jk)
-               zdenitnh4 = MIN(  ( trb(ji,jj,jk,jpno3) - rtrn ) / rdenita, zdenitnh4 ) 
+               zdenitnh4 = MAX(0., MIN(  ( trb(ji,jj,jk,jpno3) - rtrn ) / rdenita, zdenitnh4 ) ) 
                ! Update of the tracers trends
                ! ----------------------------
@@ -217 +168 @@
                tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) + zonitr - rdenita * zdenitnh4
                tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) - o2nit * zonitr
-               tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) - 2 * rno3 * zonitr + rno3 * ( rdenita - 1. ) * zdenitnh4
+               tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) - 2. * rno3 * zonitr + rno3 * ( rdenita - 1. ) * zdenitnh4
             END DO
          END DO
@@ -236 +187 @@
                ! studies (especially at Papa) have shown this uptake to be significant
                ! ---------------------------------------------------------------------
-               zbactfer = feratb *  rfact2 * 0.6_wp / rday * tgfunc(ji,jj,jk) * xlimbacl(ji,jj,jk)     &
-                  &              * trb(ji,jj,jk,jpfer) / ( xkferb + trb(ji,jj,jk,jpfer) )    &
-                  &              * zdepeff(ji,jj,jk) * zdepbac(ji,jj,jk)
+               zbactfer = feratb * 0.6_wp * xstep * tgfunc(ji,jj,jk) * xlimbacl(ji,jj,jk) * trb(ji,jj,jk,jpfer)    &
+                  &       / ( xkferb + trb(ji,jj,jk,jpfer) ) * zdepeff(ji,jj,jk) * zdepbac(ji,jj,jk)
 
                ! Only the transfer of iron from its dissolved form to particles
                ! is treated here. The GGE of bacteria supposed to be equal to 
                ! 0.33. This is hard-coded. 
-               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
+               tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zbactfer*0.18
+               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zbactfer*0.15
+               tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zbactfer*0.03
+               zfebact(ji,jj,jk)   = zbactfer * 0.18
                blim(ji,jj,jk)      = xlimbacl(ji,jj,jk)  * zdepbac(ji,jj,jk) / 1.e-6
             END DO
@@ -304 +254 @@
 
       IF( knt == nrdttrc ) THEN
-          zrfact2 = 1.e3 * rfact2r
           ALLOCATE( zw3d(jpi,jpj,jpk) )
           zfact = 1.e+3 * rfact2r  !  conversion from mol/l/kt to  mol/m3/s
@@ -321 +270 @@
           ENDIF
           IF( iom_use( "FEBACT" ) )  THEN
-               zw3d(:,:,:) = zfebact(:,:,:) * 1E9 * tmask(:,:,:) * zrfact2   ! Bacterial iron consumption
+               zw3d(:,:,:) = zfebact(:,:,:) * 1E9 * tmask(:,:,:) * zfact   ! Bacterial iron consumption
                CALL iom_put( "FEBACT" , zw3d )
           ENDIF
@@ -345 +294 @@
       !!
       !!----------------------------------------------------------------------
-      NAMELIST/nampisrem/ xremik, nitrif, xsirem, xsiremlab, xsilab, feratb, xkferb, & 
+      NAMELIST/nampisrem/ nitrif, xsirem, xsiremlab, xsilab, feratb, xkferb, & 
          &                xremikc, xremikn, xremikp
       INTEGER :: ios                 ! Local integer output status for namelist read
@@ -367 +316 @@
          WRITE(numout,*) '   Namelist parameters for remineralization, nampisrem'
          IF( ln_p4z ) THEN
-            WRITE(numout,*) '      remineralization rate of DOC              xremik    =', xremik
+            WRITE(numout,*) '      remineralization rate of DOC              xremikc   =', xremikc
          ELSE
             WRITE(numout,*) '      remineralization rate of DOC              xremikc   =', xremikc

NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zsbc.F90

@@ -483 +483 @@
       !
       sedsilfrac = 0.03     ! percentage of silica loss in the sediments
-      sedcalfrac = 0.6      ! percentage of calcite loss in the sediments
+      sedcalfrac = 0.99     ! percentage of calcite loss in the sediments
       !
    END SUBROUTINE p4z_sbc_init

NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zsed.F90

@@ -66 +66 @@
       REAL(wp) ::  zsiloss, zcaloss, zws3, zws4, zwsc, zdep
       REAL(wp) ::  zwstpoc, zwstpon, zwstpop
-      REAL(wp) ::  ztrfer, ztrpo4s, ztrdp, zwdust, zmudia, ztemp
+      REAL(wp) ::  ztrpo4s, ztrdp, zwdust, zmudia, ztemp
       REAL(wp) ::  xdiano3, xdianh4
       !
@@ -73 +73 @@
       REAL(wp), DIMENSION(jpi,jpj    ) :: zwsbio3, zwsbio4
       REAL(wp), DIMENSION(jpi,jpj    ) :: zsedcal, zsedsi, zsedc
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zsoufer, zlight
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrpo4, ztrdop, zirondep, zpdep
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) :: zsidep, zironice
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zsoufer, zlight, ztrfer
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrpo4, ztrdop, zirondep, zpdep, zsidep
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) :: zironice
       !!---------------------------------------------------------------------
       !
@@ -137 +137 @@
       IF( ln_dust ) THEN
          !                                              
-         ALLOCATE( zsidep(jpi,jpj), zpdep(jpi,jpj,jpk), zirondep(jpi,jpj,jpk) )
+         ALLOCATE( zsidep(jpi,jpj,jpk), zpdep(jpi,jpj,jpk), zirondep(jpi,jpj,jpk) )
 
          ! Iron, P and Si deposition at the surface
@@ -154 +154 @@
          ! and the solubility are hard coded
          ! ----------------------------------------------------------------
-         zsidep(:,:)   = 8.8 * 0.075 * dust(:,:) * mfrac * rfact2 / e3t_n(:,:,1) / 28.1 
-         zpdep (:,:,1) = 0.1 * 0.021 * dust(:,:) * mfrac * rfact2 / e3t_n(:,:,1) / 31. / po4r 
+         ! Si crustal abundance is about 26.9% by mass and 7.5% is soluble
+         ! (see for instance, Tegen and Kohfeld, 2006).
+         zsidep(:,:,1) = 0.269 * 0.075 * dust(:,:) * rfact2 / e3t_n(:,:,1) / 28.1 
+         ! P Crustal abundance is about 0.1% by mass and about 10% of it soluble
+         ! (Paytan and McLaughlin, 2007).
+         zpdep (:,:,1) = 0.1 * 1.e-3 * dust(:,:) * rfact2 / e3t_n(:,:,1) / 31. / po4r 
 
          ! Iron solubilization of particles in the water column
@@ -168 +172 @@
          DO jk = 2, jpkm1
             zirondep(:,:,jk) = dust(:,:) * mfrac * zwdust * rfact2 * EXP( -gdept_n(:,:,jk) / (250. * wdust) )
-            zpdep   (:,:,jk) = zirondep(:,:,jk) * 0.38 / po4r
+            zpdep   (:,:,jk) = zirondep(:,:,jk) * 1.e-3 / mfrac * 55.85 / 31.
+            zsidep  (:,:,jk) = zirondep(:,:,jk) * 0.269 / mfrac * 55.85 / 28.1
          END DO
 
          ! Solubilization of particles in the water column (Si, P, Fe)
-         tra(:,:,1,jpsil) = tra(:,:,1,jpsil) + zsidep  (:,:)
          DO jk = 1, jpkm1
             tra(:,:,jk,jppo4) = tra(:,:,jk,jppo4) + zpdep   (:,:,jk)
             tra(:,:,jk,jpfer) = tra(:,:,jk,jpfer) + zirondep(:,:,jk) 
+            tra(:,:,jk,jpsil) = tra(:,:,jk,jpsil) + zsidep  (:,:,jk)
          ENDDO
          ! 
@@ -299 +304 @@
                  zflx = (  trb(ji,jj,ikt,jpgoc) * zwsbio4(ji,jj)   &
                    &     + trb(ji,jj,ikt,jppoc) * zwsbio3(ji,jj) ) * 1E6
-                 zbureff(ji,jj) = 0.013 + 0.53 * zflx**2 / ( 7.0 + zflx )**2
+                 zbureff(ji,jj) = 0.013 + 0.13 * zflx**2 / ( 7.0 + zflx )**2
               ENDIF
             END DO
@@ -306 +311 @@
       ENDIF
 
-      ! Fraction of dSi that is remineralized in the sediments. This is 
-      ! set so that the burial in sediments equals the total input of Si
-      ! by rivers and dust (sedsilfrac)
-      ! ----------------------------------------------------------------
+      ! Fraction of dSi that is dissolved in the sediments. This fraction is  
+      ! set to a constant value in p4zsbc
+      ! --------------------------------------------------------------------
       IF( .NOT.lk_sed )  zrivsil = 1._wp - sedsilfrac
 
@@ -331 +335 @@
          ! modeled. The amount of CaCO3 that dissolves in the sediments
          ! is computed using a metamodel constructed from Archer (1996)
-         ! ------------------------------------------------------------
+         ! A minimum set to sedcalfrac is preserved. This value is defined
+         ! in p4zsbc
+         ! ---------------------------------------------------------------
          DO jj = 1, jpj
             DO ji = 1, jpi
@@ -341 +347 @@
                tra(ji,jj,ikt,jpsil) = tra(ji,jj,ikt,jpsil) + zsiloss * zrivsil 
                !
-               zfactcal = MIN( excess(ji,jj,ikt), 0.2 )
-               zfactcal = MIN( 1., 1.3 * ( 0.2 - zfactcal ) / ( 0.4 - zfactcal ) )
+               zfactcal = MAX(-0.1, MIN( excess(ji,jj,ikt), 0.2 ) )
+               zfactcal = 0.3 + 0.7 * MIN( 1., (0.1 + zfactcal) / ( 0.5 - zfactcal ) )
                zrivalk  = sedcalfrac * zfactcal
                tra(ji,jj,ikt,jptal) =  tra(ji,jj,ikt,jptal) + zcaloss * zrivalk * 2.0
@@ -394 +400 @@
                zws4 = zwsbio4(ji,jj) * zdep
                zws3 = zwsbio3(ji,jj) * zdep
+               ! Fraction that is permanently buried in the sediments
                zrivno3 = 1. - zbureff(ji,jj)
                zwstpoc = trb(ji,jj,ikt,jpgoc) * zws4 + trb(ji,jj,ikt,jppoc) * zws3
+               ! Denitrification in the sediments
                zpdenit  = MIN( 0.5 * ( trb(ji,jj,ikt,jpno3) - rtrn ) / rdenit, zdenit2d(ji,jj) * zwstpoc * zrivno3 )
+               ! Fraction that is not denitrified
                z1pdenit = zwstpoc * zrivno3 - zpdenit
+               ! Oxic remineralization of organic matter in the sediments
                zolimit = MIN( ( trb(ji,jj,ikt,jpoxy) - rtrn ) / o2ut, z1pdenit * ( 1.- nitrfac(ji,jj,ikt) ) )
+               ! The fraction that cannot be denitrified nor oxidized by O2
+               ! is released back to the water column as DOC
                tra(ji,jj,ikt,jpdoc) = tra(ji,jj,ikt,jpdoc) + z1pdenit - zolimit
+               ! Update of the tracers concentrations
                tra(ji,jj,ikt,jppo4) = tra(ji,jj,ikt,jppo4) + zpdenit + zolimit
                tra(ji,jj,ikt,jpnh4) = tra(ji,jj,ikt,jpnh4) + zpdenit + zolimit
@@ -408 +421 @@
                sdenit(ji,jj) = rdenit * zpdenit * e3t_n(ji,jj,ikt)
                zsedc(ji,jj)   = (1. - zrivno3) * zwstpoc * e3t_n(ji,jj,ikt)
+               ! PISCES-QUOTA (p5z)
                IF( ln_p5z ) THEN
                   zwstpop              = trb(ji,jj,ikt,jpgop) * zws4 + trb(ji,jj,ikt,jppop) * zws3
@@ -421 +435 @@
       ! Nitrogen fixation process : light limitation of diazotrophy
       ! Small source of iron from particulate inorganic iron (photochemistry)
+      ! This is a purely adhoc param.
       !----------------------------------------------------------------------
       DO jk = 1, jpkm1
          zlight (:,:,jk) =  ( 1.- EXP( -etot_ndcy(:,:,jk) / diazolight ) ) * ( 1. - fr_i(:,:) ) 
-         zsoufer(:,:,jk) = zlight(:,:,jk) * 2E-11 / ( 2E-11 + biron(:,:,jk) )
+         zsoufer(:,:,jk) = zlight(:,:,jk) * 1E-10 / ( 1E-10 + biron(:,:,jk) )
       ENDDO
 
@@ -438 +453 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  ! Potential nitrogen fixation dependant on temperature and iron
+                  ! Potential nitrogen fixation dependant on temperature
                   ztemp = tsn(ji,jj,jk,jp_tem)
-                  zmudia = MAX( 0.,-0.001096*ztemp**2 + 0.057*ztemp -0.637 ) * 7.625
-                  xdianh4 = trb(ji,jj,jk,jpnh4) / ( concnnh4 + trb(ji,jj,jk,jpnh4) )
-                  xdiano3 = trb(ji,jj,jk,jpno3) / ( concnno3 + trb(ji,jj,jk,jpno3) ) * (1. - xdianh4)
-                  zlim = ( 1.- xdiano3 - xdianh4 )
-                  ! Nitrogen fixation is almost fully halted when the N 
-                  ! limitation term (xdiano3+xdianh4) is > 0.9
-                  IF( zlim <= 0.1 )   zlim = 0.01
+                  zmudia = MAX( 0.,-0.001096*ztemp**2 + 0.057*ztemp -0.637 ) / rno3
+                  ! Nitrogen fixation is inhibited when enough NO3 and/or NH4
+                  zlim = ( 1.- xnanonh4(ji,jj,jk) - xnanono3(ji,jj,jk) )
                   zfact = zlim * rfact2
-                  ztrfer = biron(ji,jj,jk) / ( concfediaz + biron(ji,jj,jk) )
-                  ztrpo4(ji,jj,jk) = trb(ji,jj,jk,jppo4) / ( 1E-6 + trb(ji,jj,jk,jppo4) )
+                  ! Nitrogen fixation limitation by PO4 and Fe
+                  ztrfer(ji,jj,jk) = biron(ji,jj,jk) / ( concfediaz + biron(ji,jj,jk) )
+                  ztrpo4(ji,jj,jk) = trb(ji,jj,jk,jppo4) / ( concnnh4 + trb(ji,jj,jk,jppo4) )
                   ztrdp = ztrpo4(ji,jj,jk)
-                  nitrpot(ji,jj,jk) =  zmudia * r1_rday * zfact * MIN( ztrfer, ztrdp ) * zlight(ji,jj,jk)
+                  nitrpot(ji,jj,jk) =  zmudia * r1_rday * zfact * MIN( ztrfer(ji,jj,jk), ztrdp ) * zlight(ji,jj,jk)
                END DO
             END DO
@@ -460 +472 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  ! Potential nitrogen fixation dependant on temperature and iron
+                  ! Potential nitrogen fixation dependant on temperature
                   ztemp = tsn(ji,jj,jk,jp_tem)
                   zmudia = MAX( 0.,-0.001096*ztemp**2 + 0.057*ztemp -0.637 ) * 7.625
+                  ! Nitrogen fixation is inhibited when enough NO3 and/or NH4
                   xdianh4 = trb(ji,jj,jk,jpnh4) / ( concnnh4 + trb(ji,jj,jk,jpnh4) )
                   xdiano3 = trb(ji,jj,jk,jpno3) / ( concnno3 + trb(ji,jj,jk,jpno3) ) * (1. - xdianh4)
                   zlim = ( 1.- xdiano3 - xdianh4 )
-
-                  ! Nitrogen fixation is almost fully halted when the N 
-                  ! limitation term (xdiano3+xdianh4) is > 0.9
-                  IF( zlim <= 0.1 )   zlim = 0.01
                   zfact = zlim * rfact2
-                  ztrfer = biron(ji,jj,jk) / ( concfediaz + biron(ji,jj,jk) )
+                  ! Nitrogen fixation limitation by PO4/DOP and Fe
+                  ztrfer(ji,jj,jk) = biron(ji,jj,jk) / ( concfediaz + biron(ji,jj,jk) )
                   ztrpo4(ji,jj,jk) = trb(ji,jj,jk,jppo4) / ( 1E-6 + trb(ji,jj,jk,jppo4) )
                   ztrdop(ji,jj,jk) = trb(ji,jj,jk,jpdop) / ( 1E-6 + trb(ji,jj,jk,jpdop) ) * (1. - ztrpo4(ji,jj,jk))
                   ztrdp = ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk)
-                  nitrpot(ji,jj,jk) =  zmudia * r1_rday * zfact * MIN( ztrfer, ztrdp ) * zlight(ji,jj,jk)
+                  nitrpot(ji,jj,jk) =  zmudia * r1_rday * zfact * MIN( ztrfer(ji,jj,jk), ztrdp ) * zlight(ji,jj,jk)
                END DO
             END DO
@@ -490 +500 @@
                   zfact = nitrpot(ji,jj,jk) * nitrfix
                   ! 1/3 of the diazotrophs growth is supposed to be excreted
-                  ! as NH4. 1/3 as DOC and the rest is routed POC and GOC as 
+                  ! as NH4. 1/3 as DOC and the rest is routed to POC/GOC as 
                   ! a result of mortality by predation. Completely adhoc param 
                   tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zfact / 3.0
@@ -499 +509 @@
                   tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zfact * 1.0 / 3.0 * 1.0 / 3.0
                   tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + ( o2ut + o2nit ) * zfact * 2.0 / 3.0 + o2nit * zfact / 3.0
-                  ! Fe/c of diazotrophs is assumed to be 30umol Fe/mol C
-                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - 30E-6 * zfact * 1.0 / 3.0
-                  tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + 30E-6 * zfact * 1.0 / 3.0 * 2.0 / 3.0
-                  tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + 30E-6 * zfact * 1.0 / 3.0 * 1.0 / 3.0
-                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + 0.002 * 4E-10 * zsoufer(ji,jj,jk) * rfact2 / rday
+                  ! Fe/c of diazotrophs is assumed to be 30umol Fe/mol C at max
+                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - 30E-6 * ztrfer(ji,jj,jk) * zfact * 1.0 / 3.0
+                  tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + 30E-6 * ztrfer(ji,jj,jk) * zfact * 1.0 / 3.0 * 2.0 / 3.0
+                  tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + 30E-6 * ztrfer(ji,jj,jk) * zfact * 1.0 / 3.0 * 1.0 / 3.0
+                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + 0.005 * 4E-10 * zsoufer(ji,jj,jk) * rfact2 / rday
+                  tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) + concdnh4 / ( concdnh4 + trb(ji,jj,jk,jppo4) ) &
+                  &                     * 0.001 * trb(ji,jj,jk,jpdoc) * xstep
               END DO
             END DO 
@@ -533 +545 @@
                   tra(ji,jj,jk,jpgop) = tra(ji,jj,jk,jpgop) + 16.0 / 46.0 * zfact * 1.0 / 3.0 * 1.0 /3.0
                   tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + ( o2ut + o2nit ) * zfact * 2.0 / 3.0 + o2nit * zfact / 3.0
-                  ! Fe/c of diazotrophs is assumed to be 30umol Fe/mol C
-                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - 30E-6 * zfact * 1.0 / 3.0 
-                  tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + 30E-6 * zfact * 1.0 / 3.0 * 2.0 / 3.0
-                  tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + 30E-6 * zfact * 1.0 / 3.0 * 1.0 / 3.0
+                  ! Fe/c of diazotrophs is assumed to be 30umol Fe/mol C at max
+                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - 30E-6 * ztrfer(ji,jj,jk) * zfact * 1.0 / 3.0 
+                  tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + 30E-6 * ztrfer(ji,jj,jk) * zfact * 1.0 / 3.0 * 2.0 / 3.0
+                  tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + 30E-6 * ztrfer(ji,jj,jk) * zfact * 1.0 / 3.0 * 1.0 / 3.0
                   tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + 0.002 * 4E-10 * zsoufer(ji,jj,jk) * rfact2 / rday
               END DO
@@ -558 +570 @@
             IF( iom_use("SedSi" ) )  CALL iom_put( "SedSi",  zsedsi (:,:) * zfact ) ! Permanent burial of bSi in sediments
             IF( iom_use("SedC" ) )   CALL iom_put( "SedC",   zsedc  (:,:) * zfact ) ! Permanent burial of OC in sediments
-            IF( iom_use("Sdenit" ) ) CALL iom_put( "Sdenit", sdenit (:,:) * zfact * rno3 ) ! Denitrification in the sediments
+!            IF( iom_use("Sdenit" ) ) CALL iom_put( "Sdenit", sdenit (:,:) * zfact * rno3 ) ! Denitrification in the sediments
+            IF( iom_use("Sdenit" ) ) CALL iom_put( "Sdenit", zdenit2d(:,:) ) ! Denitrification in the sediments
          ENDIF
       ENDIF

NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zsms.F90

@@ -64 +64 @@
       REAL(wp) ::  ztra
       CHARACTER (len=25) :: charout
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jp_pisces) :: ztrbbio
       !!---------------------------------------------------------------------
       !
@@ -102 +103 @@
          END DO
       ENDIF
+
+      DO jn = jp_pcs0, jp_pcs1              !   Store the tracer concentrations before entering PISCES
+         ztrbbio(:,:,:,jn) = trb(:,:,:,jn)
+      END DO
       !
       IF( ll_sbc ) CALL p4z_sbc( kt )   ! external sources of nutrients 
@@ -147 +152 @@
             tra(:,:,:,jn) = 0._wp
          END DO
-         ! Euler-forward temporal scheme
-         IF( ln_top_euler ) THEN
-            DO jn = jp_pcs0, jp_pcs1
-               trn(:,:,:,jn) = trb(:,:,:,jn)
-            END DO
-         ENDIF
+         !
       END DO
-
+      !
+#endif
+      !
+      ! If ln_sediment is set to .true. then the sediment module is called
+      IF( ln_sediment ) THEN 
+         !
+         CALL sed_model( kt )     !  Main program of Sediment model
+         !
+      ENDIF
+      !
+      !
+      DO jn = jp_pcs0, jp_pcs1
+         tra(:,:,:,jn) = ( trb(:,:,:,jn) - ztrbbio(:,:,:,jn) ) * rfactr
+         trb(:,:,:,jn) = ztrbbio(:,:,:,jn)
+         ztrbbio(:,:,:,jn) = 0._wp
+      END DO
       !
       IF( l_trdtrc ) THEN
@@ -161 +176 @@
          END DO
       END IF
-#endif
-      !
-      ! If ln_sediment is set to .true. then the sediment module is called
-      IF( ln_sediment ) THEN 
-         !
-         CALL sed_model( kt )     !  Main program of Sediment model
-         ! Eulor forward temporal scheme
-         IF( ln_top_euler ) THEN
-            DO jn = jp_pcs0, jp_pcs1
-               trn(:,:,:,jn) = trb(:,:,:,jn)
-            END DO
-         ENDIF
-         !
-      ENDIF
       !
       IF( lrst_trc )  CALL p4z_rst( kt, 'WRITE' )  !* Write PISCES informations in restart file 
@@ -199 +200 @@
       INTEGER :: ios                 ! Local integer output status for namelist read
       !!
-      NAMELIST/nampisbio/ nrdttrc, wsbio, xkmort, ferat3, wsbio2, wsbio2max, wsbio2scale,    &
-         &                   ldocp, ldocz, lthet, no3rat3, po4rat3
+      NAMELIST/nampisbio/ nrdttrc, wsbio, xkmort, feratz, feratm, wsbio2, wsbio2max,    &
+         &                wsbio2scale, ldocp, ldocz, lthet, no3rat3, po4rat3
          !
       NAMELIST/nampisdmp/ ln_pisdmp, nn_pisdmp
@@ -226 +227 @@
          WRITE(numout,*) '      half saturation constant for mortality    xkmort      =', xkmort 
          IF( ln_p5z ) THEN
-            WRITE(numout,*) '      N/C in zooplankton                        no3rat3     =', no3rat3
-            WRITE(numout,*) '      P/C in zooplankton                        po4rat3     =', po4rat3
-         ENDIF
-         WRITE(numout,*) '      Fe/C in zooplankton                       ferat3      =', ferat3
+            WRITE(numout,*) '      N/C in zooplankton                     no3rat3     =', no3rat3
+            WRITE(numout,*) '      P/C in zooplankton                     po4rat3     =', po4rat3
+         ENDIF
+         WRITE(numout,*) '      Fe/C in microzooplankton                  feratz      =', feratz
+         WRITE(numout,*) '      Fe/C in microzooplankton                  feratz      =', feratm
          WRITE(numout,*) '      Big particles sinking speed               wsbio2      =', wsbio2
          WRITE(numout,*) '      Big particles maximum sinking speed       wsbio2max   =', wsbio2max
@@ -313 +315 @@
          ENDIF
 
+         ! Read the Fe3 consumption term by phytoplankton
+         IF( iom_varid( numrtr, 'Consfe3', ldstop = .FALSE. ) > 0 ) THEN
+            CALL iom_get( numrtr, jpdom_autoglo, 'Consfe3' , consfe3(:,:,:)  )
+         ELSE
+            consfe3(:,:,:) = 0._wp
+         ENDIF
+
+
          ! Read the cumulative total flux. If not in the restart file, it is set to 0          
          IF( iom_varid( numrtr, 'tcflxcum', ldstop = .FALSE. ) > 0 ) THEN  ! cumulative total flux of carbon
@@ -326 +336 @@
             sized(:,:,:) = 1.
          ENDIF
+         sized(:,:,:) = MAX( 1.0, sized(:,:,:) )
          IF( iom_varid( numrtr, 'sizen', ldstop = .FALSE. ) > 0 ) THEN
             CALL iom_get( numrtr, jpdom_autoglo, 'sizen' , sizen(:,:,:)  )
@@ -331 +342 @@
             sizen(:,:,:) = 1.
          ENDIF
+         sizen(:,:,:) = MAX( 1.0, sizen(:,:,:) )
 
          ! PISCES-QUOTA specific part
@@ -338 +350 @@
             IF( iom_varid( numrtr, 'sizep', ldstop = .FALSE. ) > 0 ) THEN
                CALL iom_get( numrtr, jpdom_autoglo, 'sizep' , sizep(:,:,:)  )
+               sizep(:,:,:) = MAX( 1.0, sizep(:,:,:) )
             ELSE
                sizep(:,:,:) = 1.
@@ -353 +366 @@
          CALL iom_rstput( kt, nitrst, numrtw, 'Silicalim', xksi(:,:) )    ! Si 1/2 saturation constant
          CALL iom_rstput( kt, nitrst, numrtw, 'Silicamax', xksimax(:,:) ) ! Si max concentration
+         CALL iom_rstput( kt, nitrst, numrtw, 'Consfe3', consfe3(:,:,:) ) ! Si max concentration
          CALL iom_rstput( kt, nitrst, numrtw, 'tcflxcum', t_oce_co2_flx_cum ) ! Cumulative CO2 flux
          CALL iom_rstput( kt, nitrst, numrtw, 'sizen', sizen(:,:,:) )  ! Size of nanophytoplankton
@@ -377 +391 @@
       !
       REAL(wp) ::  alkmean = 2426.     ! mean value of alkalinity ( Glodap ; for Goyet 2391. )
-      REAL(wp) ::  po4mean = 2.165     ! mean value of phosphate
-      REAL(wp) ::  no3mean = 30.90     ! mean value of nitrate
-      REAL(wp) ::  silmean = 91.51     ! mean value of silicate
+      REAL(wp) ::  po4mean = 2.174     ! mean value of phosphate
+      REAL(wp) ::  no3mean = 31.00     ! mean value of nitrate
+      REAL(wp) ::  silmean = 90.33     ! mean value of silicate
       !
       REAL(wp) :: zarea, zalksumn, zpo4sumn, zno3sumn, zsilsumn
@@ -534 +548 @@
          zwork(:,:,:) =   trn(:,:,:,jpfer) + trn(:,:,:,jpnfe) + trn(:,:,:,jpdfe)   &
             &         +   trn(:,:,:,jpbfe) + trn(:,:,:,jpsfe)                      &
-            &         + ( trn(:,:,:,jpzoo) + trn(:,:,:,jpmes) )  * ferat3    
+            &         + trn(:,:,:,jpzoo) * feratz + trn(:,:,:,jpmes) * feratm
          !
          ferbudget = glob_sum( 'p4zsms', zwork(:,:,:) * cvol(:,:,:)  )  

NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zlim.F90

@@ -88 +88 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   fvpuptk    !: Maximum potential uptake rate of picophyto
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   fvduptk    !: Maximum potential uptake rate of diatoms
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xqfuncfecp !: 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimnpn, xlimnpp, xlimnpd
 
    ! Coefficient for iron limitation following Flynn and Hipkin (1999)
@@ -127 +129 @@
       REAL(wp) ::   zconc0p, zconc0pnh4, zconc0ppo4, zconcpfe, zconcnfe, zconcdfe
       REAL(wp) ::   fanano, fananop, fananof, fadiat, fadiatp, fadiatf
-      REAL(wp) ::   fapico, fapicop, fapicof
+      REAL(wp) ::   fapico, fapicop, fapicof, zlimpo4, zlimdop
       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
-      REAL(wp) ::   zlim1f, zsizetmp
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: xlimnpn, xlimnpp, xlimnpd
+      REAL(wp) ::   zqfemn, zqfemp, zqfemd, zbactno3, zbactnh4, zbiron
+      REAL(wp) ::   znutlimtot, zlimno3, zlimnh4, zlim1f, zsizetmp
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zrassn, zrassp, zrassd
       !!---------------------------------------------------------------------
       !
@@ -143 +145 @@
          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 Chl/C ratio of each phytoplankton group
                ! ----------------------------------------------------------
@@ -190 +183 @@
                ! From Talmy et al. (2014) and Maranon et al. (2013)
                ! -------------------------------------------------------
-               xqnnmin(ji,jj,jk) = qnnmin * sizen(ji,jj,jk)**(-0.3)
+               xqnnmin(ji,jj,jk) = qnnmin * sizen(ji,jj,jk)**(-0.36)
                xqnnmax(ji,jj,jk) = qnnmax
-               xqndmin(ji,jj,jk) = qndmin * sized(ji,jj,jk)**(-0.3)
+               xqndmin(ji,jj,jk) = qndmin * sized(ji,jj,jk)**(-0.36)
                xqndmax(ji,jj,jk) = qndmax
-               xqnpmin(ji,jj,jk) = qnpmin * sizep(ji,jj,jk)**(-0.48)
-               xqnpmax(ji,jj,jk) = qnpmax * sizep(ji,jj,jk)**(-0.21)
+               xqnpmin(ji,jj,jk) = qnpmin * sizep(ji,jj,jk)**(-0.36)
+               xqnpmax(ji,jj,jk) = qnpmax
 
                ! Computation of the optimal allocation parameters
@@ -201 +194 @@
                ! Smith et al.
                ! ---------------------------------------------------
-
+               zbiron = ( 75.0 * ( 1.0 - plig(ji,jj,jk) ) + plig(ji,jj,jk) ) * biron(ji,jj,jk)
                ! Nanophytoplankton
                znutlim = MAX( trb(ji,jj,jk,jpnh4) / zconc0nnh4,    &
@@ -208 +201 @@
                znutlim = trb(ji,jj,jk,jppo4) / zconc0npo4
                fananop = MAX(0.01, MIN(0.99, 1. / ( SQRT(znutlim) + 1.) ) )
-               znutlim = biron(ji,jj,jk) / zconcnfe
+               znutlim = zbiron / zconcnfe
                fananof = MAX(0.01, MIN(0.99, 1. / ( SQRT(znutlim) + 1.) ) )
 
@@ -217 +210 @@
                znutlim = trb(ji,jj,jk,jppo4) / zconc0ppo4
                fapicop = MAX(0.01, MIN(0.99, 1. / ( SQRT(znutlim) + 1.) ) )
-               znutlim = biron(ji,jj,jk) / zconcpfe
+               znutlim = zbiron / zconcpfe
                fapicof = MAX(0.01, MIN(0.99, 1. / ( SQRT(znutlim) + 1.) ) )
 
@@ -226 +219 @@
                znutlim = trb(ji,jj,jk,jppo4) / zconc0dpo4
                fadiatp = MAX(0.01, MIN(0.99, 1. / ( SQRT(znutlim) + 1.) ) )
-               znutlim = biron(ji,jj,jk) / zconcdfe
+               znutlim = zbiron / zconcdfe
                fadiatf = MAX(0.01, MIN(0.99, 1. / ( SQRT(znutlim) + 1.) ) )
 
@@ -232 +225 @@
                ! heterotrophic 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)
+               zlimnh4 = trb(ji,jj,jk,jpnh4) / ( concbno3 + trb(ji,jj,jk,jpnh4) )
+               zlimno3 = trb(ji,jj,jk,jpno3) / ( concbno3 + trb(ji,jj,jk,jpno3) )
+               znutlimtot = ( trb(ji,jj,jk,jpnh4) + trb(ji,jj,jk,jpno3) ) / ( concbno3 + trb(ji,jj,jk,jpnh4) + trb(ji,jj,jk,jpno3) )
+               zbactnh4 = znutlimtot * 5.0 * zlimnh4 / ( zlimno3 + 5.0 * zlimnh4 + rtrn )
+               zbactno3 = znutlimtot * zlimno3 / ( zlimno3 + 5.0 * zlimnh4 + rtrn )
                !
                zlim1    = zbactno3 + zbactnh4
@@ -250 +246 @@
                ! Limitation of N based nutrients uptake (NO3 and NH4)
                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))
+               zlimnh4 = trb(ji,jj,jk,jpnh4) / ( zconc0n + trb(ji,jj,jk,jpnh4) )
+               zlimno3 = trb(ji,jj,jk,jpno3) / ( zconc0n + trb(ji,jj,jk,jpno3) )
+               znutlimtot = (1. - fanano) * ( trb(ji,jj,jk,jpnh4) + trb(ji,jj,jk,jpno3) ) / ( zfalim * zconc0n + trb(ji,jj,jk,jpnh4) + trb(ji,jj,jk,jpno3) )
+               xnanonh4(ji,jj,jk) = znutlimtot * 5.0 * zlimnh4 / ( zlimno3 + 5.0 * zlimnh4 + rtrn )
+               xnanono3(ji,jj,jk) = znutlimtot * zlimno3 / ( zlimno3 + 5.0 * zlimnh4 + rtrn )
                !
                ! Limitation of P based nutrients (PO4 and DOP)
                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.
+               zlimpo4 = trb(ji,jj,jk,jppo4) / ( trb(ji,jj,jk,jppo4) + zconc0npo4 )
+               zlimdop = trb(ji,jj,jk,jpdop) / ( trb(ji,jj,jk,jpdop) + zconc0npo4 )
+               znutlimtot = (1. - fananop) * ( trb(ji,jj,jk,jppo4) + trb(ji,jj,jk,jpdop) ) / ( zfalim * zconc0npo4 + trb(ji,jj,jk,jppo4) + trb(ji,jj,jk,jpdop) )
+               xnanopo4(ji,jj,jk) = znutlimtot * 100.0 * zlimpo4 / ( zlimdop + 100.0 * zlimpo4 + rtrn )
+               xnanodop(ji,jj,jk) = znutlimtot * zlimdop / ( zlimdop + 100.0 * zlimpo4 + rtrn )
                !
                ! Limitation of Fe uptake
                zfalim = (1.-fananof) / fananof
-               xnanofer(ji,jj,jk) = (1. - fananof) * biron(ji,jj,jk) / ( biron(ji,jj,jk) + zfalim * zconcnfe )
+               xnanofer(ji,jj,jk) = (1. - fananof) * zbiron / ( zbiron + zfalim * zconcnfe )
                !
                ! The minimum iron quota depends on the size of PSU, respiration
@@ -270 +269 @@
                zratiof   = trb(ji,jj,jk,jpnfe) * z1_trnphy
                zqfemn = xcoef1 * znanochl + xcoef2 + xcoef3 * xnanono3(ji,jj,jk)
+               xqfuncfecn(ji,jj,jk) = zqfemn + qfnopt
                !
                zration = trb(ji,jj,jk,jpnph) * z1_trnphy
                zration = MIN(xqnnmax(ji,jj,jk), MAX( xqnnmin(ji,jj,jk), zration ))
-               zzpsiuptk = xqnnmin(ji,jj,jk) * rno3 / zpsiuptk**2
-               fvnuptk(ji,jj,jk) = 1. / zzpsiuptk * xqnnmin(ji,jj,jk) / (zration + rtrn)  &
+               fvnuptk(ji,jj,jk) = 2.5 * zpsiuptk * xqnnmin(ji,jj,jk) / (zration + rtrn)  &
                &                   * MAX(0., (1. - zratchl * znanochl / 12. ) )
                !
@@ -282 +281 @@
                ! The value of the optimal quota in the formulation below
                ! has been found by solving a non linear equation
-               zlim1f = max(0., ( 1.086 - xqnnmin(ji,jj,jk) )  &
+               zlim1f = max(0., ( 1.13 - xqnnmin(ji,jj,jk) )  &
                &          / (xqnnmax(ji,jj,jk) - xqnnmin(ji,jj,jk) ) ) * xqnnmax(ji,jj,jk)
                zlim3  = MAX( 0.,( zratiof - zqfemn ) / qfnopt )
@@ -297 +296 @@
                ! Limitation of N based nutrients uptake (NO3 and NH4) 
                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))
+               zlimnh4 = trb(ji,jj,jk,jpnh4) / ( zconc0p + trb(ji,jj,jk,jpnh4) )
+               zlimno3 = trb(ji,jj,jk,jpno3) / ( zconc0p + trb(ji,jj,jk,jpno3) )
+               znutlimtot = (1. - fapico) * ( trb(ji,jj,jk,jpnh4) + trb(ji,jj,jk,jpno3) )   &
+               &            / ( zfalim * zconc0p + trb(ji,jj,jk,jpnh4) + trb(ji,jj,jk,jpno3) )
+               xpiconh4(ji,jj,jk) = znutlimtot * 5.0 * zlimnh4 / ( zlimno3 + 5.0 * zlimnh4 + rtrn )
+               xpicono3(ji,jj,jk) = znutlimtot * zlimno3 / ( zlimno3 + 5.0 * zlimnh4 + rtrn )
                !
                ! Limitation of P based nutrients uptake (PO4 and DOP)
                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.
+               zlimpo4 = trb(ji,jj,jk,jppo4) / ( trb(ji,jj,jk,jppo4) + zconc0ppo4 )
+               zlimdop = trb(ji,jj,jk,jpdop) / ( trb(ji,jj,jk,jpdop) + zconc0ppo4 )
+               znutlimtot = (1. - fapicop) * ( trb(ji,jj,jk,jppo4) + trb(ji,jj,jk,jpdop) ) / ( zfalim * zconc0ppo4 + trb(ji,jj,jk,jppo4) + trb(ji,jj,jk,jpdop) )
+               xpicopo4(ji,jj,jk) = znutlimtot * 100.0 * zlimpo4 / ( zlimdop + 100.0 * zlimpo4 + rtrn )
+               xpicodop(ji,jj,jk) = znutlimtot * zlimdop / ( zlimdop + 100.0 * zlimpo4 + rtrn )
                !
                zfalim = (1.-fapicof) / fapicof
-               xpicofer(ji,jj,jk) = (1. - fapicof) * biron(ji,jj,jk) / ( biron(ji,jj,jk) + zfalim * zconcpfe )
+               xpicofer(ji,jj,jk) = (1. - fapicof) * zbiron / ( zbiron + zfalim * zconcpfe )
                !
                ! The minimum iron quota depends on the size of PSU, respiration
                ! and the reduction of nitrate following the parameterization 
                ! proposed by Flynn and Hipkin (1999)
-               zratiof   = trb(ji,jj,jk,jppfe) * z1_trnpic
+               zratiof = trb(ji,jj,jk,jppfe) * z1_trnpic
                zqfemp = xcoef1 * zpicochl + xcoef2 + xcoef3 * xpicono3(ji,jj,jk)
+               xqfuncfecp(ji,jj,jk) = zqfemp + qfpopt
                !
                zration   = trb(ji,jj,jk,jpnpi) * z1_trnpic
                zration = MIN(xqnpmax(ji,jj,jk), MAX( xqnpmin(ji,jj,jk), zration ))
-               zzpsiuptk = xqnpmin(ji,jj,jk) * rno3 / zpsiuptk**2
-               fvpuptk(ji,jj,jk) = 1. / zzpsiuptk * xqnpmin(ji,jj,jk) / (zration + rtrn)  &
+               fvpuptk(ji,jj,jk) = 2.5 * zpsiuptk * xqnpmin(ji,jj,jk) / (zration + rtrn)  &
                &                   * MAX(0., (1. - zratchl * zpicochl / 12. ) ) 
                !
@@ -328 +331 @@
                ! The value of the optimal quota in the formulation below
                ! has been found by solving a non linear equation
-               zlim1f   = max(0., (1.367 - xqnpmin(ji,jj,jk) )  &
+               zlim1f   = max(0., (1.29 - xqnpmin(ji,jj,jk) )  &
                &          / (xqnpmax(ji,jj,jk) - xqnpmin(ji,jj,jk) ) ) * xqnpmax(ji,jj,jk)
                zlim3    = MAX( 0.,( zratiof - zqfemp ) / qfpopt )
@@ -344 +347 @@
                ! Limitation of N based nutrients uptake (NO3 and NH4)
                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))
+               zlimnh4 = trb(ji,jj,jk,jpnh4) / ( zconc1d + trb(ji,jj,jk,jpnh4) )
+               zlimno3 = trb(ji,jj,jk,jpno3) / ( zconc1d + trb(ji,jj,jk,jpno3) )
+               znutlimtot = (1.0 - fadiat) * ( trb(ji,jj,jk,jpnh4) + trb(ji,jj,jk,jpno3) )    &
+               &            / ( zfalim * zconc1d + trb(ji,jj,jk,jpnh4) + trb(ji,jj,jk,jpno3) )
+               xdiatnh4(ji,jj,jk) = znutlimtot * 5.0 * zlimnh4 / ( zlimno3 + 5.0 * zlimnh4 + rtrn )
+               xdiatno3(ji,jj,jk) = znutlimtot * zlimno3 / ( zlimno3 + 5.0 * zlimnh4 + rtrn )
                !
                ! Limitation of P based nutrients uptake (PO4 and DOP)
                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.
+               zlimpo4 = trb(ji,jj,jk,jppo4) / ( trb(ji,jj,jk,jppo4) + zconc0dpo4 )
+               zlimdop = trb(ji,jj,jk,jpdop) / ( trb(ji,jj,jk,jpdop) + zconc0dpo4 )
+               znutlimtot = (1. - fadiatp) * ( trb(ji,jj,jk,jppo4) + trb(ji,jj,jk,jpdop) ) / ( zfalim * zconc0dpo4 + trb(ji,jj,jk,jppo4) + trb(ji,jj,jk,jpdop) )
+               xdiatpo4(ji,jj,jk) = znutlimtot * 100.0 * zlimpo4 / ( zlimdop + 100.0 * zlimpo4 + rtrn )
+               xdiatdop(ji,jj,jk) = znutlimtot * zlimdop / ( zlimdop + 100.0 * zlimpo4 + rtrn )
                !
                ! Limitation of Fe uptake
                zfalim = (1.-fadiatf) / fadiatf
-               xdiatfer(ji,jj,jk) = (1. - fadiatf) * biron(ji,jj,jk) / ( biron(ji,jj,jk) + zfalim * zconcdfe )
+               xdiatfer(ji,jj,jk) = (1. - fadiatf) * zbiron / ( zbiron + zfalim * zconcdfe )
                !
                ! The minimum iron quota depends on the size of PSU, respiration
@@ -364 +371 @@
                zratiof   = trb(ji,jj,jk,jpdfe) * z1_trndia
                zqfemd = xcoef1 * zdiatchl + xcoef2 + xcoef3 * xdiatno3(ji,jj,jk)
+               xqfuncfecd(ji,jj,jk) = zqfemd + qfdopt
                !
                zration   = trb(ji,jj,jk,jpndi) * z1_trndia
                zration   = MIN(xqndmax(ji,jj,jk), MAX( xqndmin(ji,jj,jk), zration ))
-               zzpsiuptk = xqndmin(ji,jj,jk) * rno3 / zpsiuptk**2
-               fvduptk(ji,jj,jk) = 1. / zzpsiuptk * xqndmin(ji,jj,jk) / (zration + rtrn)   &
+               fvduptk(ji,jj,jk) = 2.5 * zpsiuptk * xqndmin(ji,jj,jk) / (zration + rtrn)   &
                &                   * MAX(0., (1. - zratchl * zdiatchl / 12. ) ) 
                !
@@ -376 +383 @@
                ! The value of the optimal quota in the formulation below
                ! has been found by solving a non linear equation
-               zlim1f   = max(0., (1.077 - xqndmin(ji,jj,jk) )    &
+65             zlim1f   = max(0., (1.13 - xqndmin(ji,jj,jk) )    &
                &          / (xqndmax(ji,jj,jk) - xqndmin(ji,jj,jk) ) )   &
                &          * xqndmax(ji,jj,jk)
@@ -408 +415 @@
                ! ------------------------------
                zfuptk = 0.2 + 0.12 / ( 3.0 * sizen(ji,jj,jk) + rtrn )
-               zrpho  = 1.54 * trb(ji,jj,jk,jpnch) / ( trb(ji,jj,jk,jpnph) * rno3 * 14. + rtrn )
+               ! Computed from Inomura et al. (2020) using Pavlova Lutheri
+               zrpho  = 11.55 * trb(ji,jj,jk,jpnch) / ( trb(ji,jj,jk,jpphy) * 12. + rtrn )
                zrass = MAX(0.62/4., ( 1. - zrpho - zfuptk ) * xlimnpn(ji,jj,jk) )
-               xqpnmin(ji,jj,jk) = ( 0.0 + 0.0078 + 0.62/4. * 0.0783 * xqnnmin(ji,jj,jk) ) * 16.
-               xqpnmax(ji,jj,jk) = ( zrpho * 0.0128 + zrass * 0.0783 ) * 16.
-               xqpnmax(ji,jj,jk) = xqpnmax(ji,jj,jk) * trb(ji,jj,jk,jpnph) / ( trb(ji,jj,jk,jpphy) + rtrn )  &
-               &      + (0.033 + 0.0078 ) * 16.
+               zrassn(ji,jj,jk) = zrass
+               xqpnmin(ji,jj,jk) = ( 0.0 + 0.0078 + 0.62/4. * 0.0783 ) * 16.
+               xqpnmax(ji,jj,jk) = ( zrpho * 0.0089 + zrass * 0.0783 ) * 16.
+               xqpnmax(ji,jj,jk) = xqpnmax(ji,jj,jk) + (0.033 + 0.0078 ) * 16.
                xqpnmax(ji,jj,jk) = MIN( qpnmax, xqpnmax(ji,jj,jk) )
-
 
                ! Size estimation of picophytoplankton based on total biomass
@@ -425 +432 @@
                ! N/P ratio of picophytoplankton
                ! ------------------------------
-               zfuptk = 0.2 + 0.12 / ( 0.5 * sizep(ji,jj,jk) + rtrn )
-               zrpho = 1.54 * trb(ji,jj,jk,jppch) / ( trb(ji,jj,jk,jpnpi) * rno3 * 14. + rtrn )
+               zfuptk = 0.2 + 0.12 / ( 0.8 * sizep(ji,jj,jk) + rtrn )
+               ! Computed from Inomura et al. (2020) using a synechococcus
+               zrpho = 13.4 * trb(ji,jj,jk,jppch) / ( trb(ji,jj,jk,jppic) * 12. + rtrn )
                zrass = MAX(0.4/4., ( 1. - zrpho - zfuptk ) * xlimnpp(ji,jj,jk) )
-               xqppmin(ji,jj,jk) = ( (0.0 + 0.0078 ) + 0.4/4. * 0.0517 * xqnpmin(ji,jj,jk) ) * 16.
-               xqppmax(ji,jj,jk) = ( zrpho * 0.0128 + zrass * 0.0517 ) * 16.
-               xqppmax(ji,jj,jk) = xqppmax(ji,jj,jk) * trb(ji,jj,jk,jpnpi) / ( trb(ji,jj,jk,jppic) + rtrn ) &
-               &      +  (0.033 + 0.0078 ) * 16
+               zrassp(ji,jj,jk) = zrass
+               xqppmin(ji,jj,jk) = ( (0.0 + 0.0078 ) + 0.4/4. * 0.0517 ) * 16.
+               xqppmax(ji,jj,jk) = ( zrpho * 0.0076 + zrass * 0.0517 ) * 16.
+               xqppmax(ji,jj,jk) = xqppmax(ji,jj,jk) +  (0.033 + 0.0078 ) * 16
                xqppmax(ji,jj,jk) = MIN( qppmax, xqppmax(ji,jj,jk) )
 
@@ -443 +451 @@
                ! --------------------
                zfuptk = 0.2 + 0.12 / ( 5.0 * sized(ji,jj,jk) + rtrn )
-               zrpho = 1.54 * trb(ji,jj,jk,jpdch) / ( trb(ji,jj,jk,jpndi) * rno3 * 14. + rtrn )
+               ! Computed from Inomura et al. (2020) using a synechococcus
+               zrpho = 8.08 * trb(ji,jj,jk,jpdch) / ( trb(ji,jj,jk,jpndi) * 12. + rtrn )
                zrass = MAX(0.66/4., ( 1. - zrpho - zfuptk ) * xlimnpd(ji,jj,jk) )
-
-               xqpdmin(ji,jj,jk) = ( ( 0.0 + 0.0078 ) + 0.66/4. * 0.0783 *  xqndmin(ji,jj,jk) ) * 16.
-               xqpdmax(ji,jj,jk) = ( zrpho * 0.0128 + zrass * 0.0783 ) * 16.
-               xqpdmax(ji,jj,jk) = xqpdmax(ji,jj,jk) * trb(ji,jj,jk,jpndi) / ( trb(ji,jj,jk,jpdia) + rtrn ) &
-               &      + ( 0.0078 + 0.033 ) * 16.
+               zrassd(ji,jj,jk)=zrass
+               xqpdmin(ji,jj,jk) = ( ( 0.0 + 0.0078 ) + 0.66/4. * 0.0783 ) * 16.
+               xqpdmax(ji,jj,jk) = ( zrpho * 0.0135 + zrass * 0.0783 ) * 16.
+               xqpdmax(ji,jj,jk) = xqpdmax(ji,jj,jk) + ( 0.0078 + 0.033 ) * 16.
                xqpdmax(ji,jj,jk) = MIN(qpdmax, xqpdmax(ji,jj,jk) )
 
@@ -467 +475 @@
                &        / ( 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) +  5.E-11 ) 
-               ztem1  = MAX( 0., tsn(ji,jj,jk,jp_tem) )
+               zlim3  = trb(ji,jj,jk,jpfer) / ( trb(ji,jj,jk,jpfer) +  6.E-11 ) 
+               ztem1  = MAX( 0., tsn(ji,jj,jk,jp_tem) + 1.8 )
                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 )    &
-               &                   * ztem1 / ( 1. + ztem1 ) * MAX( 1., trb(ji,jj,jk,jpphy)*1E6 )   &
+               zetot1 = MAX( 0., etot_ndcy(ji,jj,jk) - 1.) / ( 4. + etot_ndcy(ji,jj,jk) ) * 30. / ( 30. + etot_ndcy(ji,jj,jk) ) 
+
+               xfracal(ji,jj,jk) = caco3r * xlimphy(ji,jj,jk)     &
+               &                   * ztem1 / ( 0.1 + ztem1 ) * MAX( 1., trb(ji,jj,jk,jpphy)*1E6 )   &
                   &                * ( 1. + EXP(-ztem2 * ztem2 / 25. ) )         &
                   &                * zetot1 * MIN( 1., 50. / ( hmld(ji,jj) + rtrn ) )
@@ -508 +516 @@
         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
+        IF( iom_use( "RASSN"   ) ) CALL iom_put( "RASSN"  , zrassn(:,:,:) * tmask(:,:,:) )  ! Iron limitation term
+        IF( iom_use( "RASSP"   ) ) CALL iom_put( "RASSP"  , zrassp(:,:,:) * tmask(:,:,:) )  ! Iron limitation term
+        IF( iom_use( "RASSD"   ) ) CALL iom_put( "RASSD"  , zrassd(:,:,:) * tmask(:,:,:) )  ! Iron limitation term
       ENDIF
       !
@@ -640 +651 @@
          &      fvnuptk (jpi,jpj,jpk), fvduptk (jpi,jpj,jpk),       &
          &      xlimphys(jpi,jpj,jpk), xlimdias(jpi,jpj,jpk),       &
-         &      xlimpics(jpi,jpj,jpk),                              &
+         &      xlimnpp (jpi,jpj,jpk), xlimnpn (jpi,jpj,jpk),       &
+         &      xlimnpd (jpi,jpj,jpk),                              &
+         &      xlimpics(jpi,jpj,jpk), xqfuncfecp(jpi,jpj,jpk),     &
          &      fvpuptk (jpi,jpj,jpk), xlimpic (jpi,jpj,jpk),    STAT=ierr(1) )
          !

NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zmeso.F90

@@ -88 +88 @@
       REAL(wp) :: zepsherf, zepshert, zepsherq, zepsherv, zrespirc, zrespirn, zrespirp, zbasresb, zbasresi
       REAL(wp) :: zgraztotc, zgraztotn, zgraztotp, zgraztotf, zbasresn, zbasresp, zbasresf
-      REAL(wp) :: zgratmp, zgradoct, zgradont, zgrareft, zgradopt
+      REAL(wp) :: zgradoct, zgradont, zgrareft, zgradopt
       REAL(wp) :: zprcaca, zmortz, zexcess
       REAL(wp) :: zbeta, zrespz, ztortz, zgrasratp, zgrasratn, zgrasratf
@@ -104 +104 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zgrarep, zgraren, zgrapon, zgrapop
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zgradoc, zgradon, zgradop
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zgramigrem, zgramigref, zgramigpoc, zgramigpof, zstrn
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zgramigrem, zgramigref, zgramigpoc, zgramigpof
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zgramigrep, zgramigren, zgramigpop, zgramigpon
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zgramigdoc, zgramigdop, zgramigdon
@@ -190 +190 @@
                ! to avoid starvation.
                ! ----------------------------------------------------------
-               zsigma = 1.0 - zdenom**2/(0.05**2+zdenom**2)
+               zsigma = 1.0 - zdenom**3/(0.1**3+zdenom**3)
                zsigma = xsigma2 + xsigma2del * zsigma
                ! Nanophytoplankton and diatoms are the only preys considered
@@ -267 +267 @@
 
                zgraztotc  = zgrazdc + zgrazz + zgraznc + zgrazm + zgrazpoc + zgrazffep + zgrazffeg
-               zgraztotf  = zgrazdf + zgraznf + ( zgrazz + zgrazm ) * ferat3 + zgrazpof &
+               zgraztotf  = zgrazdf + zgraznf + zgrazz * feratz + zgrazm * feratm + zgrazpof &
                &            + zgrazfffp + zgrazfffg
                zgraztotn  = zgrazdn + (zgrazm + zgrazz) * no3rat3 + zgraznn + zgrazpon  &
@@ -291 +291 @@
                ! Fulton, 2012)
                ! -----------------------------------------------------------------------------------
-               zepshert  = MIN( 1., zgrasratn/ no3rat3, zgrasratp/ po4rat3, zgrasratf / ferat3)
+               zepshert  = MIN( 1., zgrasratn/ no3rat3, zgrasratp/ po4rat3, zgrasratf / feratm)
                zbeta     = MAX(0., (epsher2 - epsher2min) )
                zepsherf  = epsher2min + zbeta / ( 1.0 + 0.04E6 * 12. * zfood * zbeta )
@@ -312 +312 @@
                zexcess  = ( zgrasratp/ po4rat3 - zepshert ) / ( 1.0 - zepshert + rtrn)
                zbasresp = zbasresi * zexcess * zgrasratp
-               zexcess  = ( zgrasratf/ ferat3 - zepshert ) / ( 1.0 - zepshert + rtrn)
+               zexcess  = ( zgrasratf/ feratm - zepshert ) / ( 1.0 - zepshert + rtrn)
                zbasresf = zbasresi * zexcess * zgrasratf
 
@@ -320 +320 @@
                zgradont = (1. - unass2n) * zgraztotn - zepsherv * no3rat3 * zgraztotc - zbasresn
                zgradopt = (1. - unass2p) * zgraztotp - zepsherv * po4rat3 * zgraztotc - zbasresp
-               zgrareft = (1. - unass2c) * zgraztotf - zepsherv * ferat3 * zgraztotc - zbasresf
-               ztmp1   = ( 1. - epsher2 - unass2c ) /( 1. - 0.8 * epsher2 ) * ztortz
+               zgrareft = (1. - unass2c) * zgraztotf - zepsherv * feratm * zgraztotc - zbasresf
+               ztmp1    = (1. - epsher2 - unass2c) /( 1. - epsher2 ) * ztortz
                zgradoc(ji,jj,jk) = (zgradoct + ztmp1) * ssigma2
                zgradon(ji,jj,jk) = (zgradont + no3rat3 * ztmp1) * ssigma2
                zgradop(ji,jj,jk) = (zgradopt + po4rat3 * ztmp1) * ssigma2
-               zgratmp = 0.2 * epsher2 /( 1. - 0.8 * epsher2 ) * ztortz
 
                !  Since only semilabile DOM is represented in PISCES
@@ -331 +330 @@
                !  as dissolved inorganic compounds (ssigma2)
                !  --------------------------------------------------
-               zgrarem(ji,jj,jk) = zgratmp + ( zgradoct + ztmp1 ) * (1.0 - ssigma2)
-               zgraren(ji,jj,jk) = no3rat3 * zgratmp + ( zgradont + no3rat3 * ztmp1 ) * (1.0 - ssigma2)
-               zgrarep(ji,jj,jk) = po4rat3 * zgratmp + ( zgradopt + po4rat3 * ztmp1 ) * (1.0 - ssigma2)
-               zgraref(ji,jj,jk) = zgrareft + ferat3 * ( ztmp1 + zgratmp )
+               zgrarem(ji,jj,jk) = ( zgradoct + ztmp1 ) * (1.0 - ssigma2)
+               zgraren(ji,jj,jk) = ( zgradont + no3rat3 * ztmp1 ) * (1.0 - ssigma2)
+               zgrarep(ji,jj,jk) = ( zgradopt + po4rat3 * ztmp1 ) * (1.0 - ssigma2)
+               zgraref(ji,jj,jk) = zgrareft + feratm * ztmp1
 
                !   Defecation as a result of non assimilated products
                !   --------------------------------------------------
-               zgrapoc(ji,jj,jk)  = zgraztotc * unass2c + unass2c / ( 1. - 0.8 * epsher2 ) * ztortz
-               zgrapon(ji,jj,jk)  = zgraztotn * unass2n + no3rat3 * unass2n / ( 1. - 0.8 * epsher2 ) * ztortz
-               zgrapop(ji,jj,jk)  = zgraztotp * unass2p + po4rat3 * unass2p / ( 1. - 0.8 * epsher2 ) * ztortz
-               zgrapof(ji,jj,jk)  = zgraztotf * unass2c + ferat3  * unass2c / ( 1. - 0.8 * epsher2 ) * ztortz
+               zgrapoc(ji,jj,jk)  = zgraztotc * unass2c + unass2c / ( 1. - epsher2 ) * ztortz
+               zgrapon(ji,jj,jk)  = zgraztotn * unass2n + no3rat3 * unass2n / ( 1. - epsher2 ) * ztortz
+               zgrapop(ji,jj,jk)  = zgraztotp * unass2p + po4rat3 * unass2p / ( 1. - epsher2 ) * ztortz
+               zgrapof(ji,jj,jk)  = zgraztotf * unass2c + feratm  * unass2c / ( 1. - epsher2 ) * ztortz
 
                !  Addition of respiration to the release of inorganic nutrients
@@ -348 +347 @@
                zgraren(ji,jj,jk) = zgraren(ji,jj,jk) + zbasresn + zrespirc * no3rat3
                zgrarep(ji,jj,jk) = zgrarep(ji,jj,jk) + zbasresp + zrespirc * po4rat3
-               zgraref(ji,jj,jk) = zgraref(ji,jj,jk) + zbasresf + zrespirc * ferat3
+               zgraref(ji,jj,jk) = zgraref(ji,jj,jk) + zbasresf + zrespirc * feratm
 
                !   Update the arrays TRA which contain the biological sources and
@@ -404 +403 @@
          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
          zgramigrep(:,:)  = 0.0   ;   zgramigren(:,:) = 0.0
@@ -411 +409 @@
          zgramigdoc(:,:)  = 0.0   ;   zgramigdon(:,:) = 0.0
          zgramigdop(:,:)  = 0.0   
-
-         ! compute the day length depending on latitude and the day
-         zrum = REAL( nday_year - 80, wp ) / REAL( nyear_len(1), wp )
-         zcodel = ASIN(  SIN( zrum * rpi * 2._wp ) * SIN( rad * 23.5_wp )  )
-
-         ! day length in hours
-         zstrn(:,:) = 0.
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zargu = TAN( zcodel ) * TAN( gphit(ji,jj) * rad )
-               zargu = MAX( -1., MIN(  1., zargu ) )
-               zstrn(ji,jj) = MAX( 0.0, 24. - 2. * ACOS( zargu ) / rad / 15. )
-               zstrn(ji,jj) = MIN(0.75, MAX( 0.25, zstrn(ji,jj) / 24.) )
-            END DO
-         END DO
-
 
         ! Compute the amount of materials that will go into vertical migration
@@ -435 +417 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  zmigreltime = (1. - zstrn(ji,jj))
+                  zmigreltime = (1. - strn(ji,jj))
                   IF ( gdept_n(ji,jj,jk) <= heup(ji,jj) ) THEN
                      zgramigrem(ji,jj) = zgramigrem(ji,jj) + xfracmig * zgrarem(ji,jj,jk) * (1. - zmigreltime )    &
@@ -441 +423 @@
                      zgramigrep(ji,jj) = zgramigrep(ji,jj) + xfracmig * zgrarep(ji,jj,jk) * (1. - zmigreltime )    &
                      &                   * e3t_n(ji,jj,jk) * tmask(ji,jj,jk)
-                     zgramigrep(ji,jj) = zgramigren(ji,jj) + xfracmig * zgrarep(ji,jj,jk) * (1. - zmigreltime )    &
+                     zgramigren(ji,jj) = zgramigren(ji,jj) + xfracmig * zgraren(ji,jj,jk) * (1. - zmigreltime )    &
                      &                   * e3t_n(ji,jj,jk) * tmask(ji,jj,jk)
                      zgramigref(ji,jj) = zgramigref(ji,jj) + xfracmig * zgraref(ji,jj,jk) * (1. - zmigreltime )   &
@@ -503 +485 @@
          DEALLOCATE( zgramigrep, zgramigren, zgramigpop, zgramigpon )
          DEALLOCATE( zgramigdoc, zgramigdon, zgramigdop )
-         DEALLOCATE( zstrn )
       ! End of the ln_dvm_meso part
       ENDIF

NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zmicro.F90

@@ -110 +110 @@
                ! Proportion of nano and diatoms that are within the size range
                ! accessible to microzooplankton. 
-               zproport  = min(1.0, exp(-1.1 * MAX(0., ( sized(ji,jj,jk) - 1.8 ))**0.8 ))
-               zproport2 = sizen(ji,jj,jk)**(-0.54)
+               zproport  = min(sized(ji,jj,jk),1.8)**(-0.48)*min(1.0, exp(-1.1 * MAX(0., ( sized(ji,jj,jk) - 1.8 ))**0.8 ))
+               zproport2 = sizen(ji,jj,jk)**(-0.48)
+               zproport2 = 1.0
                !  linear mortality of mesozooplankton
                !  A michaelis menten modulation term is used to avoid extinction of 
@@ -162 +163 @@
                ! to avoid starvation.
                ! ----------------------------------------------------------
-               zsigma = 1.0 - zdenom**2/(0.05**2+zdenom**2)
+               zsigma = 1.0 - zdenom**3/(0.1**3+zdenom**3)
                zsigma = xsigma + xsigmadel * zsigma
-               zdiffpn = exp( -ABS(log(0.5 * sizep(ji,jj,jk) / (3.0 * sizen(ji,jj,jk) + rtrn )) )**2 / zsigma**2 )
+               zdiffpn = exp( -ABS(log(0.7 * sizep(ji,jj,jk) / (3.0 * sizen(ji,jj,jk) + rtrn )) )**2 / zsigma**2 )
                zdiffdn = exp( -ABS(log(3.0 * sizen(ji,jj,jk) / (5.0 * sized(ji,jj,jk) + rtrn )) )**2 / zsigma**2)
-               zdiffdp = exp( -ABS(log(0.5 * sizep(ji,jj,jk) / (5.0 * sized(ji,jj,jk) + rtrn )) )**2 / zsigma**2)
+               zdiffdp = exp( -ABS(log(0.7 * sizep(ji,jj,jk) / (5.0 * sized(ji,jj,jk) + rtrn )) )**2 / zsigma**2)
                ztmp1 = xprefn * zcompaph * ( zcompaph + zdiffdn * zcompadi + zdiffpn * zcompapi ) / ( 1.0 + zdiffdn + zdiffpn )
                ztmp2 = xprefp * zcompapi * ( zcompapi + zdiffpn * zcompaph + zdiffdp * zcompadi ) / ( 1.0 + zdiffpn + zdiffdp )
@@ -211 +212 @@
                zgraztotn = zgraznn + zgrazpn + zgrazpon + zgrazdn + zgrazz * no3rat3
                zgraztotp = zgraznp + zgrazpp + zgrazpop + zgrazdp + zgrazz * po4rat3
-               zgraztotf = zgraznf + zgrazpf + zgrazpof + zgrazdf + zgrazz * ferat3
+               zgraztotf = zgraznf + zgrazpf + zgrazpof + zgrazdf + zgrazz * feratz
                !
                ! Grazing by microzooplankton
@@ -230 +231 @@
                ! Fulton, 2012)
                ! -----------------------------------------------------------------------------------
-               zepshert  = MIN( 1., zgrasratn/ no3rat3, zgrasratp/ po4rat3, zgrasratf / ferat3)
+               zepshert  = MIN( 1., zgrasratn/ no3rat3, zgrasratp/ po4rat3, zgrasratf / feratz)
                zbeta     = MAX( 0., (epsher - epshermin) )
                ! Food density deprivation of GGE
@@ -255 +256 @@
                zexcess  = ( zgrasratp/ po4rat3 - zepshert ) / ( 1.0 - zepshert + rtrn)
                zbasresp = zbasresi * zexcess * zgrasratp
-               zexcess  = ( zgrasratf/ ferat3 - zepshert ) / ( 1.0 - zepshert + rtrn)
+               zexcess  = ( zgrasratf/ feratz - zepshert ) / ( 1.0 - zepshert + rtrn)
                zbasresf = zbasresi * zexcess * zgrasratf
 
@@ -263 +264 @@
                zgradont   = (1. - unassn) * zgraztotn - zepsherv * no3rat3 * zgraztotc - zbasresn
                zgradopt   = (1. - unassp) * zgraztotp - zepsherv * po4rat3 * zgraztotc - zbasresp
-               zgrareft   = (1. - unassc) * zgraztotf - zepsherv * ferat3 * zgraztotc - zbasresf
+               zgrareft   = (1. - unassc) * zgraztotf - zepsherv * feratz * zgraztotc - zbasresf
 
                ! Since only semilabile DOM is represented in PISCES
@@ -289 +290 @@
                zgraren = zgraren + zbasresn + zrespirc * no3rat3
                zgrarep = zgrarep + zbasresp + zrespirc * po4rat3
-               zgraref = zgraref + zbasresf + zrespirc * ferat3
+               zgraref = zgraref + zbasresf + zrespirc * feratz
 
                ! Update of the TRA arrays
@@ -330 +331 @@
                tra(ji,jj,jk,jppon) = tra(ji,jj,jk,jppon) + no3rat3 * ztortz + zgrapon - zgrazpon
                tra(ji,jj,jk,jppop) = tra(ji,jj,jk,jppop) + po4rat3 * ztortz + zgrapop - zgrazpop
-               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + ferat3 * ztortz  + zgrapof - zgrazpof
+               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + feratz * ztortz  + zgrapof - zgrazpof
                !
                ! Calcite production

NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zmort.F90

@@ -28 +28 @@
    REAL(wp), PUBLIC :: wchlp   !: Quadratic mortality rate of picophytoplankton
    REAL(wp), PUBLIC :: wchld   !: Quadratic mortality rate of diatoms
-   REAL(wp), PUBLIC :: wchldm  !: Maximum quadratic mortality rate of diatoms
    REAL(wp), PUBLIC :: mpratn  !: Linear mortality rate of nanophytoplankton
    REAL(wp), PUBLIC :: mpratp  !: Linear mortality rate of picophytoplankton
@@ -69 +68 @@
       !!---------------------------------------------------------------------
       INTEGER  :: ji, jj, jk
-      REAL(wp) :: zcompaph
+      REAL(wp) :: zcompaph, zlim1, zlim2
       REAL(wp) :: zfactfe, zfactch, zfactn, zfactp, zprcaca
       REAL(wp) :: ztortp , zrespp , zmortp
@@ -86 +85 @@
                ! blooms (Doney et al. 1996)
                ! -----------------------------------------------------
-               zrespp = wchln * 1.e6 * xstep * xdiss(ji,jj,jk) * zcompaph * trb(ji,jj,jk,jpphy)
+               zlim2   = xlimphy(ji,jj,jk) * xlimphy(ji,jj,jk)
+               zlim1   = 0.25 * ( 1. - zlim2 ) / ( 0.25 + zlim2 ) * trb(ji,jj,jk,jpphy)
+               zrespp = wchln * 1.e6 * xstep * zlim1 * xdiss(ji,jj,jk) * zcompaph
 
                ! Phytoplankton linear mortality
@@ -142 +143 @@
       !!---------------------------------------------------------------------
       INTEGER  :: ji, jj, jk
-      REAL(wp) :: zcompaph
+      REAL(wp) :: zcompaph, zlim1, zlim2
       REAL(wp) :: zfactfe, zfactch, zfactn, zfactp
       REAL(wp) :: ztortp , zrespp , zmortp 
@@ -158 +159 @@
                ! blooms (Doney et al. 1996)
                ! -----------------------------------------------------
-               zrespp = wchlp * 1.e6 * xstep * xdiss(ji,jj,jk) * zcompaph * trb(ji,jj,jk,jppic)
+               zlim2   = xlimpic(ji,jj,jk) * xlimpic(ji,jj,jk)
+               zlim1   = 0.25 * ( 1. - zlim2 ) / ( 0.25 + zlim2 ) * trb(ji,jj,jk,jppic)
+               zrespp = wchlp * 1.e6 * xstep * zlim1 * xdiss(ji,jj,jk) * zcompaph
 
                ! Phytoplankton linear mortality
@@ -229 +232 @@
                zlim2   = xlimdia(ji,jj,jk) * xlimdia(ji,jj,jk)
                zlim1   = 0.25 * ( 1. - zlim2 ) / ( 0.25 + zlim2 ) 
-               zrespp2 = 1.e6 * xstep * (  wchld + wchldm * zlim1 ) * xdiss(ji,jj,jk) * zcompadi * trb(ji,jj,jk,jpdia)
+               zrespp2 = 1.e6 * xstep * wchld * zlim1 * xdiss(ji,jj,jk) * zcompadi * trb(ji,jj,jk,jpdia)
 
                ! Phytoplankton linear mortality
@@ -291 +294 @@
       INTEGER :: ios   ! Local integer output status for namelist read
       !!
-      NAMELIST/namp5zmort/ wchln, wchlp, wchld, wchldm, mpratn, mpratp, mpratd
+      NAMELIST/namp5zmort/ wchln, wchlp, wchld, mpratn, mpratp, mpratd
       !!----------------------------------------------------------------------
 
@@ -310 +313 @@
          WRITE(numout,*) '    quadratic mortality of picophyto.         wchlp     =', wchlp
          WRITE(numout,*) '    quadratic mortality of diatoms            wchld     =', wchld
-         WRITE(numout,*) '    Additional quadratic mortality of diatoms wchldm    =', wchldm
          WRITE(numout,*) '    nanophyto. mortality rate                 mpratn    =', mpratn
          WRITE(numout,*) '    picophyto. mortality rate                 mpratp    =', mpratp

NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p5zprod.F90

@@ -72 +72 @@
       INTEGER  ::   ji, jj, jk
       REAL(wp) ::   zsilfac, znanotot, zpicotot, zdiattot, zconctemp, zconctemp2
-      REAL(wp) ::   zration, zratiop, zratiof, zmax, zsilim, ztn, zadap
-      REAL(wp) ::   zpronmax, zpropmax, zprofmax, zrat
+      REAL(wp) ::   zration, zratiop, zratiof, zmax, ztn, zadap
+      REAL(wp) ::   zpronmax, zpropmax, zprofmax, zratio
       REAL(wp) ::   zlim, zsilfac2, zsiborn, zprod, zprontot, zproptot, zprodtot
       REAL(wp) ::   zprnutmax, zdocprod, zprochln, zprochld, zprochlp
       REAL(wp) ::   zpislopen, zpislopep, zpisloped
       REAL(wp) ::   zrum, zcodel, zargu, zval, zfeup
-      REAL(wp) ::   zfact, zrfact2, zmaxsi, zratiosi, zsizetmp, zlimfac
+      REAL(wp) ::   zqfpmax, zqfnmax, zqfdmax
+      REAL(wp) ::   zfact, zrfact2, zmaxsi, zratiosi, zsizetmp, zlimfac, zsilim
       CHARACTER (len=25) :: charout
-      REAL(wp), DIMENSION(jpi,jpj    ) :: zmixnano, zmixpico, zmixdiat, zstrn
+      REAL(wp), DIMENSION(jpi,jpj    ) :: zmixnano, zmixpico, zmixdiat
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpislopeadn, zpislopeadp, zpislopeadd
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprnut, zprnutp, zprmaxp, zprmaxn, zprmaxd
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprnut, zprmaxp, zprmaxn, zprmaxd
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprbio, zprpic, zprdia, zysopt
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprchln, zprchlp, zprchld
@@ -111 +112 @@
       zysopt  (:,:,:) = 0._wp
       zrespn  (:,:,:) = 0._wp ; zrespp  (:,:,:) = 0._wp ; zrespd  (:,:,:) = 0._wp 
+      zmxl_fac(:,:,:) = 0._wp ; zmxl_chl(:,:,:) = 0._wp
+      consfe3 (:,:,:) = 0._wp
 
       ! Computation of the optimal production rates and nutrient uptake
       ! rates. Based on a Q10 description of the thermal dependency.
-      zprnut (:,:,:) = 0.6_wp * (1.0 + zpsino3 * qnnmax ) * r1_rday * tgfunc(:,:,:)
-      zprnutp(:,:,:) =  0.6_wp * (1. + zpsino3 * qnpmax ) * r1_rday * tgfunc3(:,:,:)
-      zprmaxn(:,:,:) = ( 0.6_wp * (1. + zpsino3 * qnnmax ) ) * r1_rday * tgfunc(:,:,:)
-      zprmaxd(:,:,:) = ( 0.6_wp * (1. + zpsino3 * qndmax ) ) * r1_rday * tgfunc(:,:,:)
-      zprmaxp(:,:,:) = ( 0.4_wp * (1. + zpsino3 * qnpmax ) ) * r1_rday * tgfunc3(:,:,:)
-
-      ! compute the day length depending on latitude and the day
-      ! Astronomical parameterization taken from HAMOCC3
-      zrum = REAL( nday_year - 80, wp ) / REAL( nyear_len(1), wp )
-      zcodel = ASIN(  SIN( zrum * rpi * 2._wp ) * SIN( rad * 23.5_wp )  )
-
-      ! day length in hours
-      zstrn(:,:) = 0.
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zargu = TAN( zcodel ) * TAN( gphit(ji,jj) * rad )
-            zargu = MAX( -1., MIN(  1., zargu ) )
-            zstrn(ji,jj) = MAX( 0.0, 24. - 2. * ACOS( zargu ) / rad / 15. )
-         END DO
-      END DO
+      zprnut (:,:,:) =  0.8_wp * r1_rday * tgfunc(:,:,:)
+      zprmaxn(:,:,:) =  0.8_wp * (1. + zpsino3 * qnnmax ) * r1_rday * tgfunc(:,:,:)
+      zprmaxd(:,:,:) =  0.8_wp * (1. + zpsino3 * qndmax ) * r1_rday * tgfunc(:,:,:)
+      zprmaxp(:,:,:) =  0.6_wp * (1. + zpsino3 * qnpmax ) * r1_rday * tgfunc(:,:,:)
 
       ! Impact of the day duration and light intermittency on phytoplankton growth
@@ -146 +133 @@
             DO ji = 1, jpi
                IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
-                  zval = MAX( 1., zstrn(ji,jj) )
+                  zval = MAX( 1., strn(ji,jj) )
                   IF( gdepw_n(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
                      zval = zval * MIN(1., heup_01(ji,jj) / ( hmld(ji,jj) + rtrn ))
@@ -161 +148 @@
       zprpic(:,:,:) = zprmaxp(:,:,:) * zmxl_fac(:,:,:)
 
-
       ! Maximum light intensity
-      zdaylen(:,:) = MAX(1., zstrn(:,:)) / 24.
+      zdaylen(:,:) = MAX(1., strn(:,:)) / 24.
 
       ! Computation of the P-I slope for nanos, picos and diatoms
@@ -192 +178 @@
                   ! Actual light levels are used here 
                   !  ---------------------------------------------
-                  zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen * enano(ji,jj,jk) / zmxl_chl(ji,jj,jk) )  )
-                  zprpic(ji,jj,jk) = zprpic(ji,jj,jk) * ( 1.- EXP( -zpislopep * epico(ji,jj,jk) / zmxl_chl(ji,jj,jk))  )
-                  zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediat(ji,jj,jk) / zmxl_chl(ji,jj,jk))  )
+                  zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen * enano(ji,jj,jk) / zmxl_fac(ji,jj,jk) )  )
+                  zprpic(ji,jj,jk) = zprpic(ji,jj,jk) * ( 1.- EXP( -zpislopep * epico(ji,jj,jk) / zmxl_fac(ji,jj,jk) )  )
+                  zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediat(ji,jj,jk) / zmxl_fac(ji,jj,jk) )  )
 
                   !  Computation of production function for Chlorophyll
@@ -227 +213 @@
                   ! -----------------------------------------------------------------------
                   zlim  = trb(ji,jj,jk,jpsil) / ( trb(ji,jj,jk,jpsil) + xksi1 )
-                  zsilim = MIN( zprdia(ji,jj,jk) / ( zprmaxd(ji,jj,jk) + rtrn ), xlimsi(ji,jj,jk) )
+                  zsilim = xlimdia(ji,jj,jk) * zprdia(ji,jj,jk) / ( zprmaxd(ji,jj,jk) + rtrn )
                   zsiborn = trb(ji,jj,jk,jpsil) * trb(ji,jj,jk,jpsil) * trb(ji,jj,jk,jpsil)
                   IF (gphit(ji,jj) < -30 ) THEN
@@ -240 +226 @@
                      zysopt(ji,jj,jk) = zlim * zsilfac2 * grosip * 1.0 * zmaxsi
                   ELSE
-                     zysopt(ji,jj,jk) = zlim * zsilfac2 * grosip * 1.0 * zsilim**0.75 * zmaxsi
+                     zysopt(ji,jj,jk) = zlim * zsilfac2 * grosip * 1.0 * zsilim**0.7 * zmaxsi
                   ENDIF
                ENDIF
@@ -257 +243 @@
                zprdia(ji,jj,jk)  = zprdia(ji,jj,jk) * ( 1. - fr_i(ji,jj) ) 
                zprnut(ji,jj,jk)  = zprnut(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
-               zprnutp(ji,jj,jk) = zprnutp(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
             END DO
          END DO
@@ -284 +269 @@
                   zlimfac = xlimphys(ji,jj,jk) * zprchln(ji,jj,jk) / ( zprmaxn(ji,jj,jk) + rtrn )
                   zsizetmp = 1.0 + 1.3 * ( xsizern - 1.0 ) * zlimfac**3/(0.3 + zlimfac**3)
-                  sizena(ji,jj,jk) = min(xsizern, max( sizena(ji,jj,jk), zsizetmp ) )
+                  sizena(ji,jj,jk) = MIN(xsizern, MAX( sizena(ji,jj,jk), zsizetmp ) )
                   ! Maximum potential uptake rate
                   zration = trb(ji,jj,jk,jpnph) / ( trb(ji,jj,jk,jpphy) + rtrn )
@@ -291 +276 @@
                   zprnutmax = zprnut(ji,jj,jk) * fvnuptk(ji,jj,jk) / rno3 * trb(ji,jj,jk,jpphy) * rfact2
                   ! Uptake of nitrogen
-                  zrat = 1.0 - MIN( 1., zration / (xqnnmax(ji,jj,jk) + rtrn) )
-                  zmax = MAX(0., MIN(1., zrat**2 / (0.05**2 + zrat**2) ) )
+                  zratio = 1.0 - MIN( 1., zration / (xqnnmax(ji,jj,jk) + rtrn) )
+                  zmax = MAX(0., MIN(1., zratio**2 / (0.05**2 + zratio**2) ) )
                   zpronmax = zprnutmax * zmax * MAX(0., MIN(1., ( zratiop - xqpnmin(ji,jj,jk) )   &
                   &          / ( xqpnmax(ji,jj,jk) - xqpnmin(ji,jj,jk) + rtrn ), xlimnfe(ji,jj,jk) ) )
+                  zpronmax = zpronmax * xqnnmin(ji,jj,jk) / qnnmin
                   zpronewn(ji,jj,jk) = zpronmax * xnanono3(ji,jj,jk)
                   zproregn(ji,jj,jk) = zpronmax * xnanonh4(ji,jj,jk)
                   ! Uptake of phosphorus and DOP
-                  zrat = 1.0 - MIN( 1., zratiop / (xqpnmax(ji,jj,jk) + rtrn) )
-                  zmax = MAX(0., MIN(1., zrat**2 / (0.05**2 + zrat**2) ) )
-                  zpropmax = zprnutmax * zmax * xlimnfe(ji,jj,jk) * 16. / 10.
+                  zratio = 1.0 - MIN( 1., zratiop / (xqpnmax(ji,jj,jk) + rtrn) )
+                  zmax = MAX(0., MIN(1., zratio**2 / (0.05**2 + zratio**2) ) )
+                  zpropmax = zprnutmax * zmax * xlimnfe(ji,jj,jk)
                   zpropo4n(ji,jj,jk) = zpropmax * xnanopo4(ji,jj,jk)
                   zprodopn(ji,jj,jk) = zpropmax * xnanodop(ji,jj,jk)
                   ! Uptake of iron
-                  zrat = 1.0 - MIN( 1., zratiof / qfnmax )
-                  zmax = MAX(0., MIN(1., zrat**2/ (0.05**2 + zrat**2) ) )
-                  zprofmax = zprnutmax * qfnmax * zmax 
+                  zqfnmax = xqfuncfecn(ji,jj,jk) + ( qfnmax - xqfuncfecn(ji,jj,jk) ) * xlimnpn(ji,jj,jk)
+                  zratio = 1.0 - MIN( 1., zratiof / zqfnmax )
+                  zmax = MAX(0., MIN(1., zratio**2/ (0.05**2 + zratio**2) ) )
+                  zprofmax = zprnutmax * zqfnmax * zmax 
                   zprofen(ji,jj,jk) = zprofmax * xnanofer(ji,jj,jk)    &
                   &          * (1. + 0.8 * xnanono3(ji,jj,jk) / ( rtrn  &
@@ -341 +328 @@
                   zratiop = trb(ji,jj,jk,jpppi) / ( trb(ji,jj,jk,jppic) + rtrn )
                   zratiof = trb(ji,jj,jk,jppfe) / ( trb(ji,jj,jk,jppic) + rtrn )
-                  zprnutmax = zprnutp(ji,jj,jk) * fvpuptk(ji,jj,jk) / rno3 * trb(ji,jj,jk,jppic) * rfact2
+                  zprnutmax = zprnut(ji,jj,jk) * fvpuptk(ji,jj,jk) / rno3 * trb(ji,jj,jk,jppic) * rfact2
                   ! Uptake of nitrogen
-                  zrat = 1.0 - MIN( 1., zration / (xqnpmax(ji,jj,jk) + rtrn) )
-                  zmax = MAX(0., MIN(1., zrat**2/ (0.05**2 + zrat**2) ) )
+                  zratio = 1.0 - MIN( 1., zration / (xqnpmax(ji,jj,jk) + rtrn) )
+                  zmax = MAX(0., MIN(1., zratio**2/ (0.05**2 + zratio**2) ) )
                   zpronmax = zprnutmax * zmax * MAX(0., MIN(1., ( zratiop - xqppmin(ji,jj,jk) )   &
                   &          / ( xqppmax(ji,jj,jk) - xqppmin(ji,jj,jk) + rtrn ), xlimpfe(ji,jj,jk) ) )
+                  zpronmax = zpronmax * xqnpmin(ji,jj,jk) / qnnmin
                   zpronewp(ji,jj,jk) = zpronmax * xpicono3(ji,jj,jk) 
                   zproregp(ji,jj,jk) = zpronmax * xpiconh4(ji,jj,jk)
                   ! Uptake of phosphorus
-                  zrat = 1.0 - MIN( 1., zratiop / (xqppmax(ji,jj,jk) + rtrn) )
-                  zmax = MAX(0., MIN(1., zrat**2 / (0.05**2 + zrat**2) ) )
-                  zpropmax = zprnutmax * zmax * xlimpfe(ji,jj,jk) * 16./10.
+                  zratio = 1.0 - MIN( 1., zratiop / (xqppmax(ji,jj,jk) + rtrn) )
+                  zmax = MAX(0., MIN(1., zratio**2 / (0.05**2 + zratio**2) ) )
+                  zpropmax = zprnutmax * zmax * xlimpfe(ji,jj,jk) 
                   zpropo4p(ji,jj,jk) = zpropmax * xpicopo4(ji,jj,jk)
                   zprodopp(ji,jj,jk) = zpropmax * xpicodop(ji,jj,jk)
                   ! Uptake of iron
-                  zrat = 1.0 - MIN( 1., zratiof / qfpmax )
-                  zmax = MAX(0., MIN(1., zrat**2 / (0.05**2 + zrat**2) ) )
-                  zprofmax = zprnutmax * qfpmax * zmax
+                  zqfpmax = xqfuncfecp(ji,jj,jk) + ( qfpmax - xqfuncfecp(ji,jj,jk) ) * xlimnpp(ji,jj,jk)
+                  zratio = 1.0 - MIN( 1., zratiof / zqfpmax )
+                  zmax = MAX(0., MIN(1., zratio**2 / (0.05**2 + zratio**2) ) )
+                  zprofmax = zprnutmax * zqfpmax * zmax
                   zprofep(ji,jj,jk) = zprofmax * xpicofer(ji,jj,jk)  &
                   &          * (1. + 0.8 * xpicono3(ji,jj,jk) / ( rtrn   &
@@ -395 +384 @@
                   zprnutmax = zprnut(ji,jj,jk) * fvduptk(ji,jj,jk) / rno3 * trb(ji,jj,jk,jpdia) * rfact2
                   ! Uptake of nitrogen
-                  zrat = 1.0 - MIN( 1., zration / (xqndmax(ji,jj,jk) + rtrn) )
-                  zmax = MAX(0., MIN(1., zrat**2 / (0.05**2 + zrat**2) ) )
+                  zratio = 1.0 - MIN( 1., zration / (xqndmax(ji,jj,jk) + rtrn) )
+                  zmax = MAX(0., MIN(1., zratio**2 / (0.05**2 + zratio**2) ) )
                   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) ) )
+                  zpronmax = zpronmax * xqndmin(ji,jj,jk) / qnnmin
                   zpronewd(ji,jj,jk) = zpronmax * xdiatno3(ji,jj,jk)
                   zproregd(ji,jj,jk) = zpronmax * xdiatnh4(ji,jj,jk)
                   ! Uptake of phosphorus
-                  zrat = 1.0 - MIN( 1., zratiop / (xqpdmax(ji,jj,jk) + rtrn) )
-                  zmax = MAX(0., MIN(1., zrat**2/ (0.05**2 + zrat**2) ) )
-                  zpropmax = zprnutmax * zmax * xlimdfe(ji,jj,jk) * 16./10.
+                  zratio = 1.0 - MIN( 1., zratiop / (xqpdmax(ji,jj,jk) + rtrn) )
+                  zmax = MAX(0., MIN(1., zratio**2/ (0.05**2 + zratio**2) ) )
+                  zpropmax = zprnutmax * zmax * xlimdfe(ji,jj,jk)
                   zpropo4d(ji,jj,jk) = zpropmax * xdiatpo4(ji,jj,jk)
                   zprodopd(ji,jj,jk) = zpropmax * xdiatdop(ji,jj,jk)
                   ! Uptake of iron
-                  zrat = 1.0 - MIN( 1., zratiof / qfdmax )
-                  zmax = MAX(0., MIN(1., zrat**2 / (0.05**2 + zrat**2) ) )
-                  zprofmax = zprnutmax * qfdmax * zmax
+                  zqfdmax = xqfuncfecd(ji,jj,jk) + ( qfdmax - xqfuncfecd(ji,jj,jk) ) * xlimnpd(ji,jj,jk)
+                  zratio = 1.0 - MIN( 1., zratiof / zqfdmax )
+                  zmax = MAX(0., MIN(1., zratio**2 / (0.05**2 + zratio**2) ) )
+                  zprofmax = zprnutmax * zqfdmax * zmax
                   zprofed(ji,jj,jk) = zprofmax * xdiatfer(ji,jj,jk)    &
                   &          * (1. + 0.8 * xdiatno3(ji,jj,jk) / ( rtrn   &
@@ -503 +494 @@
               zfeup = texcretn * zprofen(ji,jj,jk) + texcretd * zprofed(ji,jj,jk) + texcretp * zprofep(ji,jj,jk)
               tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zfeup
+              consfe3(ji,jj,jk)   = zfeup * 75.0 / ( rtrn + ( plig(ji,jj,jk) + 75.0 * (1.0 - plig(ji,jj,jk) ) )   &
+              &                     * trb(ji,jj,jk,jpfer) ) / rfact2
               tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) - zprmaxd(ji,jj,jk) * zysopt(ji,jj,jk) * rfact2 * trb(ji,jj,jk,jpdia)
               tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) - zprorcan(ji,jj,jk) - zprorcad(ji,jj,jk) - zprorcap(ji,jj,jk)  &
@@ -529 +522 @@
                  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) / ( rtrn + plig(ji,jj,jk) + 2.E3 * (1.0 - plig(ji,jj,jk) ) ) * lthet
+                 &       - zfeup * plig(ji,jj,jk) / ( rtrn + plig(ji,jj,jk) + 75.0 * (1.0 - plig(ji,jj,jk) ) ) * lthet
                  zpligprod1(ji,jj,jk) = zdocprod * ldocp
                  zpligprod2(ji,jj,jk) = zfeup * plig(ji,jj,jk) / ( rtrn + plig(ji,jj,jk) &
-                 &                      + 2.E3 * (1.0 - plig(ji,jj,jk) ) ) * lthet
+                 &                      + 75.0 * (1.0 - plig(ji,jj,jk) ) ) * lthet
               END DO
            END DO