MODULE sed_fluxes_mod ! !svn $Id: sed_fluxes.F 396 2009-09-11 18:53:38Z arango $ !==================================================== John C. Warner === ! Copyright (c) 2002-2009 The ROMS/TOMS Group Hernan G. Arango ! ! Licensed under a MIT/X style license ! ! See License_ROMS.txt ! !======================================================================= ! ! ! This computes sediment bed and water column exchanges: deposition, ! ! resuspension, and erosion. ! ! ! ! References: ! ! ! ! Warner, J.C., C.R. Sherwood, R.P. Signell, C.K. Harris, and H.G. ! ! Arango, 2008: Development of a three-dimensional, regional, ! ! coupled wave, current, and sediment-transport model, Computers ! ! & Geosciences, 34, 1284-1306. ! ! ! !======================================================================= ! implicit none PRIVATE PUBLIC :: sed_fluxes CONTAINS ! !*********************************************************************** SUBROUTINE sed_fluxes (ng, tile) !*********************************************************************** ! USE mod_param USE mod_forces USE mod_grid USE mod_ocean USE mod_stepping USE mod_bbl ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile ! ! Local variable declarations. ! integer :: IminS, ImaxS, JminS, JmaxS integer :: LBi, UBi, LBj, UBj, LBij, UBij ! ! Set horizontal starting and ending indices for automatic private storage ! arrays. ! IminS=BOUNDS(ng)%Istr(tile)-3 ImaxS=BOUNDS(ng)%Iend(tile)+3 JminS=BOUNDS(ng)%Jstr(tile)-3 JmaxS=BOUNDS(ng)%Jend(tile)+3 ! ! Determine array lower and upper bounds in the I- and J-directions. ! LBi=BOUNDS(ng)%LBi(tile) UBi=BOUNDS(ng)%UBi(tile) LBj=BOUNDS(ng)%LBj(tile) UBj=BOUNDS(ng)%UBj(tile) ! ! Set array lower and upper bounds for MIN(I,J)- and MAX(I,J)-directions. ! LBij=BOUNDS(ng)%LBij UBij=BOUNDS(ng)%UBij ! CALL wclock_on (ng, iNLM, 16) CALL sed_fluxes_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & nstp(ng), nnew(ng), & & GRID(ng) % Hz, & & GRID(ng) % rmask_wet, & & BBL(ng) % bustrc, & & BBL(ng) % bvstrc, & & BBL(ng) % bustrw, & & BBL(ng) % bvstrw, & & BBL(ng) % bustrcwmax, & & BBL(ng) % bvstrcwmax, & & FORCES(ng) % bustr, & & FORCES(ng) % bvstr, & & OCEAN(ng) % t, & & OCEAN(ng) % ero_flux, & & OCEAN(ng) % settling_flux, & & OCEAN(ng) % bed, & & OCEAN(ng) % bed_frac, & & OCEAN(ng) % bed_mass, & & OCEAN(ng) % bottom) CALL wclock_off (ng, iNLM, 16) RETURN END SUBROUTINE sed_fluxes ! !*********************************************************************** SUBROUTINE sed_fluxes_tile (ng, tile, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & nstp, nnew, & & Hz, & & rmask_wet, & & bustrc, bvstrc, & & bustrw, bvstrw, & & bustrcwmax, bvstrcwmax, & & bustr, bvstr, & & t, & & ero_flux, settling_flux, & & bed, bed_frac, bed_mass, & & bottom) !*********************************************************************** ! USE mod_param USE mod_scalars USE mod_sediment ! ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile integer, intent(in) :: LBi, UBi, LBj, UBj integer, intent(in) :: IminS, ImaxS, JminS, JmaxS integer, intent(in) :: nstp, nnew ! real(r8), intent(in) :: Hz(LBi:,LBj:,:) real(r8), intent(in) :: rmask_wet(LBi:,LBj:) real(r8), intent(in) :: bustrc(LBi:,LBj:) real(r8), intent(in) :: bvstrc(LBi:,LBj:) real(r8), intent(in) :: bustrw(LBi:,LBj:) real(r8), intent(in) :: bvstrw(LBi:,LBj:) real(r8), intent(in) :: bustrcwmax(LBi:,LBj:) real(r8), intent(in) :: bvstrcwmax(LBi:,LBj:) real(r8), intent(in) :: bustr(LBi:,LBj:) real(r8), intent(in) :: bvstr(LBi:,LBj:) real(r8), intent(inout) :: t(LBi:,LBj:,:,:,:) real(r8), intent(inout) :: ero_flux(LBi:,LBj:,:) real(r8), intent(inout) :: settling_flux(LBi:,LBj:,:) real(r8), intent(inout) :: bed(LBi:,LBj:,:,:) real(r8), intent(inout) :: bed_frac(LBi:,LBj:,:,:) real(r8), intent(inout) :: bed_mass(LBi:,LBj:,:,:,:) real(r8), intent(inout) :: bottom(LBi:,LBj:,:) ! ! Local variable declarations. ! logical :: EWperiodic=.FALSE. logical :: NSperiodic=.FALSE. integer :: Ksed, i, indx, ised, j, k, ks integer :: bnew real(r8) :: cff, cff1, cff2, cff3, cff4 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tau_w ! !----------------------------------------------------------------------- ! Set lower and upper tile bounds and staggered variables bounds for ! this horizontal domain partition. Notice that if tile=-1, it will ! set the values for the global grid. !----------------------------------------------------------------------- ! integer :: Istr, IstrR, IstrT, IstrU, Iend, IendR, IendT integer :: Jstr, JstrR, JstrT, JstrV, Jend, JendR, JendT ! Istr =BOUNDS(ng)%Istr (tile) IstrR=BOUNDS(ng)%IstrR(tile) IstrT=BOUNDS(ng)%IstrT(tile) IstrU=BOUNDS(ng)%IstrU(tile) Iend =BOUNDS(ng)%Iend (tile) IendR=BOUNDS(ng)%IendR(tile) IendT=BOUNDS(ng)%IendT(tile) Jstr =BOUNDS(ng)%Jstr (tile) JstrR=BOUNDS(ng)%JstrR(tile) JstrT=BOUNDS(ng)%JstrT(tile) JstrV=BOUNDS(ng)%JstrV(tile) Jend =BOUNDS(ng)%Jend (tile) JendR=BOUNDS(ng)%JendR(tile) JendT=BOUNDS(ng)%JendT(tile) bnew=nstp ! !----------------------------------------------------------------------- ! Compute sediment deposition, resuspension, and erosion. !----------------------------------------------------------------------- ! DO j=Jstr-1,Jend+1 DO i=Istr-1,Iend+1 tau_w(i,j)=SQRT(bustrcwmax(i,j)*bustrcwmax(i,j)+ & & bvstrcwmax(i,j)*bvstrcwmax(i,j)) tau_w(i,j)=tau_w(i,j)*rmask_wet(i,j) END DO END DO ! !----------------------------------------------------------------------- ! Sediment deposition and resuspension near the bottom. !----------------------------------------------------------------------- ! ! The deposition and resuspension of sediment on the bottom "bed" ! is due to precepitation settling_flux, already computed, and the ! resuspension (erosion, hence called ero_flux). The resuspension is ! applied to the bottom-most grid box value qc(:,1) so the total mass ! is conserved. Restrict "ero_flux" so that "bed" cannot go negative ! after both fluxes are applied. ! J_LOOP : DO j=Jstr,Jend DO k=1,N(ng) DO i=Istr,Iend Hz_inv(i,k)=1.0_r8/Hz(i,j,k) END DO END DO ! SED_LOOP: DO ised=1,NST indx=idsed(ised) DO i=Istr,Iend ! ! Calculate critical shear stress in Pa ! cff=1.0_r8/tau_ce(ised,ng) ! ! Compute erosion, ero_flux (kg/m2). ! cff1=(1.0_r8-bed(i,j,1,iporo))*bed_frac(i,j,1,ised) cff2=dt(ng)*Erate(ised,ng)*cff1 cff3=Srho(ised,ng)*cff1 cff4=bed_mass(i,j,1,bnew,ised) ero_flux(i,j,ised)= & & MIN(MAX(0.0_r8,cff2*(cff*tau_w(i,j)-1.0_r8)), & & MIN(cff3*bottom(i,j,iactv),cff4)+ & & settling_flux(i,j,ised)) ! ! Update global tracer variables (m Tunits for nnew indx, Tuints for 3) ! for erosive flux. ! t(i,j,1,nnew,indx)=t(i,j,1,nnew,indx)+ero_flux(i,j,ised) t(i,j,1,3,indx)=t(i,j,1,3,indx)+ & & ero_flux(i,j,ised)*Hz_inv(i,1) END DO END DO SED_LOOP END DO J_LOOP RETURN END SUBROUTINE sed_fluxes_tile END MODULE sed_fluxes_mod