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 – NEMO

Changeset 10375


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

Create generic routines for vertical sinking of particles

Location:
NEMO/trunk
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • NEMO/trunk/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_top_cfg

    r9516 r10375  
    9090/ 
    9191!----------------------------------------------------------------------- 
     92&namtrc_snk      !  sedimentation of particles 
     93!----------------------------------------------------------------------- 
     94/ 
     95!----------------------------------------------------------------------- 
    9296&namtrc_dmp      !   passive tracer newtonian damping    
    9397!----------------------------------------------------------------------- 

  • NEMO/trunk/cfgs/ORCA2_OFF_PISCES/EXPREF/namelist_top_cfg

    r9516 r10375  
    8989/ 
    9090!----------------------------------------------------------------------- 
     91&namtrc_snk      !  sedimentation of particles 
     92!----------------------------------------------------------------------- 
     93/ 
     94!----------------------------------------------------------------------- 
    9195&namtrc_dmp      !   passive tracer newtonian damping    
    9296!----------------------------------------------------------------------- 

  • NEMO/trunk/cfgs/SHARED/namelist_pisces_ref

    r10363 r10375  
    5151   wsbio2max  =  50.      ! Big particles maximum sinking speed 
    5252   wsbio2scale =  5000.    ! Big particles length scale of sinking 
    53    niter1max  =  1        ! Maximum number of iterations for POC 
    54    niter2max  =  2        ! Maximum number of iterations for GOC 
    5553!                         !  ln_ligand enabled 
    5654   wfep       =  0.2      ! FeP sinking speed  
    5755   ldocp      =  1.E-4    ! Phyto ligand production per unit doc  
    5856   ldocz      =  1.E-4    ! Zoo ligand production per unit doc  
    59    lthet      =  0.5      ! Proportional loss of ligands due to Fe uptake  
     57   lthet      =  1.0      ! Proportional loss of ligands due to Fe uptake  
    6058!                         !  ln_p5z enabled 
    6159   no3rat3    =  0.182    ! N/C ratio in zooplankton 
    … …  
    310308   xlam1        =  0.005  ! scavenging rate of Iron 
    311309   xlamdust     =  150.0  ! Scavenging rate of dust 
    312    ligand       =  0.6E-9 ! Ligands concentration  
     310   ligand       =  0.7E-9 ! Ligands concentration  
    313311   kfep         =  0.01   ! Nanoparticle formation rate constant 
    314312/ 
    … …  
    322320   xsilab    =  0.5       ! Fraction of labile biogenic silica 
    323321   feratb    =  10.E-6    ! Fe/C quota in bacteria 
    324    xkferb    =  2.5E-10   ! Half-saturation constant for bacteria Fe/C 
     322   xkferb    =  3E-10     ! Half-saturation constant for bacteria Fe/C 
    325323!                         ! ln_p5z 
    326324   xremikc   =  0.25      ! remineralization rate of DOC 
    … …  
    372370   ln_ironsed  =  .true.   ! boolean for Fe input from sediments 
    373371   ln_ironice  =  .true.   ! boolean for Fe input from sea ice 
    374    ln_hydrofe  =  .false.  ! boolean for from hydrothermal vents 
     372   ln_hydrofe  =  .true.   ! boolean for from hydrothermal vents 
    375373   sedfeinput  =  2.e-9    ! Coastal release of Iron 
    376374   distcoast   =  5.e3     ! Distance off the coast for Iron from sediments 

  • NEMO/trunk/cfgs/SHARED/namelist_top_ref

    r9526 r10375  
    9595/ 
    9696!----------------------------------------------------------------------- 
     97&namtrc_snk      !  Sedimentation of particles 
     98!----------------------------------------------------------------------- 
     99   nitermax      =  2   !  number of iterations for sedimentation 
     100/ 
     101!----------------------------------------------------------------------- 
    97102&namtrc_dmp      !   passive tracer newtonian damping                   (ln_trcdmp=T) 
    98103!----------------------------------------------------------------------- 

  • 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      !!---------------------------------------------------------------------- 

  • NEMO/trunk/src/TOP/PISCES/P4Z/p4zsms.F90

    r10341 r10375  
    7474            CALL p4z_che                              ! initialize the chemical constants 
    7575            CALL ahini_for_at(hi)   !  set PH at kt=nit000 
    76             t_oce_co2_flx_cum = 0._wp 
    7776        ELSE 
    7877            CALL p4z_rst( nittrc000, 'READ' )  !* read or initialize all required fields 
    … …  
    189188      !! 
    190189      NAMELIST/nampisbio/ nrdttrc, wsbio, xkmort, ferat3, wsbio2, wsbio2max, wsbio2scale,    & 
    191          &                   niter1max, niter2max, wfep, ldocp, ldocz, lthet,  & 
     190         &                   wfep, ldocp, ldocz, lthet,  & 
    192191         &                   no3rat3, po4rat3 
    193192         ! 
    … …  
    223222         WRITE(numout,*) '      Big particles maximum sinking speed       wsbio2max   =', wsbio2max 
    224223         WRITE(numout,*) '      Big particles sinking speed length scale  wsbio2scale =', wsbio2scale 
    225          WRITE(numout,*) '      Maximum number of iterations for POC      niter1max   =', niter1max 
    226          WRITE(numout,*) '      Maximum number of iterations for GOC      niter2max   =', niter2max 
    227224         IF( ln_ligand ) THEN 
    228225            WRITE(numout,*) '      FeP sinking speed                              wfep   =', wfep 

  • NEMO/trunk/src/TOP/PISCES/sms_pisces.F90

    r10362 r10375  
    3636 
    3737   !!*  Biological parameters  
    38    INTEGER  ::   niter1max, niter2max   !: Maximum number of iterations for sinking 
    3938   REAL(wp) ::   rno3              !: ??? 
    4039   REAL(wp) ::   o2ut              !: ??? 

  • NEMO/trunk/src/TOP/trcini.F90

    r10222 r10375  
    205205      USE trcldf , ONLY:  trc_ldf_ini 
    206206      USE trcrad , ONLY:  trc_rad_ini 
     207      USE trcsink, ONLY:  trc_sink_ini 
    207208      ! 
    208209      INTEGER :: ierr 
    … …  
    214215                       !                          ! vertical diffusion: always implicit time stepping scheme 
    215216                       CALL  trc_rad_ini          ! positivity of passive tracers  
     217                       CALL  trc_sink_ini         ! Vertical sedimentation of particles 
    216218      ! 
    217219   END SUBROUTINE trc_ini_trp 
Note: See TracChangeset for help on using the changeset viewer.