Index: branches/2016/dev_r7012_ROBUST5_CNRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zagg.F90
===================================================================
--- branches/2016/dev_r7012_ROBUST5_CNRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zagg.F90	(revision 7162)
+++ branches/2016/dev_r7012_ROBUST5_CNRS/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zagg.F90	(revision 7176)
@@ -8,4 +8,5 @@
    !!             3.4  !  2011-06  (O. Aumont, C. Ethe) Change aggregation formula
    !!             3.5  !  2012-07  (O. Aumont) Introduce potential time-splitting
+   !!             3.6  !  2015-05  (O. Aumont) PISCES quota
    !!----------------------------------------------------------------------
    !!----------------------------------------------------------------------
@@ -39,7 +40,10 @@
       INTEGER, INTENT(in) :: kt, knt
       INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   zagg1, zagg2, zagg3, zagg4
-      REAL(wp) ::   zagg , zaggfe, zaggdoc, zaggdoc2, zaggdoc3
-      REAL(wp) ::   zfact
+      REAL(wp) ::   zagg, zagg1, zagg2, zagg3, zagg4
+      REAL(wp) ::   zaggpoc1, zaggpoc2, zaggpoc3, zaggpoc4
+      REAL(wp) ::   zaggpoc , zaggfe, zaggdoc, zaggdoc2, zaggdoc3
+      REAL(wp) ::   zaggpon , zaggdon, zaggdon2, zaggdon3
+      REAL(wp) ::   zaggpop, zaggdop, zaggdop2, zaggdop3
+      REAL(wp) ::   zaggtmp, zfact, zmax
       CHARACTER (len=25) :: charout
       !!---------------------------------------------------------------------
@@ -50,48 +54,121 @@
       !  Exchange between organic matter compartments due to coagulation/disaggregation
       !  ---------------------------------------------------
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               !
-               zfact = xstep * xdiss(ji,jj,jk)
-               !  Part I : Coagulation dependent on turbulence
-               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)
+      IF( ln_p4z ) THEN
+         !
+         DO jk = 1, jpkm1
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  !
+                  zfact = xstep * xdiss(ji,jj,jk)
+                  !  Part I : Coagulation dependent on turbulence
+                  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)
 
-               ! Part II : Differential settling
+                  ! Part II : Differential settling
 
-               !  Aggregation of small into large particles
-               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)
+                  !  Aggregation of small into large particles
+                  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)
 
-               zagg   = zagg1 + zagg2 + zagg3 + zagg4
-               zaggfe = zagg * trb(ji,jj,jk,jpsfe) / ( trb(ji,jj,jk,jppoc) + rtrn )
+                  zagg   = zagg1 + zagg2 + zagg3 + zagg4
+                  zaggfe = zagg * trb(ji,jj,jk,jpsfe) / ( trb(ji,jj,jk,jppoc) + rtrn )
 
-               ! Aggregation of DOC to POC : 
-               ! 1st term is shear aggregation of DOC-DOC
-               ! 2nd term is shear aggregation of DOC-POC
-               ! 3rd term is differential settling of DOC-POC
-               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)
-               ! transfer of DOC to GOC : 
-               ! 1st term is shear aggregation
-               ! 2nd term is differential settling 
-               zaggdoc2 = ( 3.53E3 * zfact + 0.1 * xstep ) * trb(ji,jj,jk,jpgoc) * 0.3 * trb(ji,jj,jk,jpdoc)
-               ! tranfer of DOC to POC due to brownian motion
-               zaggdoc3 =  114. * 0.3 * trb(ji,jj,jk,jpdoc) *xstep * 0.3 * trb(ji,jj,jk,jpdoc)
+                  ! Aggregation of DOC to POC : 
+                  ! 1st term is shear aggregation of DOC-DOC
+                  ! 2nd term is shear aggregation of DOC-POC
+                  ! 3rd term is differential settling of DOC-POC
+                  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)
+                  ! transfer of DOC to GOC : 
+                  ! 1st term is shear aggregation
+                  ! 2nd term is differential settling 
+                  zaggdoc2 = ( 3.53E3 * zfact + 0.1 * xstep ) * trb(ji,jj,jk,jpgoc) * 0.3 * trb(ji,jj,jk,jpdoc)
+                  ! tranfer of DOC to POC due to brownian motion
+                  zaggdoc3 =  114. * 0.3 * trb(ji,jj,jk,jpdoc) *xstep * 0.3 * trb(ji,jj,jk,jpdoc)
 
-               !  Update the trends
-               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) - zagg + zaggdoc + zaggdoc3
-               tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zagg + zaggdoc2
-               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) - zaggfe
-               tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zaggfe
-               tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - zaggdoc - zaggdoc2 - zaggdoc3
-               !
-               conspoc(ji,jj,jk) = conspoc(ji,jj,jk) - zagg + zaggdoc + zaggdoc3
-               prodgoc(ji,jj,jk) = prodgoc(ji,jj,jk) + zagg + zaggdoc2
-               !
+                  !  Update the trends
+                  tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) - zagg + zaggdoc + zaggdoc3
+                  tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zagg + zaggdoc2
+                  tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) - zaggfe
+                  tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zaggfe
+                  tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - zaggdoc - zaggdoc2 - zaggdoc3
+                  !
+                  conspoc(ji,jj,jk) = conspoc(ji,jj,jk) - zagg + zaggdoc + zaggdoc3
+                  prodgoc(ji,jj,jk) = prodgoc(ji,jj,jk) + zagg + zaggdoc2
+                  !
+               END DO
             END DO
          END DO
-      END DO
+      ELSE    ! ln_p5z
+        !
+         DO jk = 1, jpkm1
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  !
+                  zfact = xstep * xdiss(ji,jj,jk)
+                  !  Part I : Coagulation dependent on turbulence
+                  zaggtmp = 25.9  * zfact * trb(ji,jj,jk,jppoc)
+                  zaggpoc1 = zaggtmp * trb(ji,jj,jk,jppoc)
+                  zaggtmp = 4452. * zfact * trb(ji,jj,jk,jpgoc)
+                  zaggpoc2 = zaggtmp * trb(ji,jj,jk,jppoc)
+
+                  ! Part II : Differential settling
+   
+                  !  Aggregation of small into large particles
+                  zaggtmp =  47.1 * xstep * trb(ji,jj,jk,jpgoc)
+                  zaggpoc3 = zaggtmp * trb(ji,jj,jk,jppoc)
+                  zaggtmp =  3.3  * xstep * trb(ji,jj,jk,jppoc)
+                  zaggpoc4 = zaggtmp * trb(ji,jj,jk,jppoc)
+
+                  zaggpoc   = zaggpoc1 + zaggpoc2 + zaggpoc3 + zaggpoc4
+                  zaggpon = zaggpoc * trb(ji,jj,jk,jppon) / ( trb(ji,jj,jk,jppoc) + rtrn)
+                  zaggpop = zaggpoc * trb(ji,jj,jk,jppop) / ( trb(ji,jj,jk,jppoc) + rtrn)
+                  zaggfe = zaggpoc * trb(ji,jj,jk,jpsfe) / ( trb(ji,jj,jk,jppoc)  + rtrn )
+
+                  ! Aggregation of DOC to POC : 
+                  ! 1st term is shear aggregation of DOC-DOC
+                  ! 2nd term is shear aggregation of DOC-POC
+                  ! 3rd term is differential settling of DOC-POC
+                  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) )
+                  zaggdoc  = zaggtmp * 0.3 * trb(ji,jj,jk,jpdoc)
+                  zaggdon  = zaggtmp * 0.3 * trb(ji,jj,jk,jpdon)
+                  zaggdop  = zaggtmp * 0.3 * trb(ji,jj,jk,jpdop)
+
+                  ! transfer of DOC to GOC : 
+                  ! 1st term is shear aggregation
+                  ! 2nd term is differential settling 
+                  zaggtmp = ( 3.53E3 * zfact + 0.1 * xstep ) * trb(ji,jj,jk,jpgoc)
+                  zaggdoc2 = zaggtmp * 0.3 * trb(ji,jj,jk,jpdoc)
+                  zaggdon2 = zaggtmp * 0.3 * trb(ji,jj,jk,jpdon)
+                  zaggdop2 = zaggtmp * 0.3 * trb(ji,jj,jk,jpdop)
+
+                  ! tranfer of DOC to POC due to brownian motion
+                  zaggtmp = ( 114. * 0.3 * trb(ji,jj,jk,jpdoc) ) * xstep
+                  zaggdoc3 =  zaggtmp * 0.3 * trb(ji,jj,jk,jpdoc)
+                  zaggdon3 =  zaggtmp * 0.3 * trb(ji,jj,jk,jpdon)
+                  zaggdop3 =  zaggtmp * 0.3 * trb(ji,jj,jk,jpdop)
+
+                  !  Update the trends
+                  tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) - zaggpoc + zaggdoc + zaggdoc3
+                  tra(ji,jj,jk,jppon) = tra(ji,jj,jk,jppon) - zaggpon + zaggdon + zaggdon3
+                  tra(ji,jj,jk,jppop) = tra(ji,jj,jk,jppop) - zaggpop + zaggdop + zaggdop3
+                  tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zaggpoc + zaggdoc2
+                  tra(ji,jj,jk,jpgon) = tra(ji,jj,jk,jpgon) + zaggpon + zaggdon2
+                  tra(ji,jj,jk,jpgop) = tra(ji,jj,jk,jpgop) + zaggpop + zaggdop2
+                  tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) - zaggfe
+                  tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zaggfe
+                  tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - zaggdoc - zaggdoc2 - zaggdoc3
+                  tra(ji,jj,jk,jpdon) = tra(ji,jj,jk,jpdon) - zaggdon - zaggdon2 - zaggdon3
+                  tra(ji,jj,jk,jpdop) = tra(ji,jj,jk,jpdop) - zaggdop - zaggdop2 - zaggdop3
+                  !
+                  conspoc(ji,jj,jk) = conspoc(ji,jj,jk) - zaggpoc + zaggdoc + zaggdoc3
+                  prodgoc(ji,jj,jk) = prodgoc(ji,jj,jk) + zaggpoc + zaggdoc2
+                  !
+               END DO
+            END DO
+         END DO
+         !
+      ENDIF
       !
       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
