Index: /trunk/NEMO/TOP_SRC/SMS/p4zbio.kriest.h
===================================================================
--- /trunk/NEMO/TOP_SRC/SMS/p4zbio.kriest.h	(revision 618)
+++ /trunk/NEMO/TOP_SRC/SMS/p4zbio.kriest.h	(revision 618)
@@ -0,0 +1,93 @@
+
+
+CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+CCC p4zbio : PISCES MODEL  - Kriest parameterization
+CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+
+      DO jk = 1,jpkm1
+        DO jj = 1,jpj
+          DO ji = 1,jpi
+C    
+C     Evolution of DOC
+C     ----------------
+C     
+          trn(ji,jj,jk,jpdoc) = trn(ji,jj,jk,jpdoc)
+     &      +orem(ji,jj,jk)+excret2*prorca2(ji,jj,jk)
+     &      +excret*prorca(ji,jj,jk)-olimi(ji,jj,jk)-denitr(ji,jj,jk)
+     &      +grarem(ji,jj,jk)*(1.-sigma1)+grarem2(ji,jj,jk)
+     &      *(1.-sigma2)-xaggdoc(ji,jj,jk)
+          END DO
+        END DO
+      END DO
+
+
+      DO jk = 1,jpkm1
+        DO jj = 1,jpj
+          DO ji = 1,jpi
+C     
+C     Evolution of Detritus
+C     ---------------------
+C  
+            znumpoc=trn(ji,jj,jk,jpnum)/(trn(ji,jj,jk,jppoc)+rtrn)
+            znumdoc=znumpoc
+
+            trn(ji,jj,jk,jppoc) = trn(ji,jj,jk,jppoc)
+     &          -grazpoc(ji,jj,jk)
+     &          +grapoc(ji,jj,jk)+grapoc2(ji,jj,jk)
+     &          -grazm(ji,jj,jk)+tortz2(ji,jj,jk)
+     &          +respz(ji,jj,jk)+respz2(ji,jj,jk)
+     &          +respp(ji,jj,jk)+respp2(ji,jj,jk)
+     &          +tortp2(ji,jj,jk)+tortz(ji,jj,jk)
+     &          +tortp(ji,jj,jk)-orem(ji,jj,jk)
+     &          +xaggdoc(ji,jj,jk)-grazffe(ji,jj,jk)
+C    
+C     Evolution of number of aggregates
+C     ---------------------------------
+C
+            trn(ji,jj,jk,jpnum) = trn(ji,jj,jk,jpnum)
+     &          -xagg(ji,jj,jk)-(orem(ji,jj,jk)
+     &          +grazpoc(ji,jj,jk))*znumpoc
+     &          +(tortp(ji,jj,jk))*xkr_nnano
+     &          +(tortp2(ji,jj,jk)+respp(ji,jj,jk)+tortz(ji,jj,jk)
+     &          +grapoc(ji,jj,jk)-grazm(ji,jj,jk)
+     &          +respz(ji,jj,jk))*xkr_ndiat
+     &          +(grapoc2(ji,jj,jk)+tortz2(ji,jj,jk)
+     &          +respz2(ji,jj,jk))*xkr_nmeso
+     &          +respp2(ji,jj,jk)*xkr_naggr
+     &          +xaggdoc(ji,jj,jk)*znumdoc
+     &          -grazffe(ji,jj,jk)*znumpoc*wsbio4(ji,jj,jk)
+     &          /(wsbio3(ji,jj,jk)+rtrn)
+C
+            trn(ji,jj,jk,jpnum) = max(trn(ji,jj,jk,jpnum),
+     &          trn(ji,jj,jk,jppoc)/xkr_massp/xnumm(jk))
+C
+            trn(ji,jj,jk,jpnum) = min(trn(ji,jj,jk,jpnum),
+     &          trn(ji,jj,jk,jppoc)/xkr_massp/1.1)
+C
+
+          END DO
+        END DO
+      END DO
+
+
+
+      DO jk = 1,jpkm1
+        DO jj = 1,jpj
+          DO ji = 1,jpi
+C
+C     Evolution of biogenic Iron
+C     --------------------------
+C
+          trn(ji,jj,jk,jpsfe) = trn(ji,jj,jk,jpsfe)
+     &     +unass*(grazpf(ji,jj,jk)+grazsf(ji,jj,jk))
+     &     -(1.-unass2)*grazpof(ji,jj,jk)-(1.-unass)*grazmf(ji,jj,jk)
+     &     -(1.-unass2)*grazfff(ji,jj,jk)+unass2*(graznf(ji,jj,jk)
+     &     +grazf(ji,jj,jk)+ferat3*grazz(ji,jj,jk))+ferat3
+     &     *(tortz2(ji,jj,jk)+respz2(ji,jj,jk)+tortz(ji,jj,jk)
+     &     +respz(ji,jj,jk))-ofer(ji,jj,jk)+(respnf(ji,jj,jk)
+     &     +tortnf(ji,jj,jk))+tortdf(ji,jj,jk)+respdf(ji,jj,jk)
+     &     +xbactfer(ji,jj,jk)+xscave(ji,jj,jk)*zdenom1(ji,jj,jk)
+          END DO
+        END DO
+      END DO
+
Index: /trunk/NEMO/TOP_SRC/SMS/p4zbio.std.h
===================================================================
--- /trunk/NEMO/TOP_SRC/SMS/p4zbio.std.h	(revision 618)
+++ /trunk/NEMO/TOP_SRC/SMS/p4zbio.std.h	(revision 618)
@@ -0,0 +1,80 @@
+CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+CCC p4zbio : PISCES MODEL  - Standard parameterization
+CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+
+      DO jk = 1,jpkm1
+        DO jj = 1,jpj
+          DO ji = 1,jpi
+C    
+C     Evolution of DOC
+C     ----------------
+C     
+          trn(ji,jj,jk,jpdoc) = trn(ji,jj,jk,jpdoc)
+     &      +orem(ji,jj,jk)+excret2*prorca2(ji,jj,jk)
+     &      +excret*prorca(ji,jj,jk)-olimi(ji,jj,jk)-denitr(ji,jj,jk)
+     &      +grarem(ji,jj,jk)*(1.-sigma1)+grarem2(ji,jj,jk)
+     &      *(1.-sigma2)-xaggdoc(ji,jj,jk)-xaggdoc2(ji,jj,jk)
+          END DO
+        END DO
+      END DO
+
+      DO jk = 1,jpkm1
+        DO jj = 1,jpj
+          DO ji = 1,jpi
+C     
+C     Evolution of Detritus
+C     ---------------------
+C     
+          trn(ji,jj,jk,jppoc) = trn(ji,jj,jk,jppoc)
+     &     -grazpoc(ji,jj,jk)+grapoc(ji,jj,jk)-grazm(ji,jj,jk)
+     &     +respz(ji,jj,jk)-xagg(ji,jj,jk)+xaggdoc(ji,jj,jk)
+     &     +(1.-0.5*zfracal(ji,jj,jk))*(tortp(ji,jj,jk)
+     &     +respp(ji,jj,jk))+0.5*tortp2(ji,jj,jk)
+     &     +orem2(ji,jj,jk)+tortz(ji,jj,jk)-orem(ji,jj,jk)
+C    
+C     Evolution of rapid Detritus
+C     ---------------------
+C    
+          trn(ji,jj,jk,jpgoc) = trn(ji,jj,jk,jpgoc)
+     &     +grapoc2(ji,jj,jk)+respp2(ji,jj,jk)+xagg(ji,jj,jk)
+     &     +tortz2(ji,jj,jk)+respz2(ji,jj,jk)-orem2(ji,jj,jk)
+     &     +0.5*zfracal(ji,jj,jk)*(respp(ji,jj,jk)+tortp(ji,jj,jk))
+     &     +0.5*tortp2(ji,jj,jk)+xaggdoc2(ji,jj,jk)-grazffe(ji,jj,jk)
+C
+          END DO
+        END DO
+      END DO
+
+
+
+      DO jk = 1,jpkm1
+        DO jj = 1,jpj
+          DO ji = 1,jpi
+C
+C     Evolution of small biogenic Iron
+C     --------------------------
+C
+          trn(ji,jj,jk,jpsfe) = trn(ji,jj,jk,jpsfe)
+     &     +unass*(grazpf(ji,jj,jk)+grazsf(ji,jj,jk))
+     &     -grazpof(ji,jj,jk)-(1.-unass)*grazmf(ji,jj,jk)
+     &     +(1.-0.5*zfracal(ji,jj,jk))*(tortnf(ji,jj,jk)
+     &     +respnf(ji,jj,jk))+0.5*tortdf(ji,jj,jk)+ferat3*
+     &     (tortz(ji,jj,jk)+respz(ji,jj,jk))-ofer(ji,jj,jk)
+     &     +ofer2(ji,jj,jk)-xaggfe(ji,jj,jk)
+     &     +xscave(ji,jj,jk)*zdenom1(ji,jj,jk)
+C
+C     Evolution of big biogenic Iron
+C     --------------------------
+C
+          trn(ji,jj,jk,jpbfe) = trn(ji,jj,jk,jpbfe)
+     &     +unass2*(graznf(ji,jj,jk)+grazf(ji,jj,jk)+grazfff(ji,jj,jk)
+     &     +grazpof(ji,jj,jk)+ferat3*grazz(ji,jj,jk))+ferat3*
+     &     (tortz2(ji,jj,jk)+respz2(ji,jj,jk))-ofer2(ji,jj,jk)
+     &     +0.5*zfracal(ji,jj,jk)*(respnf(ji,jj,jk)+tortnf(ji,jj,jk))
+     &     +0.5*tortdf(ji,jj,jk)+respdf(ji,jj,jk)+xaggfe(ji,jj,jk)
+     &     +xbactfer(ji,jj,jk)-grazfff(ji,jj,jk)+xscave(ji,jj,jk)
+     &     *zdenom2(ji,jj,jk)
+          END DO
+        END DO
+      END DO
+
Index: /trunk/NEMO/TOP_SRC/SMS/p4zsink.kriest.h
===================================================================
--- /trunk/NEMO/TOP_SRC/SMS/p4zsink.kriest.h	(revision 618)
+++ /trunk/NEMO/TOP_SRC/SMS/p4zsink.kriest.h	(revision 618)
@@ -0,0 +1,211 @@
+CCCCCC PISCES MODEL: Kriest parameterization
+CDIR$ LIST
+CC----------------------------------------------------------------------
+CC local declarations
+CC ==================
+      INTEGER jksed, ji, jj, jk
+      REAL xagg1,xagg2,xagg3,xagg4,xagg5,xaggsi,xaggsh
+      REAL znum(jpi,jpj,jpk)
+      REAL xnum,xeps,xfm,xgm,xsm
+      REAL xdiv,xdiv1,xdiv2,xdiv3,xdiv4,xdiv5
+      REAL zval1, zval2, zval3, zval4
+      REAL zstep
+#if defined key_trc_dia3d
+      REAL zrfact2
+#endif
+      REAL sinking(jpi,jpj,jpk),sinking2(jpi,jpj,jpk)
+      REAL sinkfer(jpi,jpj,jpk)
+      REAL sinkcal(jpi,jpj,jpk),sinksil(jpi,jpj,jpk)
+C
+C
+C    Time step duration for biology
+C    ------------------------------
+C
+       zstep=rfact2/rjjss
+       
+C
+C
+C     Initialisation of variables used to compute Sinking Speed
+C     ---------------------------------------------------------
+C
+       znum(:,:,:) = 0.
+       jksed = 10
+       zval1 = 1. + xkr_zeta
+       zval2 = 1. + xkr_zeta + xkr_eta
+       zval3 = 1. + xkr_eta
+C
+C     Computation of the vertical sinking speed : Kriest et Evans, 2000
+C     -----------------------------------------------------------------
+C    
+       do jk=1,jpk-1
+         do jj=1,jpj
+           do ji=1,jpi
+             IF (tmask(ji,jj,jk).NE.0) THEN
+                 xnum = trn(ji,jj,jk,jppoc) / (trn(ji,jj,jk,jpnum)+rtrn)
+     &                     / xkr_massp
+C -------------- To avoid sinking speed over 50 m/day -------
+                 xnum = min( xnumm(jk), xnum )
+                 xnum = max( 1.1, xnum )
+                 znum(ji,jj,jk) = xnum
+C------------------------------------------------------------
+                 xeps = ( zval1 * xnum - 1. )/ ( xnum - 1. )
+                 xfm  = xkr_frac**( 1. - xeps )
+                 xgm  = xkr_frac**( zval1 - xeps )
+                 xdiv = max(1E-4,abs(xeps-zval2))*sign(1.,(xeps-zval2))
+                 xdiv1=(xeps-zval3)
+                 wsbio3(ji,jj,jk)= xkr_wsbio_min * ( xeps-zval1 ) / xdiv
+     &                           - xkr_wsbio_max * xgm * xkr_eta / xdiv
+                 wsbio4(ji,jj,jk)= xkr_wsbio_min * ( xeps-1. ) / xdiv1
+     &                           - xkr_wsbio_max * xfm * xkr_eta / xdiv1
+                 IF( xnum == 1.1) THEN
+                     wsbio3(ji,jj,jk) = wsbio4(ji,jj,jk)
+                 ENDIF
+             ENDIF
+           end do
+         end do
+       end do
+C
+       wscal(:,:,:)=max(wsbio3(:,:,:),50.)
+C
+C
+C   INITIALIZE TO ZERO ALL THE SINKING ARRAYS
+C   -----------------------------------------
+C
+         sinking=0.
+         sinking2=0.
+         sinkcal=0.
+         sinkfer=0.
+         sinksil=0.
+C
+C   Compute the sedimentation term using p4zsink2 for all
+C   the sinking particles
+C   -----------------------------------------------------
+C
+         CALL p4zsink2(wsbio3,sinking,jppoc)
+         CALL p4zsink2(wsbio4,sinking2,jpnum)
+         CALL p4zsink2(wsbio3,sinkfer,jpsfe)
+         CALL p4zsink2(wscal,sinksil,jpdsi)
+         CALL p4zsink2(wscal,sinkcal,jpcal)
+
+C
+C  Exchange between organic matter compartments due to
+C  coagulation/disaggregation
+C  ---------------------------------------------------
+C
+         zval1 = 1. + xkr_zeta
+         zval2 = 1. + xkr_eta
+         zval3 = 3. + xkr_eta
+         zval4 = 4. + xkr_eta
+
+         DO jk = 1,jpkm1
+           DO jj = 1,jpj
+             DO ji = 1,jpi
+               IF (tmask(ji,jj,jk).NE.0.) THEN
+C
+                   xnum=trn(ji,jj,jk,jppoc)/(trn(ji,jj,jk,jpnum)+rtrn)
+     &                         /xkr_massp
+C -------------- To avoid sinking speed over 50 m/day -------
+                   xnum=min(xnumm(jk),xnum)
+                   xnum=max(1.1,xnum)
+C------------------------------------------------------------
+                   xeps =(zval1*xnum-1.)/(xnum-1.)
+                   xdiv =max(1E-4,abs(xeps-zval3))*sign(1.,(xeps-zval3))
+                   xdiv1=max(1E-4,abs(xeps-4.))*sign(1.,(xeps-4.))
+                   xdiv2=(xeps-2.)
+                   xdiv3=(xeps-3.)
+                   xdiv4=(xeps-zval2)
+                   xdiv5=(2*xeps-zval4)
+                   xfm=xkr_frac**(1.-xeps)
+                   xsm=xkr_frac**xkr_eta
+C
+C    Part I : Coagulation dependant on turbulence
+C    ----------------------------------------------
+C
+                   xagg1=(0.163*trn(ji,jj,jk,jpnum)**2
+     &           *2.*( (xfm-1.)*(xfm*xkr_mass_max**3-xkr_mass_min**3)
+     &           *(xeps-1)/xdiv1 + 3.*(xfm*xkr_mass_max-xkr_mass_min)
+     &           *(xfm*xkr_mass_max**2-xkr_mass_min**2)
+     &           *(xeps-1.)**2/(xdiv2*xdiv3)))
+#    if defined key_off_degrad
+     &                 *facvol(ji,jj,jk)
+#    endif
+
+C
+                   xagg2=(2*0.163*trn(ji,jj,jk,jpnum)**2*xfm*
+     &                   ((xkr_mass_max**3+3.*(xkr_mass_max**2
+     &                    *xkr_mass_min*(xeps-1.)/xdiv2
+     &                    +xkr_mass_max*xkr_mass_min**2*(xeps-1.)/xdiv3)
+     &                    +xkr_mass_min**3*(xeps-1)/xdiv1)
+     &                    -xfm*xkr_mass_max**3*(1.+3.*((xeps-1.)/
+     &                    (xeps-2.)+(xeps-1.)/xdiv3)+(xeps-1.)/xdiv1)))
+#    if defined key_off_degrad
+     &                 *facvol(ji,jj,jk)
+#    endif
+C
+                   xagg3=(0.163*trn(ji,jj,jk,jpnum)**2*xfm**2*8.
+     &                 *xkr_mass_max**3)
+#    if defined key_off_degrad
+     &                 *facvol(ji,jj,jk)
+#    endif
+C
+                   xaggsh=(xagg1+xagg2+xagg3)*rfact2*zdiss(ji,jj,jk)
+     &                 /1000.
+C
+C    Aggregation of small into large particles
+C    Part II : Differential settling
+C    ----------------------------------------------
+C
+                   xagg4=(2.*3.141*0.125*trn(ji,jj,jk,jpnum)**2*
+     &                 xkr_wsbio_min*(xeps-1.)**2
+     &                 *(xkr_mass_min**2*((1.-xsm*xfm)/(xdiv3*xdiv4)
+     &                 -(1.-xfm)/(xdiv*(xeps-1.)))-
+     &                 ((xfm*xfm*xkr_mass_max**2*xsm-xkr_mass_min**2)
+     &                 *xkr_eta)/(xdiv*xdiv3*xdiv5)))
+#    if defined key_off_degrad
+     &                 *facvol(ji,jj,jk)
+#    endif
+C
+                   xagg5=(2.*3.141*0.125*trn(ji,jj,jk,jpnum)**2
+     &                 *(xeps-1.)*xfm*xkr_wsbio_min
+     &                 *(xsm*(xkr_mass_min**2-xfm*xkr_mass_max**2)
+     &                 /xdiv3-(xkr_mass_min**2-xfm*xsm*xkr_mass_max**2)
+     &                 /xdiv))
+#    if defined key_off_degrad
+     &                 *facvol(ji,jj,jk)
+#    endif
+C
+                   xaggsi=(xagg4+xagg5)*zstep/10.
+C
+                   xagg(ji,jj,jk)=0.5 * xkr_stick*(xaggsh+xaggsi)
+C
+C     Aggregation of DOC to small particles
+C     --------------------------------------
+C
+                   xaggdoc(ji,jj,jk)=(0.4*trn(ji,jj,jk,jpdoc)
+     &                 +1018.*trn(ji,jj,jk,jppoc))*zstep
+     &                 *zdiss(ji,jj,jk)*trn(ji,jj,jk,jpdoc)
+#    if defined key_off_degrad
+     &                 *facvol(ji,jj,jk)
+#    endif
+C
+
+               ENDIF
+             END DO
+           END DO
+         END DO
+C
+#    if defined key_trc_dia3d
+         zrfact2 = 1.e3*rfact2r
+         trc2d(:,:,5) = sinking(:,:,jksed+1)*zrfact2
+         trc2d(:,:,6) = sinking2(:,:,jksed+1)*zrfact2
+         trc2d(:,:,7) = sinkfer(:,:,jksed+1)*zrfact2
+         trc2d(:,:,9) = sinksil(:,:,jksed+1)*zrfact2
+         trc2d(:,:,10) = sinkcal(:,:,jksed+1)*zrfact2
+         trc3d(:,:,:,12) = sinking(:,:,:)*zrfact2
+         trc3d(:,:,:,13) = sinking2(:,:,:)*zrfact2
+         trc3d(:,:,:,14) = sinksil(:,:,:)*zrfact2
+         trc3d(:,:,:,15) = sinkcal(:,:,:)*zrfact2
+         trc3d(:,:,:,16) = znum(:,:,:)
+         trc3d(:,:,:,17) = wsbio3(:,:,:)
+         trc3d(:,:,:,18) = wsbio4(:,:,:)
+#    endif
Index: /trunk/NEMO/TOP_SRC/SMS/p4zsink.std.h
===================================================================
--- /trunk/NEMO/TOP_SRC/SMS/p4zsink.std.h	(revision 618)
+++ /trunk/NEMO/TOP_SRC/SMS/p4zsink.std.h	(revision 618)
@@ -0,0 +1,143 @@
+
+CDIR$ LIST
+CC----------------------------------------------------------------------
+CC local declarations
+CC ==================
+      INTEGER jksed, ji, jj, jk
+      REAL xagg1,xagg2,xagg3,xagg4
+      REAL zfact,zstep,wsmax
+#if defined key_trc_dia3d
+      REAL zrfact2
+#endif
+      REAL sinking(jpi,jpj,jpk),sinking2(jpi,jpj,jpk)
+      REAL sinkfer(jpi,jpj,jpk),sinkfer2(jpi,jpj,jpk)
+      REAL sinkcal(jpi,jpj,jpk),sinksil(jpi,jpj,jpk)
+C
+C    Time step duration for biology
+C    ------------------------------
+C
+       zstep=rfact2/rjjss
+C
+C    Sinking speeds of detritus is increased with depth as shown
+C    by data and from the coagulation theory
+C    -----------------------------------------------------------
+C
+       jksed=10
+C
+       DO jk=1,jpk-1
+         DO jj=1,jpj
+           DO ji=1,jpi
+       zfact=max(0.,fsdepw(ji,jj,jk+1)-hmld(ji,jj))/2000.
+       wsbio4(ji,jj,jk)=wsbio2+(200.-wsbio2)*zfact
+             END DO
+           END DO
+         END DO
+C
+C      LIMIT THE VALUES OF THE SINKING SPEEDS 
+C      TO AVOID NUMERICAL INSTABILITIES
+C
+      wsbio3(:,:,:)=wsbio
+
+      Do jk=1,jpk-1
+        DO jj=1,jpj
+          DO ji=1,jpi
+       wsmax=0.8*fse3t(ji,jj,jk)/zstep
+       wsbio4(ji,jj,jk)=min(wsbio4(ji,jj,jk),wsmax)
+       wsbio3(ji,jj,jk)=min(wsbio3(ji,jj,jk),wsmax)
+          END DO
+        END DO
+       END DO
+
+      wscal(:,:,:)=wsbio4(:,:,:)
+C
+C
+C   INITIALIZE TO ZERO ALL THE SINKING ARRAYS
+C   -----------------------------------------
+C
+         sinking=0.
+         sinking2=0.
+         sinkcal=0.
+         sinkfer=0.
+         sinksil=0.
+         sinkfer2=0.
+C
+C   Compute the sedimentation term using p4zsink2 for all
+C   the sinking particles
+C   -----------------------------------------------------
+C
+         CALL p4zsink2(wsbio3,sinking,jppoc)
+         CALL p4zsink2(wsbio3,sinkfer,jpsfe)
+         CALL p4zsink2(wsbio4,sinking2,jpgoc)
+         CALL p4zsink2(wsbio4,sinkfer2,jpbfe)
+         CALL p4zsink2(wsbio4,sinksil,jpdsi)
+         CALL p4zsink2(wscal,sinkcal,jpcal)
+C
+C  Exchange between organic matter compartments due to
+C  coagulation/disaggregation
+C  ---------------------------------------------------
+C
+         DO jk = 1,jpkm1
+           DO jj = 1,jpj
+             DO ji = 1,jpi
+C
+        zfact=zstep*zdiss(ji,jj,jk)
+C
+C    Part I : Coagulation dependent on turbulence
+C    ----------------------------------------------
+C
+         xagg1=940.*zfact*trn(ji,jj,jk,jppoc)*trn(ji,jj,jk,jppoc)
+#    if defined key_off_degrad
+     &     *facvol(ji,jj,jk)
+#    endif
+
+         xagg2=1.054E4*zfact*trn(ji,jj,jk,jppoc)*trn(ji,jj,jk,jpgoc)
+#    if defined key_off_degrad
+     &     *facvol(ji,jj,jk)
+#    endif
+C
+C    Aggregation of small into large particles
+C    Part II : Differential settling
+C    ----------------------------------------------
+C
+         xagg3=0.66*zstep*trn(ji,jj,jk,jppoc)*trn(ji,jj,jk,jppoc)
+#    if defined key_off_degrad
+     &     *facvol(ji,jj,jk)
+#    endif
+
+         xagg4=0.*zstep*trn(ji,jj,jk,jppoc)*trn(ji,jj,jk,jpgoc)
+#    if defined key_off_degrad
+     &     *facvol(ji,jj,jk)
+#    endif
+C
+         xagg(ji,jj,jk)=xagg1+xagg2+xagg3+xagg4
+         xaggfe(ji,jj,jk)=xagg(ji,jj,jk)*trn(ji,jj,jk,jpsfe)/
+     &     (trn(ji,jj,jk,jppoc)+rtrn)
+C
+C     Aggregation of DOC to small particles
+C     --------------------------------------
+C
+         xaggdoc(ji,jj,jk)=(80*trn(ji,jj,jk,jpdoc)+698.
+     &     *trn(ji,jj,jk,jppoc))*zfact*trn(ji,jj,jk,jpdoc)
+#    if defined key_off_degrad
+     &     *facvol(ji,jj,jk)
+#    endif
+
+         xaggdoc2(ji,jj,jk)=1.05E4*zfact*trn(ji,jj,jk,jpgoc)
+     &     *trn(ji,jj,jk,jpdoc)
+#    if defined key_off_degrad
+     &     *facvol(ji,jj,jk)
+#    endif
+C
+             END DO
+           END DO
+         END DO
+
+#    if defined key_trc_dia3d
+         zrfact2 = 1.e3*rfact2r
+         trc2d(:,:,5) = sinking(:,:,jksed+1)*zrfact2
+         trc2d(:,:,6) = sinking2(:,:,jksed+1)*zrfact2
+         trc2d(:,:,7) = sinkfer(:,:,jksed+1)*zrfact2
+         trc2d(:,:,8) = sinkfer2(:,:,jksed+1)*zrfact2
+         trc2d(:,:,9) = sinksil(:,:,jksed+1)*zrfact2
+         trc2d(:,:,10) = sinkcal(:,:,jksed+1)*zrfact2
+#    endif
