Index: NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zsed.F90
===================================================================
--- NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zsed.F90	(revision 13233)
+++ NEMO/branches/2019/dev_r11708_aumont_PISCES_QUOTA/src/TOP/PISCES/P4Z/p4zsed.F90	(revision 14276)
@@ -66,5 +66,5 @@
       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,7 +73,7 @@
       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,5 +137,5 @@
       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,6 +154,10 @@
          ! 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,12 +172,13 @@
          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,5 +304,5 @@
                  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,8 +311,7 @@
       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,5 +335,7 @@
          ! 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,6 +347,6 @@
                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,10 +400,17 @@
                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,4 +421,5 @@
                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,8 +435,9 @@
       ! 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,18 +453,15 @@
             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,20 +472,18 @@
             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,5 +500,5 @@
                   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,9 +509,11 @@
                   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,8 +545,8 @@
                   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,5 +570,6 @@
             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
