sed_bed.F

Go to the documentation of this file.

#include "cppdefs.h"
 
      MODULE sed_bed_mod
 
#if defined NONLINEAR && defined SEDIMENT && !defined COHESIVE_BED
!
!git $Id$
!==================================================== John C. Warner ===
!  Copyright (c) 2002-2025 The ROMS Group           Hernan G. Arango   !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.md                                               !
!=======================================================================
!                                                                      !
!  This routine computes sediment bed layer stratigraphy.              !
!                                                                      !
!  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_bed
!
      CONTAINS
!
!***********************************************************************
      SUBROUTINE sed_bed (ng, tile)
!***********************************************************************
!
      USE mod_param
      USE mod_forces
      USE mod_grid
      USE mod_ocean
      USE mod_sedbed
      USE mod_stepping
# ifdef BBL_MODEL
      USE mod_bbl
# endif
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, tile
!
!  Local variable declarations.
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__
!
# include "tile.h"
!
# ifdef PROFILE
      CALL wclock_on (ng, inlm, 16, __line__, myfile)
# endif
      CALL sed_bed_tile (ng, tile,                                      &
     &                   lbi, ubi, lbj, ubj,                            &
     &                   imins, imaxs, jmins, jmaxs,                    &
     &                   nstp(ng), nnew(ng),                            &
# ifdef WET_DRY
     &                   grid(ng) % rmask_wet,                          &
# endif
# ifdef BBL_MODEL
     &                   bbl(ng) % bustrc,                              &
     &                   bbl(ng) % bvstrc,                              &
     &                   bbl(ng) % bustrw,                              &
     &                   bbl(ng) % bvstrw,                              &
     &                   bbl(ng) % bustrcwmax,                          &
     &                   bbl(ng) % bvstrcwmax,                          &
# else
     &                   forces(ng) % bustr,                            &
     &                   forces(ng) % bvstr,                            &
# endif
     &                   ocean(ng) % t,                                 &
# ifdef SUSPLOAD
     &                   sedbed(ng) % ero_flux,                         &
     &                   sedbed(ng) % settling_flux,                    &
# endif
# if defined SED_MORPH
     &                   sedbed(ng) % bed_thick,                        &
# endif
     &                   sedbed(ng) % bed,                              &
     &                   sedbed(ng) % bed_frac,                         &
     &                   sedbed(ng) % bed_mass,                         &
     &                   sedbed(ng) % bottom)
# ifdef PROFILE
      CALL wclock_off (ng, inlm, 16, __line__, myfile)
# endif
!
      RETURN
      END SUBROUTINE sed_bed
!
!***********************************************************************
      SUBROUTINE sed_bed_tile (ng, tile,                                &
     &                         LBi, UBi, LBj, UBj,                      &
     &                         IminS, ImaxS, JminS, JmaxS,              &
     &                         nstp, nnew,                              &
# ifdef WET_DRY
     &                         rmask_wet,                               &
# endif
# ifdef BBL_MODEL
     &                         bustrc, bvstrc,                          &
     &                         bustrw, bvstrw,                          &
     &                         bustrcwmax, bvstrcwmax,                  &
# else
     &                         bustr, bvstr,                            &
# endif
     &                         t,                                       &
# ifdef SUSPLOAD
     &                         ero_flux, settling_flux,                 &
# endif
# if defined SED_MORPH
     &                         bed_thick,                               &
# endif
     &                         bed, bed_frac, bed_mass,                 &
     &                         bottom)
!***********************************************************************
!
      USE mod_param
      USE mod_scalars
      USE mod_sediment
!
      USE bc_3d_mod, ONLY : bc_r3d_tile
      USE exchange_2d_mod, ONLY : exchange_r2d_tile
# ifdef DISTRIBUTE
      USE mp_exchange_mod, ONLY : mp_exchange3d, mp_exchange4d
# endif
!
!  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
!
# ifdef ASSUMED_SHAPE
#  ifdef WET_DRY
      real(r8), intent(in) :: rmask_wet(LBi:,LBj:)
#  endif
#  ifdef BBL_MODEL
      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:)
#  else
      real(r8), intent(in) :: bustr(LBi:,LBj:)
      real(r8), intent(in) :: bvstr(LBi:,LBj:)
#  endif
#  if defined SED_MORPH
      real(r8), intent(inout):: bed_thick(LBi:,LBj:,:)
#  endif
      real(r8), intent(inout) :: t(LBi:,LBj:,:,:,:)
#  ifdef SUSPLOAD
      real(r8), intent(inout) :: ero_flux(LBi:,LBj:,:)
      real(r8), intent(inout) :: settling_flux(LBi:,LBj:,:)
#  endif
      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:,:)
# else
#  ifdef WET_DRY
      real(r8), intent(in) :: rmask_wet(LBi:UBi,LBj:UBj)
#  endif
#  ifdef BBL_MODEL
      real(r8), intent(in) :: bustrc(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: bvstrc(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: bustrw(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: bvstrw(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: bustrcwmax(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: bvstrcwmax(LBi:UBi,LBj:UBj)
#  else
      real(r8), intent(in) :: bustr(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: bvstr(LBi:UBi,LBj:UBj)
#  endif
#  if defined SED_MORPH
      real(r8), intent(inout):: bed_thick(LBi:UBi,LBj:UBj,3)
#  endif
      real(r8), intent(inout) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
#  ifdef SUSPLOAD
      real(r8), intent(inout) :: ero_flux(LBi:UBi,LBj:UBj,NST)
      real(r8), intent(inout) :: settling_flux(LBi:UBi,LBj:UBj,NST)
#  endif
      real(r8), intent(inout) :: bed(LBi:UBi,LBj:UBj,Nbed,MBEDP)
      real(r8), intent(inout) :: bed_frac(LBi:UBi,LBj:UBj,Nbed,NST)
      real(r8), intent(inout) :: bed_mass(LBi:UBi,LBj:UBj,Nbed,1:2,NST)
      real(r8), intent(inout) :: bottom(LBi:UBi,LBj:UBj,MBOTP)
# endif
!
!  Local variable declarations.
!
      integer :: Ksed, i, ised, j, k, ks
      integer :: bnew
 
      real(r8), parameter :: eps = 1.0e-14_r8
 
      real(r8) :: cff, cff1, cff2, cff3
      real(r8) :: thck_avail, thck_to_add
 
      real(r8), dimension(IminS:ImaxS,NST) :: dep_mass
      real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tau_w
 
# include "set_bounds.h"
 
# ifdef BEDLOAD
      bnew=nnew
# else
      bnew=nstp
# endif
!
!-----------------------------------------------------------------------
! Compute sediment bed layer stratigraphy.
!-----------------------------------------------------------------------
!
# if defined BEDLOAD_MPM || defined SUSPLOAD
#  ifdef BBL_MODEL
      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))
#   ifdef WET_DRY
          tau_w(i,j)=tau_w(i,j)*rmask_wet(i,j)
#   endif
        END DO
      END DO
#  else
        DO j=jstrm1,jendp1
          DO i=istrm1,iendp1
          tau_w(i,j)=0.5_r8*sqrt((bustr(i,j)+bustr(i+1,j))*             &
     &                           (bustr(i,j)+bustr(i+1,j))+             &
     &                           (bvstr(i,j)+bvstr(i,j+1))*             &
     &                           (bvstr(i,j)+bvstr(i,j+1)))
#   ifdef WET_DRY
          tau_w(i,j)=tau_w(i,j)*rmask_wet(i,j)
#   endif
        END DO
      END DO
#  endif
# endif
!
!-----------------------------------------------------------------------
!  Update bed properties according to ero_flux and dep_flux.
!-----------------------------------------------------------------------
!
# ifdef SUSPLOAD
      j_loop : DO j=jstr,jend
        sed_loop: DO ised=1,nst
!
!  The deposition and resuspension of sediment on the bottom "bed"
!  is due to precepitation flux FC(:,0), 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.
!
          DO i=istr,iend
            dep_mass(i,ised)=0.0_r8
 
#  ifdef SED_MORPH
!
! Apply morphology factor.
!
            ero_flux(i,j,ised)=ero_flux(i,j,ised)*morph_fac(ised,ng)
            settling_flux(i,j,ised)=settling_flux(i,j,ised)*            &
     &                              morph_fac(ised,ng)
!
#  endif
            IF ((ero_flux(i,j,ised)-settling_flux(i,j,ised)).lt.        &
     &           0.0_r8) THEN
!
!  If first time step of deposit, then store deposit material in
!  temporary array, dep_mass.
!
              IF ((time(ng).gt.(bed(i,j,1,iaged)+1.1_r8*dt(ng))).and.   &
     &            (bed(i,j,1,ithck).gt.newlayer_thick(ng))) THEN
                dep_mass(i,ised)=settling_flux(i,j,ised)-               &
     &                           ero_flux(i,j,ised)
              END IF
                bed(i,j,1,iaged)=time(ng)
            END IF
!
!  Update bed mass arrays.
!
            bed_mass(i,j,1,nnew,ised)=max(bed_mass(i,j,1,bnew,ised)-    &
     &                                   (ero_flux(i,j,ised)-           &
     &                                    settling_flux(i,j,ised)),     &
     &                                    0.0_r8)
            DO k=2,nbed
              bed_mass(i,j,k,nnew,ised)=bed_mass(i,j,k,nstp,ised)
            END DO
          END DO
        END DO sed_loop
!
!  If first time step of deposit, create new layer and combine bottom
!  two bed layers.
!
        DO i=istr,iend
          cff=0.0_r8
!
!  Determine if deposition ocurred here.
!
          IF (nbed.gt.1) THEN
            DO ised=1,nst
              cff=cff+dep_mass(i,ised)
            END DO
            IF (cff.gt.0.0_r8) THEN
!
!  Combine bottom layers.
!
              bed(i,j,nbed,iporo)=0.5_r8*(bed(i,j,nbed-1,iporo)+        &
     &                                    bed(i,j,nbed,iporo))
              bed(i,j,nbed,iaged)=0.5_r8*(bed(i,j,nbed-1,iaged)+        &
     &                                    bed(i,j,nbed,iaged))
              DO ised=1,nst
                bed_mass(i,j,nbed,nnew,ised)=                           &
     &                             bed_mass(i,j,nbed-1,nnew,ised)+      &
     &                             bed_mass(i,j,nbed  ,nnew,ised)
              END DO
!
!  Push layers down.
!
              DO k=nbed-1,2,-1
                bed(i,j,k,iporo)=bed(i,j,k-1,iporo)
                bed(i,j,k,iaged)=bed(i,j,k-1,iaged)
                DO ised =1,nst
                  bed_mass(i,j,k,nnew,ised)=bed_mass(i,j,k-1,nnew,ised)
                END DO
              END DO
!
!  Set new top layer parameters.
!
              DO ised=1,nst
                bed_mass(i,j,2,nnew,ised)=max(bed_mass(i,j,2,nnew,ised)-&
     &                                    dep_mass(i,ised),0.0_r8)
                bed_mass(i,j,1,nnew,ised)=dep_mass(i,ised)
              END DO
            END IF
          END IF !NBED=1
!
! Recalculate thickness and fractions for all layers.
!
          DO k=1,nbed
            cff3=0.0_r8
            DO ised=1,nst
              cff3=cff3+bed_mass(i,j,k,nnew,ised)
            END DO
            IF (cff3.eq.0.0_r8) THEN
              cff3=eps
            END IF
            bed(i,j,k,ithck)=0.0_r8
            DO ised=1,nst
              bed_frac(i,j,k,ised)=bed_mass(i,j,k,nnew,ised)/cff3
              bed(i,j,k,ithck)=max(bed(i,j,k,ithck)+                    &
     &                         bed_mass(i,j,k,nnew,ised)/               &
     &                         (srho(ised,ng)*                          &
     &                          (1.0_r8-bed(i,j,k,iporo))),0.0_r8)
            END DO
          END DO
        END DO
      END DO j_loop
!
!  End of Suspended Sediment only section.
!
# endif
!
!  Ensure top bed layer thickness is greater or equal than active layer
!  thickness. If need to add sed to top layer, then entrain from lower
!  levels. Create new layers at bottom to maintain Nbed.
!
      j_loop2 : DO j=jstr,jend
        DO i=istr,iend
!
!  Calculate active layer thickness, bottom(i,j,iactv).
!
          bottom(i,j,iactv)=max(0.0_r8,                                 &
     &                          0.007_r8*                               &
     &                          (tau_w(i,j)-bottom(i,j,itauc))*rho0)+   &
     &                          6.0_r8*bottom(i,j,isd50)
!
# ifdef SED_MORPH
!
! Apply morphology factor.
!
          bottom(i,j,iactv)=max(bottom(i,j,iactv)*morph_fac(1,ng),      &
     &                          bottom(i,j,iactv))
# endif
!
          IF (bottom(i,j,iactv).gt.bed(i,j,1,ithck)) THEN
            IF (nbed.eq.1) THEN
              bottom(i,j,iactv)=bed(i,j,1,ithck)
            ELSE
              thck_to_add=bottom(i,j,iactv)-bed(i,j,1,ithck)
              thck_avail=0.0_r8
              ksed=1                                        ! initialize
              DO k=2,nbed
                IF (thck_avail.lt.thck_to_add) THEN
                  thck_avail=thck_avail+bed(i,j,k,ithck)
                  ksed=k
                END IF
              END DO
!
!  Catch here if there was not enough bed material.
!
              IF (thck_avail.lt.thck_to_add) THEN
                bottom(i,j,iactv)=bed(i,j,1,ithck)+thck_avail
                thck_to_add=thck_avail
              END IF
!
!  Update bed mass of top layer and fractional layer.
!
              cff2=max(thck_avail-thck_to_add,0.0_r8)/                  &
     &             max(bed(i,j,ksed,ithck),eps)
              DO ised=1,nst
                cff1=0.0_r8
                DO k=1,ksed
                  cff1=cff1+bed_mass(i,j,k,nnew,ised)
                END DO
                cff3=cff2*bed_mass(i,j,ksed,nnew,ised)
                bed_mass(i,j,1   ,nnew,ised)=cff1-cff3
                bed_mass(i,j,ksed,nnew,ised)=cff3
              END DO
!
!  Update thickness of fractional layer ksource_sed.
!
              bed(i,j,ksed,ithck)=max(thck_avail-thck_to_add,0.0_r8)
!
!  Upate bed fraction of top layer.
!
              cff3=0.0_r8
              DO ised=1,nst
                cff3=cff3+bed_mass(i,j,1,nnew,ised)
              END DO
              IF (cff3.eq.0.0_r8) THEN
                cff3=eps
              END IF
              DO ised=1,nst
                bed_frac(i,j,1,ised)=bed_mass(i,j,1,nnew,ised)/cff3
              END DO
!
!  Upate bed thickness of top layer.
!
              bed(i,j,1,ithck)=bottom(i,j,iactv)
!
!  Pull all layers closer to the surface.
!
              DO k=ksed,nbed
                ks=ksed-2
                bed(i,j,k-ks,ithck)=bed(i,j,k,ithck)
                bed(i,j,k-ks,iporo)=bed(i,j,k,iporo)
                bed(i,j,k-ks,iaged)=bed(i,j,k,iaged)
                DO ised=1,nst
                  bed_frac(i,j,k-ks,ised)=bed_frac(i,j,k,ised)
                  bed_mass(i,j,k-ks,nnew,ised)=bed_mass(i,j,k,nnew,ised)
                END DO
              END DO
!
!  Add new layers onto the bottom. Split what was in the bottom layer to
!  fill these new empty cells. ("ks" is the number of new layers).
!
              ks=ksed-2
              cff=1.0_r8/real(ks+1,r8)
              DO k=nbed,nbed-ks,-1
                bed(i,j,k,ithck)=bed(i,j,nbed-ks,ithck)*cff
                bed(i,j,k,iaged)=bed(i,j,nbed-ks,iaged)
                DO ised=1,nst
                  bed_frac(i,j,k,ised)=bed_frac(i,j,nbed-ks,ised)
                  bed_mass(i,j,k,nnew,ised)=                            &
     &                             bed_mass(i,j,nbed-ks,nnew,ised)*cff
                END DO
              END DO
            END IF  ! Nbed > 1
          END IF  ! increase top bed layer
        END DO
      END DO j_loop2
!
!-----------------------------------------------------------------------
! Store old bed thickness.
!-----------------------------------------------------------------------
!
# if defined SED_MORPH
      DO j=jstrr,jendr
        DO i=istrr,iendr
            bed_thick(i,j,nnew)=0.0_r8
            DO k=1,nbed
              bed_thick(i,j,nnew)=bed_thick(i,j,nnew)+                  &
     &                            bed(i,j,k,ithck)
            END DO
          END DO
        END DO
        IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
          CALL exchange_r2d_tile (ng, tile,                             &
     &                            lbi, ubi, lbj, ubj,                   &
     &                            bed_thick(:,:,nnew))
        END IF
# endif
!
!-----------------------------------------------------------------------
!  Apply periodic or gradient boundary conditions to property arrays.
!-----------------------------------------------------------------------
!
      DO ised=1,nst
        CALL bc_r3d_tile (ng, tile,                                     &
     &                    lbi, ubi, lbj, ubj, 1, nbed,                  &
     &                    bed_frac(:,:,:,ised))
        CALL bc_r3d_tile (ng, tile,                                     &
     &                    lbi, ubi, lbj, ubj, 1, nbed,                  &
     &                    bed_mass(:,:,:,nnew,ised))
      END DO
# ifdef DISTRIBUTE
      CALL mp_exchange4d (ng, tile, inlm, 2,                            &
     &                    lbi, ubi, lbj, ubj, 1, nbed, 1, nst,          &
     &                    nghostpoints,                                 &
     &                    ewperiodic(ng), nsperiodic(ng),               &
     &                    bed_frac,                                     &
     &                    bed_mass(:,:,:,nnew,:))
# endif
 
      DO i=1,mbedp
        CALL bc_r3d_tile (ng, tile,                                     &
     &                    lbi, ubi, lbj, ubj, 1, nbed,                  &
     &                    bed(:,:,:,i))
      END DO
# ifdef DISTRIBUTE
      CALL mp_exchange4d (ng, tile, inlm, 1,                            &
     &                    lbi, ubi, lbj, ubj, 1, nbed, 1, mbedp,        &
     &                    nghostpoints,                                 &
     &                    ewperiodic(ng), nsperiodic(ng),               &
     &                    bed)
# endif
!
      RETURN
      END SUBROUTINE sed_bed_tile
#endif
      END MODULE sed_bed_mod