source: trunk/NEMO/TOP_SRC/PISCES/p4zsink2.F90 @ 935

MODULE p4zsink2
   !!======================================================================
   !!                         ***  MODULE p4zsink2  ***
   !! TOP :   PISCES Compute 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
   !!----------------------------------------------------------------------
#if defined key_pisces
   !!----------------------------------------------------------------------
   !!   'key_pisces'                                       PISCES bio-model
   !!----------------------------------------------------------------------
   !!   p4z_sink2       :  Compute vertical flux of particulate matter due to gravitational sinking
   !!----------------------------------------------------------------------
   USE oce_trc         !
   USE trp_trc
   USE sms

   IMPLICIT NONE
   PRIVATE

   PUBLIC   p4z_sink2    ! called in p4zbio.F90

   !!* Substitution
#  include "domzgr_substitute.h90"
   !!----------------------------------------------------------------------
   !! NEMO/TOP 2.0 , LOCEAN-IPSL (2007) 
   !! $Header:$ 
   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------

CONTAINS

   SUBROUTINE p4z_sink2( wstmp, sinktemp, jn )
      !!---------------------------------------------------------------------
      !!                     ***  ROUTINE p4z_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   )                         ::   jn         ! tracer index index      
      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   wstmp      ! ???
      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   sinktemp   ! ???
      !!
      INTEGER  ::   ji, jj, jk, jnt
      REAL(wp) ::   zigma,zew,zstep,zign
      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ztraz, zakz
      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zkz  , zwstmp2
      !!---------------------------------------------------------------------

      zstep  = rfact2 / 2.

      ztraz(:,:,:) = 0.e0
      zkz  (:,:,:) = 0.e0
      zakz (:,:,:) = 0.e0

      DO jk = 1, jpkm1
# if defined key_off_degrad
         zwstmp2(:,:,jk+1)=-wstmp(:,:,jk)/rjjss*tmask(:,:,jk+1)*facvol(:,:,jk)
# else
         zwstmp2(:,:,jk+1)=-wstmp(:,:,jk)/rjjss*tmask(:,:,jk+1)

# endif
      END DO
 
      zwstmp2(:,:,1) = 0.e0
!
! Vertical advective flux
!-------------------------------

      DO jnt = 1, 2

! ... first guess of the slopes interior values

         DO jk = 2, jpkm1
            ztraz(:,:,jk) = (trn(:,:,jk-1,jn) - trn(:,:,jk,jn)) *tmask(:,:,jk)
         END DO

         ztraz(:,:,1  ) = 0.0
         ztraz(:,:,jpk) = 0.0
!
! slopes
         DO jk=2,jpkm1
            DO jj = 1,jpj
               DO ji = 1, jpi
                  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
            END DO
         END DO
         !
         ! Slopes limitation
         DO jk = 2, jpkm1
            DO jj = 1, jpj
               DO ji = 1, jpi
                  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
            END DO
         END DO
         
         ! vertical advective flux
         DO jk = 1, jpkm1
            DO jj = 1, jpj      
               DO ji = 1, jpi    
                  zigma = zwstmp2(ji,jj,jk+1)*zstep/fse3w(ji,jj,jk+1)
                  zew   = zwstmp2(ji,jj,jk+1)
                  sinktemp(ji,jj,jk+1) = -zew*(trn(ji,jj,jk,jn)                &
                     &                 -0.5*(1+zigma)*zakz(ji,jj,jk))*zstep
               END DO
            END DO
         END DO
         !
         ! Boundary conditions
         sinktemp(:,:,1  ) = 0.e0
         sinktemp(:,:,jpk) = 0.e0
         
         DO jk=1,jpkm1
            DO jj = 1,jpj
               DO ji = 1, jpi
                  trn(ji,jj,jk,jn) = trn(ji,jj,jk,jn)               &
                     &        + (sinktemp(ji,jj,jk)-sinktemp(ji,jj,jk+1))     &
                     &        /fse3t(ji,jj,jk)
               END DO
            END DO
         END DO

      ENDDO

      DO jk=1,jpkm1
         DO jj = 1,jpj
            DO ji = 1, jpi
               ! 
               trb(ji,jj,jk,jn) = trb(ji,jj,jk,jn)      &
               &        + 2.*(sinktemp(ji,jj,jk)-sinktemp(ji,jj,jk+1)) &
               &        /fse3t(ji,jj,jk)
!
          ENDDO
        ENDDO
      ENDDO
!
        trn(:,:,:,jn)   = trb(:,:,:,jn)
        sinktemp(:,:,:) = 2. * sinktemp(:,:,:)

      !
   END SUBROUTINE p4z_sink2

#else
   !!======================================================================
   !!  Dummy module :                                   No PISCES bio-model
   !!======================================================================
CONTAINS
   SUBROUTINE p4z_sink2( wstmp, sinktemp, jn )         ! Empty routine
      INTEGER, INTENT( in ) ::   jn
      REAL   , INTENT( in ) ::   wstmp,sinktemp    
      WRITE(*,*) 'p4z_sink2: You should not have seen this print! error?', jn, wstmp, sinktemp
   END SUBROUTINE p4z_sink2
#endif 

   !!======================================================================
END MODULE  p4zsink2