Index: /NEMO/trunk/src/TOP/TRP/trcsink.F90
===================================================================
--- /NEMO/trunk/src/TOP/TRP/trcsink.F90	(revision 10376)
+++ /NEMO/trunk/src/TOP/TRP/trcsink.F90	(revision 10376)
@@ -0,0 +1,235 @@
+MODULE trcsink
+   !!======================================================================
+   !!                         ***  MODULE trcsink  ***
+   !! TOP :  vertical flux of particulate matter due to gravitational sinking
+   !!======================================================================
+   !! History :   1.0  !  2004     (O. Aumont) Original code
+   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
+   !!             3.4  !  2011-06  (O. Aumont, C. Ethe) Change aggregation formula
+   !!             3.5  !  2012-07  (O. Aumont) Introduce potential time-splitting
+   !!             4.0  !  2018-12  (O. Aumont) Generalize the PISCES code to make it usable by any model
+   !!----------------------------------------------------------------------
+   !!   trc_sink       :  Compute vertical flux of particulate matter due to gravitational sinking
+   !!----------------------------------------------------------------------
+   USE oce_trc         !  shared variables between ocean and passive tracers
+   USE trc             !  passive tracers common variables 
+   USE lib_mpp
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC trc_sink
+   PUBLIC trc_sink_ini
+
+   INTEGER, PUBLIC :: nitermax      !: Maximum number of iterations for sinking
+
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
+   !! $Id: trcsink.F90 10069 2018-08-28 14:12:24Z nicolasmartin $ 
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   !!----------------------------------------------------------------------
+   !!   'standard sinking parameterisation'                  ???
+   !!----------------------------------------------------------------------
+
+   SUBROUTINE trc_sink ( kt, pwsink, psinkflx, jp_tra, rsfact )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE trc_sink  ***
+      !!
+      !! ** Purpose :   Compute vertical flux of particulate matter due to 
+      !!                gravitational sinking
+      !!
+      !! ** Method  : - ???
+      !!---------------------------------------------------------------------
+      INTEGER , INTENT(in)  :: kt
+      INTEGER , INTENT(in)  :: jp_tra    ! tracer index index      
+      REAL(wp), INTENT(in)  :: rsfact    ! time step duration
+      REAL(wp), INTENT(in)   , DIMENSION(jpi,jpj,jpk) :: pwsink
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: psinkflx
+      INTEGER  ::   ji, jj, jk
+      INTEGER, DIMENSION(jpi, jpj) ::   iiter
+      REAL(wp) ::   zfact, zwsmax, zmax
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwsink
+      !!---------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('trc_sink')
+      !
+      !
+      ! OA This is (I hope) a temporary solution for the problem that may 
+      ! OA arise in specific situation where the CFL criterion is broken 
+      ! OA for vertical sedimentation of particles. To avoid this, a time
+      ! OA splitting algorithm has been coded. A specific maximum
+      ! OA iteration number is provided and may be specified in the namelist 
+      ! OA This is to avoid very large iteration number when explicit free
+      ! OA surface is used (for instance). When niter?max is set to 1, 
+      ! OA this computation is skipped. The crude old threshold method is 
+      ! OA then applied. This also happens when niter exceeds nitermax.
+      IF( nitermax == 1 ) THEN
+         iiter(:,:) = 1
+      ELSE
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               iiter(ji,jj) = 1
+               DO jk = 1, jpkm1
+                  IF( tmask(ji,jj,jk) == 1.0 ) THEN
+                      zwsmax =  0.5 * e3t_n(ji,jj,jk) * rday / rsfact
+                      iiter(ji,jj) =  MAX( iiter(ji,jj), INT( pwsink(ji,jj,jk) / zwsmax ) )
+                  ENDIF
+               END DO
+            END DO
+         END DO
+         iiter(ji,jj) = MIN( iiter(ji,jj), nitermax )
+      ENDIF
+
+      DO jk = 1,jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( tmask(ji,jj,jk) == 1 ) THEN
+                 zwsmax = 0.5 * e3t_n(ji,jj,jk) * rday / rsfact
+                 zwsink(ji,jj,jk) = MIN( pwsink(ji,jj,jk), zwsmax * REAL( iiter(ji,jj), wp ) )
+               ENDIF
+            END DO
+         END DO
+      END DO
+
+      !  Initializa to zero all the sinking arrays 
+      !  -----------------------------------------
+      psinkflx(:,:,:) = 0.e0
+
+      !   Compute the sedimentation term using trc_sink2 for the considered sinking particle
+      !   -----------------------------------------------------
+      CALL trc_sink2( zwsink, psinkflx, jp_tra, iiter, rsfact )
+      !
+      IF( ln_timing )   CALL timing_stop('trc_sink')
+      !
+   END SUBROUTINE trc_sink
+
+   SUBROUTINE trc_sink2( pwsink, psinkflx, jp_tra, kiter, rsfact )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE trc_sink2  ***
+      !!
+      !! ** Purpose :   Compute the sedimentation terms for the various sinking
+      !!     particles. The scheme used to compute the trends is based
+      !!     on MUSCL.
+      !!
+      !! ** Method  : - this ROUTINE compute not exactly the advection but the
+      !!      transport term, i.e.  div(u*tra).
+      !!---------------------------------------------------------------------
+      INTEGER,  INTENT(in   )                         ::   jp_tra    ! tracer index index      
+      REAL(wp), INTENT(in   )                         ::   rsfact    ! duration of time step
+      INTEGER,  INTENT(in   ), DIMENSION(jpi,jpj)     ::   kiter     ! number of iterations for time-splitting 
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pwsink    ! sinking speed
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   psinkflx  ! sinking fluxe
+      !
+      INTEGER  ::   ji, jj, jk, jn
+      REAL(wp) ::   zigma,zew,zign, zflx, zstep
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztraz, zakz, zwsink2, ztrb 
+      !!---------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('trc_sink2')
+      !
+      ztraz(:,:,:) = 0.e0
+      zakz (:,:,:) = 0.e0
+      ztrb (:,:,:) = trb(:,:,:,jp_tra)
+
+      DO jk = 1, jpkm1
+         zwsink2(:,:,jk+1) = -pwsink(:,:,jk) / rday * tmask(:,:,jk+1) 
+      END DO
+      zwsink2(:,:,1) = 0.e0
+
+
+      ! Vertical advective flux
+      DO jn = 1, 2
+         !  first guess of the slopes interior values
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               !
+               zstep = rsfact / REAL( kiter(ji,jj), wp ) / 2.
+               !              
+               DO jk = 2, jpkm1
+                  ztraz(ji,jj,jk) = ( trb(ji,jj,jk-1,jp_tra) - trb(ji,jj,jk,jp_tra) ) * tmask(ji,jj,jk)
+               END DO
+               ztraz(ji,jj,1  ) = 0.0
+               ztraz(ji,jj,jpk) = 0.0
+
+               ! slopes
+               DO jk = 2, jpkm1
+                  zign = 0.25 + SIGN( 0.25, ztraz(ji,jj,jk) * ztraz(ji,jj,jk+1) )
+                  zakz(ji,jj,jk) = ( ztraz(ji,jj,jk) + ztraz(ji,jj,jk+1) ) * zign
+               END DO
+         
+               ! Slopes limitation
+               DO jk = 2, jpkm1
+                  zakz(ji,jj,jk) = SIGN( 1., zakz(ji,jj,jk) ) *        &
+                     &             MIN( ABS( zakz(ji,jj,jk) ), 2. * ABS(ztraz(ji,jj,jk+1)), 2. * ABS(ztraz(ji,jj,jk) ) )
+               END DO
+         
+               ! vertical advective flux
+               DO jk = 1, jpkm1
+                  zigma = zwsink2(ji,jj,jk+1) * zstep / e3w_n(ji,jj,jk+1)
+                  zew   = zwsink2(ji,jj,jk+1)
+                  psinkflx(ji,jj,jk+1) = -zew * ( trb(ji,jj,jk,jp_tra) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * zstep
+               END DO
+               !
+               ! Boundary conditions
+               psinkflx(ji,jj,1  ) = 0.e0
+               psinkflx(ji,jj,jpk) = 0.e0
+         
+               DO jk=1,jpkm1
+                  zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / e3t_n(ji,jj,jk)
+                  trb(ji,jj,jk,jp_tra) = trb(ji,jj,jk,jp_tra) + zflx
+               END DO
+            END DO
+         END DO
+      END DO
+
+      DO jk = 1,jpkm1
+         DO jj = 1,jpj
+            DO ji = 1, jpi
+               zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / e3t_n(ji,jj,jk)
+               ztrb(ji,jj,jk) = ztrb(ji,jj,jk) + 2. * zflx
+            END DO
+         END DO
+      END DO
+
+      trb(:,:,:,jp_tra) = ztrb(:,:,:)
+      psinkflx(:,:,:)   = 2. * psinkflx(:,:,:)
+      !
+      IF( ln_timing )  CALL timing_stop('trc_sink2')
+      !
+   END SUBROUTINE trc_sink2
+
+  SUBROUTINE trc_sink_ini
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE trc_sink_ini ***
+      !!
+      !! ** Purpose :   read  namelist options 
+      !!----------------------------------------------------------------------
+      INTEGER ::   ios   ! Local integer output status for namelist read
+      !!
+      NAMELIST/namtrc_snk/ nitermax
+      !!----------------------------------------------------------------------
+      !
+      REWIND( numnat_ref )              ! namtrc_rad in reference namelist 
+      READ  ( numnat_ref, namtrc_snk, IOSTAT = ios, ERR = 907)
+907   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtrc_snk in reference namelist', lwp )
+      REWIND( numnat_cfg )              ! namtrc_rad in configuration namelist 
+      READ  ( numnat_cfg, namtrc_snk, IOSTAT = ios, ERR = 908 )
+908   IF( ios > 0 )   CALL ctl_nam ( ios , 'namtrc_snk in configuration namelist', lwp )
+      IF(lwm) WRITE( numont, namtrc_snk )
+
+      IF(lwp) THEN                     !   ! Control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'trc_sink : Sedimentation of particles '
+         WRITE(numout,*) '~~~~~~~ '
+         WRITE(numout,*) '   Namelist namtrc_snk : sedimentation of particles'
+         WRITE(numout,*) '   Maximum number of iterations   nitermax = ', nitermax
+         WRITE(numout,*)
+      ENDIF
+
+   END SUBROUTINE trc_sink_ini
+
+   !!======================================================================
+END MODULE trcsink
