New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 10375 for NEMO/trunk/src/TOP/PISCES/P4Z/p4zsink.F90 – NEMO

Ignore:
Timestamp:
2018-12-07T17:27:26+01:00 (5 years ago)
Author:
aumont
Message:

Create generic routines for vertical sinking of particles

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/TOP/PISCES/P4Z/p4zsink.F90

    r10362 r10375  
    1616   USE trc             !  passive tracers common variables  
    1717   USE sms_pisces      !  PISCES Source Minus Sink variables 
     18   USE trcsink         !  General routine to compute sedimentation 
    1819   USE prtctl_trc      !  print control for debugging 
    1920   USE iom             !  I/O manager 
     
    5960      !!--------------------------------------------------------------------- 
    6061      INTEGER, INTENT(in) :: kt, knt 
    61       INTEGER  ::   ji, jj, jk, jit 
    62       INTEGER  ::   iiter1, iiter2 
    63       REAL(wp) ::   zagg1, zagg2, zagg3, zagg4 
    64       REAL(wp) ::   zagg , zaggfe, zaggdoc, zaggdoc2, zaggdoc3 
    65       REAL(wp) ::   zfact, zwsmax, zmax 
     62      INTEGER  ::   ji, jj, jk 
    6663      CHARACTER (len=25) :: charout 
     64      REAL(wp) :: zmax, zfact 
    6765      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d 
    6866      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) :: zw2d 
     
    7068      ! 
    7169      IF( ln_timing )   CALL timing_start('p4z_sink') 
    72  
    7370 
    7471      ! Initialization of some global variables 
     
    9794 
    9895      ! 
    99       ! OA This is (I hope) a temporary solution for the problem that may  
    100       ! OA arise in specific situation where the CFL criterion is broken  
    101       ! OA for vertical sedimentation of particles. To avoid this, a time 
    102       ! OA splitting algorithm has been coded. A specific maximum 
    103       ! OA iteration number is provided and may be specified in the namelist  
    104       ! OA This is to avoid very large iteration number when explicit free 
    105       ! OA surface is used (for instance). When niter?max is set to 1,  
    106       ! OA this computation is skipped. The crude old threshold method is  
    107       ! OA then applied. This also happens when niter exceeds nitermax. 
    108       IF( MAX( niter1max, niter2max ) == 1 ) THEN 
    109         iiter1 = 1 
    110         iiter2 = 1 
    111       ELSE 
    112         iiter1 = 1 
    113         iiter2 = 1 
    114         DO jk = 1, jpkm1 
    115           DO jj = 1, jpj 
    116              DO ji = 1, jpi 
    117                 IF( tmask(ji,jj,jk) == 1) THEN 
    118                    zwsmax =  0.5 * e3t_n(ji,jj,jk) / xstep 
    119                    iiter1 =  MAX( iiter1, INT( wsbio3(ji,jj,jk) / zwsmax ) ) 
    120                    iiter2 =  MAX( iiter2, INT( wsbio4(ji,jj,jk) / zwsmax ) ) 
    121                 ENDIF 
    122              END DO 
    123           END DO 
    124         END DO 
    125         IF( lk_mpp ) THEN 
    126            CALL mpp_max( iiter1 ) 
    127            CALL mpp_max( iiter2 ) 
    128         ENDIF 
    129         iiter1 = MIN( iiter1, niter1max ) 
    130         iiter2 = MIN( iiter2, niter2max ) 
    131       ENDIF 
    132  
    133       DO jk = 1,jpkm1 
    134          DO jj = 1, jpj 
    135             DO ji = 1, jpi 
    136                IF( tmask(ji,jj,jk) == 1 ) THEN 
    137                  zwsmax = 0.5 * e3t_n(ji,jj,jk) / xstep 
    138                  wsbio3(ji,jj,jk) = MIN( wsbio3(ji,jj,jk), zwsmax * REAL( iiter1, wp ) ) 
    139                  wsbio4(ji,jj,jk) = MIN( wsbio4(ji,jj,jk), zwsmax * REAL( iiter2, wp ) ) 
    140                ENDIF 
    141             END DO 
    142          END DO 
    143       END DO 
    144  
    14596      !  Initializa to zero all the sinking arrays  
    14697      !   ----------------------------------------- 
     
    154105      !   Compute the sedimentation term using p4zsink2 for all the sinking particles 
    155106      !   ----------------------------------------------------- 
    156       DO jit = 1, iiter1 
    157         CALL p4z_sink2( wsbio3, sinking , jppoc, iiter1 ) 
    158         CALL p4z_sink2( wsbio3, sinkfer , jpsfe, iiter1 ) 
    159       END DO 
    160  
    161       DO jit = 1, iiter2 
    162         CALL p4z_sink2( wsbio4, sinking2, jpgoc, iiter2 ) 
    163         CALL p4z_sink2( wsbio4, sinkfer2, jpbfe, iiter2 ) 
    164         CALL p4z_sink2( wsbio4, sinksil , jpgsi, iiter2 ) 
    165         CALL p4z_sink2( wsbio4, sinkcal , jpcal, iiter2 ) 
    166       END DO 
     107      CALL trc_sink( kt, wsbio3, sinking , jppoc, rfact2 ) 
     108      CALL trc_sink( kt, wsbio3, sinkfer , jpsfe, rfact2 ) 
     109      CALL trc_sink( kt, wsbio4, sinking2, jpgoc, rfact2 ) 
     110      CALL trc_sink( kt, wsbio4, sinkfer2, jpbfe, rfact2 ) 
     111      CALL trc_sink( kt, wsbio4, sinksil , jpgsi, rfact2 ) 
     112      CALL trc_sink( kt, wsbio4, sinkcal , jpcal, rfact2 ) 
    167113 
    168114      IF( ln_p5z ) THEN 
     
    174120         !   Compute the sedimentation term using p4zsink2 for all the sinking particles 
    175121         !   ----------------------------------------------------- 
    176          DO jit = 1, iiter1 
    177            CALL p4z_sink2( wsbio3, sinkingn , jppon, iiter1 ) 
    178            CALL p4z_sink2( wsbio3, sinkingp , jppop, iiter1 ) 
    179          END DO 
    180  
    181          DO jit = 1, iiter2 
    182            CALL p4z_sink2( wsbio4, sinking2n, jpgon, iiter2 ) 
    183            CALL p4z_sink2( wsbio4, sinking2p, jpgop, iiter2 ) 
    184          END DO 
     122         CALL trc_sink( kt, wsbio3, sinkingn , jppon, rfact2 ) 
     123         CALL trc_sink( kt, wsbio3, sinkingp , jppop, rfact2 ) 
     124         CALL trc_sink( kt, wsbio4, sinking2n, jpgon, rfact2 ) 
     125         CALL trc_sink( kt, wsbio4, sinking2p, jpgop, rfact2 ) 
    185126      ENDIF 
    186127 
    187128      IF( ln_ligand ) THEN 
    188129         wsfep (:,:,:) = wfep 
    189          DO jk = 1,jpkm1 
    190             DO jj = 1, jpj 
    191                DO ji = 1, jpi 
    192                   IF( tmask(ji,jj,jk) == 1 ) THEN 
    193                     zwsmax = 0.5 * e3t_n(ji,jj,jk) / xstep 
    194                     wsfep(ji,jj,jk) = MIN( wsfep(ji,jj,jk), zwsmax * REAL( iiter1, wp ) ) 
    195                   ENDIF 
    196                END DO 
    197             END DO 
    198          END DO 
    199130         ! 
    200131         sinkfep(:,:,:) = 0.e0 
    201          DO jit = 1, iiter1 
    202            CALL p4z_sink2( wsfep, sinkfep , jpfep, iiter1 ) 
    203          END DO 
     132         CALL trc_sink( kt, wsfep, sinkfep , jpfep, rfact2 ) 
    204133      ENDIF 
    205134 
     
    281210   END SUBROUTINE p4z_sink_init 
    282211 
    283  
    284    SUBROUTINE p4z_sink2( pwsink, psinkflx, jp_tra, kiter ) 
    285       !!--------------------------------------------------------------------- 
    286       !!                     ***  ROUTINE p4z_sink2  *** 
    287       !! 
    288       !! ** Purpose :   Compute the sedimentation terms for the various sinking 
    289       !!     particles. The scheme used to compute the trends is based 
    290       !!     on MUSCL. 
    291       !! 
    292       !! ** Method  : - this ROUTINE compute not exactly the advection but the 
    293       !!      transport term, i.e.  div(u*tra). 
    294       !!--------------------------------------------------------------------- 
    295       INTEGER , INTENT(in   )                         ::   jp_tra    ! tracer index index       
    296       INTEGER , INTENT(in   )                         ::   kiter     ! number of iterations for time-splitting  
    297       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pwsink    ! sinking speed 
    298       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   psinkflx  ! sinking fluxe 
    299       ! 
    300       INTEGER  ::   ji, jj, jk, jn 
    301       REAL(wp) ::   zigma,zew,zign, zflx, zstep 
    302       REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztraz, zakz, zwsink2, ztrb  
    303       !!--------------------------------------------------------------------- 
    304       ! 
    305       IF( ln_timing )   CALL timing_start('p4z_sink2') 
    306       ! 
    307       zstep = rfact2 / REAL( kiter, wp ) / 2. 
    308  
    309       ztraz(:,:,:) = 0.e0 
    310       zakz (:,:,:) = 0.e0 
    311       ztrb (:,:,:) = trb(:,:,:,jp_tra) 
    312  
    313       DO jk = 1, jpkm1 
    314          zwsink2(:,:,jk+1) = -pwsink(:,:,jk) / rday * tmask(:,:,jk+1)  
    315       END DO 
    316       zwsink2(:,:,1) = 0.e0 
    317  
    318  
    319       ! Vertical advective flux 
    320       DO jn = 1, 2 
    321          !  first guess of the slopes interior values 
    322          DO jk = 2, jpkm1 
    323             ztraz(:,:,jk) = ( trb(:,:,jk-1,jp_tra) - trb(:,:,jk,jp_tra) ) * tmask(:,:,jk) 
    324          END DO 
    325          ztraz(:,:,1  ) = 0.0 
    326          ztraz(:,:,jpk) = 0.0 
    327  
    328          ! slopes 
    329          DO jk = 2, jpkm1 
    330             DO jj = 1,jpj 
    331                DO ji = 1, jpi 
    332                   zign = 0.25 + SIGN( 0.25, ztraz(ji,jj,jk) * ztraz(ji,jj,jk+1) ) 
    333                   zakz(ji,jj,jk) = ( ztraz(ji,jj,jk) + ztraz(ji,jj,jk+1) ) * zign 
    334                END DO 
    335             END DO 
    336          END DO 
    337           
    338          ! Slopes limitation 
    339          DO jk = 2, jpkm1 
    340             DO jj = 1, jpj 
    341                DO ji = 1, jpi 
    342                   zakz(ji,jj,jk) = SIGN( 1., zakz(ji,jj,jk) ) *        & 
    343                      &             MIN( ABS( zakz(ji,jj,jk) ), 2. * ABS(ztraz(ji,jj,jk+1)), 2. * ABS(ztraz(ji,jj,jk) ) ) 
    344                END DO 
    345             END DO 
    346          END DO 
    347           
    348          ! vertical advective flux 
    349          DO jk = 1, jpkm1 
    350             DO jj = 1, jpj       
    351                DO ji = 1, jpi     
    352                   zigma = zwsink2(ji,jj,jk+1) * zstep / e3w_n(ji,jj,jk+1) 
    353                   zew   = zwsink2(ji,jj,jk+1) 
    354                   psinkflx(ji,jj,jk+1) = -zew * ( trb(ji,jj,jk,jp_tra) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * zstep 
    355                END DO 
    356             END DO 
    357          END DO 
    358          ! 
    359          ! Boundary conditions 
    360          psinkflx(:,:,1  ) = 0.e0 
    361          psinkflx(:,:,jpk) = 0.e0 
    362           
    363          DO jk=1,jpkm1 
    364             DO jj = 1,jpj 
    365                DO ji = 1, jpi 
    366                   zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / e3t_n(ji,jj,jk) 
    367                   trb(ji,jj,jk,jp_tra) = trb(ji,jj,jk,jp_tra) + zflx 
    368                END DO 
    369             END DO 
    370          END DO 
    371  
    372       ENDDO 
    373  
    374       DO jk = 1,jpkm1 
    375          DO jj = 1,jpj 
    376             DO ji = 1, jpi 
    377                zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / e3t_n(ji,jj,jk) 
    378                ztrb(ji,jj,jk) = ztrb(ji,jj,jk) + 2. * zflx 
    379             END DO 
    380          END DO 
    381       END DO 
    382  
    383       trb(:,:,:,jp_tra) = ztrb(:,:,:) 
    384       psinkflx(:,:,:)   = 2. * psinkflx(:,:,:) 
    385       ! 
    386       IF( ln_timing )  CALL timing_stop('p4z_sink2') 
    387       ! 
    388    END SUBROUTINE p4z_sink2 
    389  
    390  
    391212   INTEGER FUNCTION p4z_sink_alloc() 
    392213      !!---------------------------------------------------------------------- 
Note: See TracChangeset for help on using the changeset viewer.