Index: branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/PISCES/p4zsink.F90
===================================================================
--- branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/PISCES/p4zsink.F90	(revision 2038)
+++ branches/DEV_r2006_merge_TRA_TRC/NEMO/TOP_SRC/PISCES/p4zsink.F90	(revision 2104)
@@ -19,5 +19,6 @@
    PRIVATE
 
-   PUBLIC   p4z_sink    ! called in p4zbio.F90
+   PUBLIC   p4z_sink         ! called in p4zbio.F90
+   PUBLIC   p4z_sink_init    ! called in trcsms_pisces.F90
 
    !! * Shared module variables
@@ -31,7 +32,4 @@
      sinkcal, sinksil,    &    !: CaCO3 and BSi sinking fluxes
      sinkfer                   !: Small BFe sinking flux
-
-   REAL(wp) ::   &
-     xstep , xstep2            !: Time step duration for biology
 
    INTEGER  :: &
@@ -106,12 +104,6 @@
       !!---------------------------------------------------------------------
 
-      IF( ( kt * jnt ) == nittrc000  )  THEN 
-          CALL p4z_sink_init   ! Initialization (first time-step only)
-          xstep  = rfact2 / rday      ! Time step duration for biology
-          xstep2 = rfact2 /  2.
-      ENDIF
-
-!     Initialisation of variables used to compute Sinking Speed
-!     ---------------------------------------------------------
+      !     Initialisation of variables used to compute Sinking Speed
+      !     ---------------------------------------------------------
 
        znum3d(:,:,:) = 0.e0
@@ -120,6 +112,6 @@
        zval3 = 1. + xkr_eta
 
-!     Computation of the vertical sinking speed : Kriest et Evans, 2000
-!     -----------------------------------------------------------------
+     !     Computation of the vertical sinking speed : Kriest et Evans, 2000
+     !     -----------------------------------------------------------------
 
       DO jk = 1, jpkm1
@@ -128,9 +120,9 @@
                IF( tmask(ji,jj,jk) /= 0.e0 ) THEN
                   znum = trn(ji,jj,jk,jppoc) / ( trn(ji,jj,jk,jpnum) + rtrn ) / xkr_massp
-! -------------- To avoid sinking speed over 50 m/day -------
+                  ! -------------- To avoid sinking speed over 50 m/day -------
                   znum  = MIN( xnumm(jk), znum )
                   znum  = MAX( 1.1      , znum )
                   znum3d(ji,jj,jk) = znum
-!------------------------------------------------------------
+                  !------------------------------------------------------------
                   zeps  = ( zval1 * znum - 1. )/ ( znum - 1. )
                   zfm   = xkr_frac**( 1. - zeps )
@@ -150,7 +142,6 @@
       wscal(:,:,:) = MAX( wsbio3(:,:,:), 50. )
 
-
-!   INITIALIZE TO ZERO ALL THE SINKING ARRAYS
-!   -----------------------------------------
+      !   INITIALIZE TO ZERO ALL THE SINKING ARRAYS
+      !   -----------------------------------------
 
       sinking (:,:,:) = 0.e0
@@ -160,7 +151,6 @@
       sinksil (:,:,:) = 0.e0
 
-!   Compute the sedimentation term using p4zsink2 for all
-!   the sinking particles
-!   -----------------------------------------------------
+     !   Compute the sedimentation term using p4zsink2 for all the sinking particles
+     !   -----------------------------------------------------
 
       CALL p4z_sink2( wsbio3, sinking , jppoc )
@@ -170,7 +160,6 @@
       CALL p4z_sink2( wscal , sinkcal , jpcal )
 
-!  Exchange between organic matter compartments due to
-!  coagulation/disaggregation
-!  ---------------------------------------------------
+     !  Exchange between organic matter compartments due to coagulation/disaggregation
+     !  ---------------------------------------------------
 
       zval1 = 1. + xkr_zeta
@@ -185,8 +174,8 @@
 
                   znum = trn(ji,jj,jk,jppoc)/(trn(ji,jj,jk,jpnum)+rtrn) / xkr_massp
-! -------------- To avoid sinking speed over 50 m/day -------
+                  !-------------- To avoid sinking speed over 50 m/day -------
                   znum  = min(xnumm(jk),znum)
                   znum  = MAX( 1.1,znum)
-!------------------------------------------------------------
+                  !------------------------------------------------------------
                   zeps  = ( zval1 * znum - 1.) / ( znum - 1.)
                   zdiv  = MAX( 1.e-4, ABS( zeps - zval3) ) * SIGN( 1., zeps - zval3 )
@@ -199,6 +188,6 @@
                   zsm   = xkr_frac**xkr_eta
 
-!    Part I : Coagulation dependant on turbulence
-!    ----------------------------------------------
+                  !    Part I : Coagulation dependant on turbulence
+                  !    ----------------------------------------------
 
                   zagg1 = ( 0.163 * trn(ji,jj,jk,jpnum)**2               &
@@ -232,7 +221,7 @@
                   zaggsh = ( zagg1 + zagg2 + zagg3 ) * rfact2 * xdiss(ji,jj,jk) / 1000.
 
-!    Aggregation of small into large particles
-!    Part II : Differential settling
-!    ----------------------------------------------
+                 !    Aggregation of small into large particles
+                 !    Part II : Differential settling
+                 !    ----------------------------------------------
 
                   zagg4 = (  2.*3.141*0.125*trn(ji,jj,jk,jpnum)**2*                       &
@@ -261,6 +250,6 @@
                   zagg = 0.5 * xkr_stick * ( zaggsh + zaggsi )
 
-!     Aggregation of DOC to small particles
-!     --------------------------------------
+                  !     Aggregation of DOC to small particles
+                  !     --------------------------------------
 
                   zaggdoc = ( 0.4 * trn(ji,jj,jk,jpdoc)               &
@@ -473,5 +462,5 @@
       REAL(wp) ::   zagg1, zagg2, zagg3, zagg4
       REAL(wp) ::   zagg , zaggfe, zaggdoc, zaggdoc2
-      REAL(wp) ::   zfact, zwsmax
+      REAL(wp) ::   zfact, zwsmax, zstep
 #if defined key_diatrc
       REAL(wp) ::   zrfact2
@@ -481,16 +470,11 @@
       !!---------------------------------------------------------------------
 
-      IF( ( kt * jnt ) == nittrc000  )  THEN
-          xstep  = rfact2 / rday      ! Timestep duration for biology
-          xstep2 = rfact2 /  2.
-      ENDIF
-
-!    Sinking speeds of detritus is increased with depth as shown
-!    by data and from the coagulation theory
-!    -----------------------------------------------------------
+      !    Sinking speeds of detritus is increased with depth as shown
+      !    by data and from the coagulation theory
+      !    -----------------------------------------------------------
       DO jk = 1, jpkm1
          DO jj = 1, jpj
             DO ji=1,jpi
-               zfact = MAX( 0., fsdepw(ji,jj,jk+1)-hmld(ji,jj) ) / 4000.
+               zfact = MAX( 0., fsdepw(ji,jj,jk+1) - hmld(ji,jj) ) / 4000.
                wsbio4(ji,jj,jk) = wsbio2 + ( 200.- wsbio2 ) * zfact
             END DO
@@ -498,15 +482,13 @@
       END DO
 
-!      LIMIT THE VALUES OF THE SINKING SPEEDS 
-!      TO AVOID NUMERICAL INSTABILITIES
-
+      ! limit the values of the sinking speeds to avoid numerical instabilities  
       wsbio3(:,:,:) = wsbio
-!
-! OA Below, this is garbage. the ideal would be to find a time-splitting
-! OA algorithm that does not increase the computing cost by too much
-! OA In ROMS, I have included a time-splitting procedure. But it is 
-! OA too expensive as the loop is computed globally. Thus, a small e3t
-! OA at one place determines the number of subtimesteps globally
-! OA AWFULLY EXPENSIVE !! Not able to find a better approach. Damned !!
+      !
+      ! OA Below, this is garbage. the ideal would be to find a time-splitting 
+      ! OA algorithm that does not increase the computing cost by too much
+      ! OA In ROMS, I have included a time-splitting procedure. But it is 
+      ! OA too expensive as the loop is computed globally. Thus, a small e3t
+      ! OA at one place determines the number of subtimesteps globally
+      ! OA AWFULLY EXPENSIVE !! Not able to find a better approach. Damned !!
 
       DO jk = 1,jpkm1
@@ -522,6 +504,6 @@
       wscal(:,:,:) = wsbio4(:,:,:)
 
-!   INITIALIZE TO ZERO ALL THE SINKING ARRAYS
-!   -----------------------------------------
+      !  Initializa to zero all the sinking arrays 
+      !   -----------------------------------------
 
       sinking (:,:,:) = 0.e0
@@ -532,7 +514,6 @@
       sinkfer2(:,:,:) = 0.e0
 
-!   Compute the sedimentation term using p4zsink2 for all
-!   the sinking particles
-!   -----------------------------------------------------
+      !   Compute the sedimentation term using p4zsink2 for all the sinking particles
+      !   -----------------------------------------------------
 
       CALL p4z_sink2( wsbio3, sinking , jppoc )
@@ -543,45 +524,33 @@
       CALL p4z_sink2( wscal , sinkcal , jpcal )
 
-!  Exchange between organic matter compartments due to
-!  coagulation/disaggregation
-!  ---------------------------------------------------
+      !  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)
+# if defined key_degrad
+               zstep = xstep * facvol(ji,jj,jk)
+# else
+               zstep = xstep 
+# endif
+               zfact = zstep * xdiss(ji,jj,jk)
                !  Part I : Coagulation dependent on turbulence
-# if defined key_degrad
-               zagg1 = 940.* zfact * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jppoc) * facvol(ji,jj,jk)
-               zagg2 = 1.054e4 * zfact * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc) * facvol(ji,jj,jk)
-# else
                zagg1 = 940.* zfact * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jppoc)
                zagg2 = 1.054e4 * zfact * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc)
-# endif
 
                ! Part II : Differential settling
 
                !  Aggregation of small into large particles
-# if defined key_degrad
-               zagg3 = 0.66 * xstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc) * facvol(ji,jj,jk)
-               zagg4 = 0.e0 * xstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jppoc) * facvol(ji,jj,jk)
-# else
-               zagg3 = 0.66 * xstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc)
-               zagg4 = 0.e0 * xstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jppoc)
-# endif
+               zagg3 = 0.66 * zstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc)
+               zagg4 = 0.e0 * zstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jppoc)
+
                zagg   = zagg1 + zagg2 + zagg3 + zagg4
                zaggfe = zagg * trn(ji,jj,jk,jpsfe) / ( trn(ji,jj,jk,jppoc) + rtrn )
 
                ! Aggregation of DOC to small particles
-#if defined key_degrad
-               zaggdoc = ( 80.* trn(ji,jj,jk,jpdoc) + 698. * trn(ji,jj,jk,jppoc) )       &
-                  &      * facvol(ji,jj,jk)  * zfact * trn(ji,jj,jk,jpdoc)
-               zaggdoc2 = 1.05e4 * zfact * trn(ji,jj,jk,jpgoc)   &
-                  &      * facvol(ji,jj,jk) * trn(ji,jj,jk,jpdoc)
-#else
-               zaggdoc = ( 80.* trn(ji,jj,jk,jpdoc) + 698. * trn(ji,jj,jk,jppoc) )    &
-                  &      *  zfact * trn(ji,jj,jk,jpdoc)
+               zaggdoc = ( 80.* trn(ji,jj,jk,jpdoc) + 698. * trn(ji,jj,jk,jppoc) ) *  zfact * trn(ji,jj,jk,jpdoc) 
                zaggdoc2 = 1.05e4 * zfact * trn(ji,jj,jk,jpgoc) * trn(ji,jj,jk,jpdoc)
-#endif
+
                !  Update the trends
                tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) - zagg + zaggdoc
@@ -623,4 +592,10 @@
    END SUBROUTINE p4z_sink
 
+   SUBROUTINE p4z_sink_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE p4z_sink_init  ***
+      !!----------------------------------------------------------------------
+   END SUBROUTINE p4z_sink_init
+
 #endif
 
@@ -641,9 +616,11 @@
       !!
       INTEGER  ::   ji, jj, jk, jn
-      REAL(wp) ::   zigma,zew,zign, zflx
+      REAL(wp) ::   zigma,zew,zign, zflx, zstep
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::  ztraz, zakz
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zwsink2
       !!---------------------------------------------------------------------
 
+
+      zstep = rfact2 / 2.
 
       ztraz(:,:,:) = 0.e0
@@ -693,7 +670,7 @@
             DO jj = 1, jpj      
                DO ji = 1, jpi    
-                  zigma = zwsink2(ji,jj,jk+1) * xstep2 / fse3w(ji,jj,jk+1)
+                  zigma = zwsink2(ji,jj,jk+1) * zstep / fse3w(ji,jj,jk+1)
                   zew   = zwsink2(ji,jj,jk+1)
-                  psinkflx(ji,jj,jk+1) = -zew * ( trn(ji,jj,jk,jp_tra) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * xstep2
+                  psinkflx(ji,jj,jk+1) = -zew * ( trn(ji,jj,jk,jp_tra) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * zstep
                END DO
             END DO
