source: trunk/NEMO/TOP_SRC/SMS/p4zsink.kriest.h @ 700

CC----------------------------------------------------------------------
CC  TOP 1.0 , LOCEAN-IPSL (2005) 
CC $Id$
CC This software is governed by CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
CC----------------------------------------------------------------------

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