Index: NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/P4Z/p4zsed.F90
===================================================================
--- NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/P4Z/p4zsed.F90	(revision 14384)
+++ NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/P4Z/p4zsed.F90	(revision 14385)
@@ -66,13 +66,12 @@
       REAL(wp) ::  zsiloss, zcaloss, zws3, zws4, zwsc, zdep
       REAL(wp) ::  zwstpoc, zwstpon, zwstpop
-      REAL(wp) ::  ztrfer, ztrpo4s, ztrdp, zwdust, zmudia, ztemp
+      REAL(wp) ::  ztrfer, ztrpo4, ztrdop, zmudia, ztemp
       REAL(wp) ::  xdiano3, xdianh4
+      REAL(wp) ::  zsoufer, zlight
       !
       CHARACTER (len=25) :: charout
-      REAL(wp), DIMENSION(jpi,jpj    ) :: zdenit2d, zbureff, zwork
+      REAL(wp), DIMENSION(jpi,jpj    ) :: zdenit2d, zbureff
       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
       !!---------------------------------------------------------------------
       !
@@ -81,26 +80,27 @@
 
       ! Allocate temporary workspace
-      ALLOCATE( ztrpo4(jpi,jpj,jpk) )
-      IF( ln_p5z )    ALLOCATE( ztrdop(jpi,jpj,jpk) )
-
       zdenit2d(:,:) = 0.e0
       zbureff (:,:) = 0.e0
-      zwork   (:,:) = 0.e0
       zsedsi  (:,:) = 0.e0
       zsedcal (:,:) = 0.e0
       zsedc   (:,:) = 0.e0
 
+      ! OA: Warning, the following part is necessary to avoid CFL problems 
+      ! above the sediments. Vertical sinking speed is limited using the 
+      ! typical CFL criterion
+      ! --------------------------------------------------------------------
+      DO_2D( 1, 1, 1, 1 )
+         ikt  = mbkt(ji,jj)
+         zdep = e3t(ji,jj,ikt,Kmm) / xstep
+         zwsbio4(ji,jj) = MIN( 0.99 * zdep, wsbio4(ji,jj,ikt) )
+         zwsbio3(ji,jj) = MIN( 0.99 * zdep, wsbio3(ji,jj,ikt) )
+      END_2D
+
       IF( .NOT.lk_sed ) THEN
-         ! OA: Warning, the following part is necessary to avoid CFL problems above the sediments
-         ! --------------------------------------------------------------------
-         DO_2D( 1, 1, 1, 1 )
-            ikt  = mbkt(ji,jj)
-            zdep = e3t(ji,jj,ikt,Kmm) / xstep
-            zwsbio4(ji,jj) = MIN( 0.99 * zdep, wsbio4(ji,jj,ikt) )
-            zwsbio3(ji,jj) = MIN( 0.99 * zdep, wsbio3(ji,jj,ikt) )
-         END_2D
-
-         ! Computation of the sediment denitrification proportion: The metamodel from midlleburg (2006) is being used
-         ! Computation of the fraction of organic matter that is permanently buried from Dunne's model
+
+         ! Computation of the sediment denitrification proportion: The metamodel 
+         ! from Middleburg (2006) is used
+         ! Computation of the fraction of organic matter that is permanently 
+         ! buried from Dunne's model (2007)
          ! -------------------------------------------------------
          DO_2D( 1, 1, 1, 1 )
@@ -125,9 +125,10 @@
       ENDIF
 
-      ! This loss is scaled at each bottom grid cell for equilibrating the total budget of silica in the ocean.
-      ! Thus, the amount of silica lost in the sediments equal the supply at the surface (dust+rivers)
-      ! ------------------------------------------------------
+      ! 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
 
+      ! Loss of bSi and CaCO3 to the sediments
       DO_2D( 1, 1, 1, 1 )
          ikt  = mbkt(ji,jj)
@@ -142,4 +143,10 @@
       !
       IF( .NOT.lk_sed ) THEN
+         ! Dissolution of CaCO3 and bSi in the sediments. This is 
+         ! instantaneous since here sediments are not explicitly 
+         ! 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 p4zbc
+         ! ---------------------------------------------------------------             
          DO_2D( 1, 1, 1, 1 )
             ikt  = mbkt(ji,jj)
@@ -160,4 +167,5 @@
       ENDIF
       !
+      ! Loss of particulate organic carbon and Fe to the sediments
       DO_2D( 1, 1, 1, 1 )
          ikt  = mbkt(ji,jj)
@@ -171,4 +179,5 @@
       END_2D
       !
+      ! Loss of particulate organic N and P to the sediments (p5z)
       IF( ln_p5z ) THEN
          DO_2D( 1, 1, 1, 1 )
@@ -185,6 +194,9 @@
 
       IF( .NOT.lk_sed ) THEN
+         ! Degradation of organic matter in the sediments. The metamodel of 
+         ! Middleburg (2006) is used here to mimic the diagenetic reactions. 
          ! The 0.5 factor in zpdenit is to avoid negative NO3 concentration after
          ! denitrification in the sediments. Not very clever, but simpliest option.
+         ! ------------------------------------------------------------------------
          DO_2D( 1, 1, 1, 1 )
             ikt  = mbkt(ji,jj)
@@ -192,10 +204,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 = tr(ji,jj,ikt,jpgoc,Kbb) * zws4 + tr(ji,jj,ikt,jppoc,Kbb) * zws3
+            ! Denitrification in the sediments            
             zpdenit  = MIN( 0.5 * ( tr(ji,jj,ikt,jpno3,Kbb) - 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( ( tr(ji,jj,ikt,jpoxy,Kbb) - 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           
             tr(ji,jj,ikt,jpdoc,Krhs) = tr(ji,jj,ikt,jpdoc,Krhs) + z1pdenit - zolimit
+            ! Update of the tracers concentrations            
             tr(ji,jj,ikt,jppo4,Krhs) = tr(ji,jj,ikt,jppo4,Krhs) + zpdenit + zolimit
             tr(ji,jj,ikt,jpnh4,Krhs) = tr(ji,jj,ikt,jpnh4,Krhs) + zpdenit + zolimit
@@ -206,4 +225,5 @@
             sdenit(ji,jj) = rdenit * zpdenit * e3t(ji,jj,ikt,Kmm)
             zsedc(ji,jj)   = (1. - zrivno3) * zwstpoc * e3t(ji,jj,ikt,Kmm)
+            ! PISCES-QUOTA (p5z)            
             IF( ln_p5z ) THEN
                zwstpop              = tr(ji,jj,ikt,jpgop,Kbb) * zws4 + tr(ji,jj,ikt,jppop,Kbb) * zws3
@@ -215,18 +235,37 @@
        ENDIF
 
-
-      ! Nitrogen fixation process
-      ! Small source iron from particulate inorganic iron
-      !-----------------------------------
-      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) )
-      ENDDO
+      ! Nitrogen fixation process : light limitation of diazotrophy
+      ! Small source of iron from particulate inorganic iron (photochemistry)
+      ! This is a purely adhoc param.
+      !----------------------------------------------------------------------
+
+      ! Diazotrophy (nitrogen fixation) is modeled according to an empirical
+      ! formulation. This is described in Aumont et al. (2015). Limitation 
+      ! by P and Fe is computed. Inhibition by high N concentrations is imposed.
+      ! Diazotrophy sensitivity to temperature is parameterized as in 
+      ! Ye et al. (2012)  
+      ! ------------------------------------------------------------------------    
       IF( ln_p4z ) THEN
+         ! PISCES part              
          DO_3D( 1, 1, 1, 1, 1, jpkm1 )
-            !                      ! Potential nitrogen fixation dependant on temperature and iron
+            ! Potential nitrogen fixation dependant on temperature and iron
             ztemp = ts(ji,jj,jk,jp_tem,Kmm)
-            zmudia = MAX( 0.,-0.001096*ztemp**2 + 0.057*ztemp -0.637 ) * 7.625
-            !       Potential nitrogen fixation dependant on temperature and iron
+            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
+            ! Nitrogen fixation limitation by PO4 and Fe            
+            ztrfer = biron(ji,jj,jk) / ( concfediaz + biron(ji,jj,jk) )
+            ztrpo4 = tr(ji,jj,jk,jppo4,Kbb) / ( 1E-6 + tr(ji,jj,jk,jppo4,Kbb) )
+            zlight =  ( 1.- EXP( -etot_ndcy(ji,jj,jk) / diazolight ) ) * ( 1. - fr_i(ji,jj) ) 
+            nitrpot(ji,jj,jk) =  zmudia * r1_rday * zfact * MIN( ztrfer, ztrpo4 ) * zlight
+         END_3D
+      ELSE      
+         ! PISCES-QUOTA part
+         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+            !  Potential nitrogen fixation dependant on temperature and iron
+            ztemp = ts(ji,jj,jk,jp_tem,Kmm)
+            zmudia = MAX( 0.,-0.001096*ztemp**2 + 0.057*ztemp -0.637 ) / rno3
+            ! Nitrogen fixation is inhibited when enough NO3 and/or NH4
             xdianh4 = tr(ji,jj,jk,jpnh4,Kbb) / ( concnnh4 + tr(ji,jj,jk,jpnh4,Kbb) )
             xdiano3 = tr(ji,jj,jk,jpno3,Kbb) / ( concnno3 + tr(ji,jj,jk,jpno3,Kbb) ) * (1. - xdianh4)
@@ -234,25 +273,10 @@
             IF( zlim <= 0.1 )   zlim = 0.01
             zfact = zlim * rfact2
+            ! Nitrogen fixation limitation by PO4/DOP and Fe            
             ztrfer = biron(ji,jj,jk) / ( concfediaz + biron(ji,jj,jk) )
-            ztrpo4(ji,jj,jk) = tr(ji,jj,jk,jppo4,Kbb) / ( 1E-6 + tr(ji,jj,jk,jppo4,Kbb) )
-            ztrdp = ztrpo4(ji,jj,jk)
-            nitrpot(ji,jj,jk) =  zmudia * r1_rday * zfact * MIN( ztrfer, ztrdp ) * zlight(ji,jj,jk)
-         END_3D
-      ELSE       ! p5z
-         DO_3D( 1, 1, 1, 1, 1, jpkm1 )
-            !                      ! Potential nitrogen fixation dependant on temperature and iron
-            ztemp = ts(ji,jj,jk,jp_tem,Kmm)
-            zmudia = MAX( 0.,-0.001096*ztemp**2 + 0.057*ztemp -0.637 ) * 7.625
-            !       Potential nitrogen fixation dependant on temperature and iron
-            xdianh4 = tr(ji,jj,jk,jpnh4,Kbb) / ( concnnh4 + tr(ji,jj,jk,jpnh4,Kbb) )
-            xdiano3 = tr(ji,jj,jk,jpno3,Kbb) / ( concnno3 + tr(ji,jj,jk,jpno3,Kbb) ) * (1. - xdianh4)
-            zlim = ( 1.- xdiano3 - xdianh4 )
-            IF( zlim <= 0.1 )   zlim = 0.01
-            zfact = zlim * rfact2
-            ztrfer = biron(ji,jj,jk) / ( concfediaz + biron(ji,jj,jk) )
-            ztrpo4(ji,jj,jk) = tr(ji,jj,jk,jppo4,Kbb) / ( 1E-6 + tr(ji,jj,jk,jppo4,Kbb) )
-            ztrdop(ji,jj,jk) = tr(ji,jj,jk,jpdop,Kbb) / ( 1E-6 + tr(ji,jj,jk,jpdop,Kbb) ) * (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)
+            ztrpo4 = tr(ji,jj,jk,jppo4,Kbb) / ( 1E-6 + tr(ji,jj,jk,jppo4,Kbb) )
+            ztrdop = tr(ji,jj,jk,jpdop,Kbb) / ( 1E-6 + tr(ji,jj,jk,jpdop,Kbb) ) * (1. - ztrpo4)
+            zlight =  ( 1.- EXP( -etot_ndcy(ji,jj,jk) / diazolight ) ) * ( 1. - fr_i(ji,jj) ) 
+            nitrpot(ji,jj,jk) =  zmudia * r1_rday * zfact * MIN( ztrfer, ztrpo4 + ztrdop ) * zlight
          END_3D
       ENDIF
@@ -261,6 +285,10 @@
       ! ----------------------------------------
       IF( ln_p4z ) THEN
+         ! PISCES part              
          DO_3D( 1, 1, 1, 1, 1, jpkm1 )
             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 to POC/GOC as 
+            ! a result of mortality by predation. Completely adhoc param            
             tr(ji,jj,jk,jpnh4,Krhs) = tr(ji,jj,jk,jpnh4,Krhs) + zfact / 3.0
             tr(ji,jj,jk,jptal,Krhs) = tr(ji,jj,jk,jptal,Krhs) + rno3 * zfact / 3.0
@@ -270,23 +298,33 @@
             tr(ji,jj,jk,jpgoc,Krhs) = tr(ji,jj,jk,jpgoc,Krhs) + zfact * 1.0 / 3.0 * 1.0 / 3.0
             tr(ji,jj,jk,jpoxy,Krhs) = tr(ji,jj,jk,jpoxy,Krhs) + ( o2ut + o2nit ) * zfact * 2.0 / 3.0 + o2nit * zfact / 3.0
+            ! Fe/c of diazotrophs is assumed to be 30umol Fe/mol C at max            
+            zlight  =  ( 1.- EXP( -etot_ndcy(ji,jj,jk) / diazolight ) ) * ( 1. - fr_i(ji,jj) ) 
+            zsoufer = zlight * 2E-11 / ( 2E-11 + biron(ji,jj,jk) )
             tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) - 30E-6 * zfact * 1.0 / 3.0
             tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) + 30E-6 * zfact * 1.0 / 3.0 * 2.0 / 3.0
             tr(ji,jj,jk,jpbfe,Krhs) = tr(ji,jj,jk,jpbfe,Krhs) + 30E-6 * zfact * 1.0 / 3.0 * 1.0 / 3.0
-            tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + 0.002 * 4E-10 * zsoufer(ji,jj,jk) * rfact2 / rday
+            tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + 0.002 * 4E-10 * zsoufer * rfact2 / rday
             tr(ji,jj,jk,jppo4,Krhs) = tr(ji,jj,jk,jppo4,Krhs) + concdnh4 / ( concdnh4 + tr(ji,jj,jk,jppo4,Kbb) ) &
             &                     * 0.001 * tr(ji,jj,jk,jpdoc,Kbb) * xstep
          END_3D
       ELSE    ! p5z
+         ! PISCES-QUOTA part              
          DO_3D( 1, 1, 1, 1, 1, jpkm1 )
             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 
+            ! a result of mortality by predation. Completely adhoc param             
             tr(ji,jj,jk,jpnh4,Krhs) = tr(ji,jj,jk,jpnh4,Krhs) + zfact / 3.0
             tr(ji,jj,jk,jptal,Krhs) = tr(ji,jj,jk,jptal,Krhs) + rno3 * zfact / 3.0
+            ! N/P ratio of diazotrophs is supposed to be 46    
+            ztrpo4 = tr(ji,jj,jk,jppo4,Kbb) / ( 1E-6 + tr(ji,jj,jk,jppo4,Kbb) )
+            ztrdop = tr(ji,jj,jk,jpdop,Kbb) / ( 1E-6 + tr(ji,jj,jk,jpdop,Kbb) ) * (1. - ztrpo4)
+        
             tr(ji,jj,jk,jppo4,Krhs) = tr(ji,jj,jk,jppo4,Krhs) - 16.0 / 46.0 * zfact * ( 1.0 - 1.0 / 3.0 ) &
-            &                     * ztrpo4(ji,jj,jk) / (ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk) + rtrn)
+            &                     * ztrpo4 / (ztrpo4 + ztrdop + rtrn)
             tr(ji,jj,jk,jpdon,Krhs) = tr(ji,jj,jk,jpdon,Krhs) + zfact * 1.0 / 3.0
             tr(ji,jj,jk,jpdoc,Krhs) = tr(ji,jj,jk,jpdoc,Krhs) + zfact * 1.0 / 3.0
             tr(ji,jj,jk,jpdop,Krhs) = tr(ji,jj,jk,jpdop,Krhs) + 16.0 / 46.0 * zfact / 3.0  &
-            &                     - 16.0 / 46.0 * zfact * ztrdop(ji,jj,jk)   &
-            &                     / (ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk) + rtrn)
+            &                     - 16.0 / 46.0 * zfact * ztrdop / (ztrpo4 + ztrdop + rtrn)
             tr(ji,jj,jk,jppoc,Krhs) = tr(ji,jj,jk,jppoc,Krhs) + zfact * 1.0 / 3.0 * 2.0 / 3.0
             tr(ji,jj,jk,jppon,Krhs) = tr(ji,jj,jk,jppon,Krhs) + zfact * 1.0 / 3.0 * 2.0 /3.0
@@ -296,8 +334,11 @@
             tr(ji,jj,jk,jpgop,Krhs) = tr(ji,jj,jk,jpgop,Krhs) + 16.0 / 46.0 * zfact * 1.0 / 3.0 * 1.0 /3.0
             tr(ji,jj,jk,jpoxy,Krhs) = tr(ji,jj,jk,jpoxy,Krhs) + ( o2ut + o2nit ) * zfact * 2.0 / 3.0 + o2nit * zfact / 3.0
+            ! Fe/c of diazotrophs is assumed to be 30umol Fe/mol C at max           
+            zlight  =  ( 1.- EXP( -etot_ndcy(ji,jj,jk) / diazolight ) ) * ( 1. - fr_i(ji,jj) ) 
+            zsoufer = zlight * 2E-11 / ( 2E-11 + biron(ji,jj,jk) )
             tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) - 30E-6 * zfact * 1.0 / 3.0 
             tr(ji,jj,jk,jpsfe,Krhs) = tr(ji,jj,jk,jpsfe,Krhs) + 30E-6 * zfact * 1.0 / 3.0 * 2.0 / 3.0
             tr(ji,jj,jk,jpbfe,Krhs) = tr(ji,jj,jk,jpbfe,Krhs) + 30E-6 * zfact * 1.0 / 3.0 * 1.0 / 3.0
-            tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + 0.002 * 4E-10 * zsoufer(ji,jj,jk) * rfact2 / rday
+            tr(ji,jj,jk,jpfer,Krhs) = tr(ji,jj,jk,jpfer,Krhs) + 0.002 * 4E-10 * zsoufer * rfact2 / rday
          END_3D
          !
@@ -319,6 +360,4 @@
       ENDIF
       !
-      IF( ln_p5z )    DEALLOCATE( ztrpo4, ztrdop )
-      !
       IF( ln_timing )  CALL timing_stop('p4z_sed')
       !
