sed_settling_mod Module Reference

Functions/Subroutines
subroutine, public	sed_settling (ng, tile)
subroutine	sed_settling_tile (ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, nstp, nnew, hz, z_w, settling_flux, t)

Function/Subroutine Documentation

sed_settling()

subroutine, public sed_settling_mod::sed_settling	(	integer, intent(in)	ng,
		integer, intent(in)	tile )

Definition at line 44 of file sed_settling.F.

!***********************************************************************
!
      USE mod_param
      USE mod_forces
      USE mod_grid
      USE mod_ocean
      USE mod_sedbed
      USE mod_stepping
!
!  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_settling_tile (ng, tile,                                 &
     &                        lbi, ubi, lbj, ubj,                       &
     &                        imins, imaxs, jmins, jmaxs,               &
     &                        nstp(ng), nnew(ng),                       &
     &                        grid(ng) % Hz,                            &
     &                        grid(ng) % z_w,                           &
     &                        sedbed(ng) % settling_flux,               &
     &                        ocean(ng) % t)
# ifdef PROFILE
      CALL wclock_off (ng, inlm, 16, __line__, myfile)
# endif
!
      RETURN

References mod_grid::grid, mod_param::inlm, mod_stepping::nnew, mod_stepping::nstp, mod_ocean::ocean, sed_settling_tile(), mod_sedbed::sedbed, wclock_off(), and wclock_on().

Referenced by sediment_mod::sediment().

Here is the call graph for this function:

Here is the caller graph for this function:

sed_settling_tile()

subroutine sed_settling_mod::sed_settling_tile ( integer, intent(in) ng, integer, intent(in) tile, integer, intent(in) lbi, integer, intent(in) ubi, integer, intent(in) lbj, integer, intent(in) ubj, integer, intent(in) imins, integer, intent(in) imaxs, integer, intent(in) jmins, integer, intent(in) jmaxs, integer, intent(in) nstp, integer, intent(in) nnew, real(r8), dimension(lbi:,lbj:,:), intent(in) hz, real(r8), dimension(lbi:,lbj:,0:), intent(in) z_w, real(r8), dimension(lbi:,lbj:,:), intent(inout) settling_flux, real(r8), dimension(lbi:,lbj:,:,:,:), intent(inout) t )

private

Definition at line 84 of file sed_settling.F.

!***********************************************************************
!
      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
!
# ifdef ASSUMED_SHAPE
      real(r8), intent(in) :: Hz(LBi:,LBj:,:)
      real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
      real(r8), intent(inout) :: settling_flux(LBi:,LBj:,:)
      real(r8), intent(inout) :: t(LBi:,LBj:,:,:,:)
# else
      real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
      real(r8), intent(inout) :: settling_flux(LBi:UBi,LBj:UBj,NST)
      real(r8), intent(inout) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
# endif
!
!  Local variable declarations.
!
      integer :: i, indx, ised, j, k, ks
 
      real(r8) :: cff, cu, cffL, cffR, dltL, dltR
 
      integer, dimension(IminS:ImaxS,N(ng)) :: ksource
 
      real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FC
 
      real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv
      real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv2
      real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv3
      real(r8), dimension(IminS:ImaxS,N(ng)) :: qc
      real(r8), dimension(IminS:ImaxS,N(ng)) :: qR
      real(r8), dimension(IminS:ImaxS,N(ng)) :: qL
      real(r8), dimension(IminS:ImaxS,N(ng)) :: WR
      real(r8), dimension(IminS:ImaxS,N(ng)) :: WL
 
# include "set_bounds.h"
!
!-----------------------------------------------------------------------
!  Add sediment vertical sinking (settling) term.
!-----------------------------------------------------------------------
!
!  Compute inverse thicknesses to avoid repeated divisions.
!
      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
        DO k=1,n(ng)-1
          DO i=istr,iend
            hz_inv2(i,k)=1.0_r8/(hz(i,j,k)+hz(i,j,k+1))
          END DO
        END DO
        DO k=2,n(ng)-1
          DO i=istr,iend
            hz_inv3(i,k)=1.0_r8/(hz(i,j,k-1)+hz(i,j,k)+hz(i,j,k+1))
          END DO
        END DO
!
!  Copy concentration of suspended sediment into scratch array "qc"
!  (q-central, restrict it to be positive) which is hereafter
!  interpreted as a set of grid-box averaged values for sediment
!  concentration.
!
        sed_loop: DO ised=1,nst
          indx=idsed(ised)
          DO k=1,n(ng)
            DO i=istr,iend
              qc(i,k)=t(i,j,k,nnew,indx)*hz_inv(i,k)
            END DO
          END DO
!
!-----------------------------------------------------------------------
!  Vertical sinking of suspended sediment.
!-----------------------------------------------------------------------
!
!  Reconstruct vertical profile of suspended sediment "qc" in terms
!  of a set of parabolic segments within each grid box. Then, compute
!  semi-Lagrangian flux due to sinking.
!
          DO k=n(ng)-1,1,-1
            DO i=istr,iend
              fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
            END DO
          END DO
          DO k=2,n(ng)-1
            DO i=istr,iend
              dltr=hz(i,j,k)*fc(i,k)
              dltl=hz(i,j,k)*fc(i,k-1)
              cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
              cffr=cff*fc(i,k)
              cffl=cff*fc(i,k-1)
!
!  Apply PPM monotonicity constraint to prevent oscillations within the
!  grid box.
!
              IF ((dltr*dltl).le.0.0_r8) THEN
                dltr=0.0_r8
                dltl=0.0_r8
              ELSE IF (abs(dltr).gt.abs(cffl)) THEN
                dltr=cffl
              ELSE IF (abs(dltl).gt.abs(cffr)) THEN
                dltl=cffr
              END IF
!
!  Compute right and left side values (qR,qL) of parabolic segments
!  within grid box Hz(k); (WR,WL) are measures of quadratic variations.
!
!  NOTE: Although each parabolic segment is monotonic within its grid
!        box, monotonicity of the whole profile is not guaranteed,
!        because qL(k+1)-qR(k) may still have different sign than
!        qc(k+1)-qc(k).  This possibility is excluded, after qL and qR
!        are reconciled using WENO procedure.
!
              cff=(dltr-dltl)*hz_inv3(i,k)
              dltr=dltr-cff*hz(i,j,k+1)
              dltl=dltl+cff*hz(i,j,k-1)
              qr(i,k)=qc(i,k)+dltr
              ql(i,k)=qc(i,k)-dltl
              wr(i,k)=(2.0_r8*dltr-dltl)**2
              wl(i,k)=(dltr-2.0_r8*dltl)**2
            END DO
          END DO
          cff=1.0e-14_r8
          DO k=2,n(ng)-2
            DO i=istr,iend
              dltl=max(cff,wl(i,k  ))
              dltr=max(cff,wr(i,k+1))
              qr(i,k)=(dltr*qr(i,k)+dltl*ql(i,k+1))/(dltr+dltl)
              ql(i,k+1)=qr(i,k)
            END DO
          END DO
          DO i=istr,iend
            fc(i,n(ng))=0.0_r8              ! no-flux boundary condition
# if defined LINEAR_CONTINUATION
            ql(i,n(ng))=qr(i,n(ng)-1)
            qr(i,n(ng))=2.0_r8*qc(i,n(ng))-ql(i,n(ng))
# elif defined NEUMANN
            ql(i,n(ng))=qr(i,n(ng)-1)
            qr(i,n(ng))=1.5_r8*qc(i,n(ng))-0.5_r8*ql(i,n(ng))
# else
            qr(i,n(ng))=qc(i,n(ng))         ! default strictly monotonic
            ql(i,n(ng))=qc(i,n(ng))         ! conditions
            qr(i,n(ng)-1)=qc(i,n(ng))
# endif
# if defined LINEAR_CONTINUATION
            qr(i,1)=ql(i,2)
            ql(i,1)=2.0_r8*qc(i,1)-qr(i,1)
# elif defined NEUMANN
            qr(i,1)=ql(i,2)
            ql(i,1)=1.5_r8*qc(i,1)-0.5_r8*qr(i,1)
# else
            ql(i,2)=qc(i,1)                 ! bottom grid boxes are
            qr(i,1)=qc(i,1)                 ! re-assumed to be
            ql(i,1)=qc(i,1)                 ! piecewise constant.
# endif
          END DO
!
!  Apply monotonicity constraint again, since the reconciled interfacial
!  values may cause a non-monotonic behavior of the parabolic segments
!  inside the grid box.
!
          DO k=1,n(ng)
            DO i=istr,iend
              dltr=qr(i,k)-qc(i,k)
              dltl=qc(i,k)-ql(i,k)
              cffr=2.0_r8*dltr
              cffl=2.0_r8*dltl
              IF ((dltr*dltl).lt.0.0_r8) THEN
                dltr=0.0_r8
                dltl=0.0_r8
              ELSE IF (abs(dltr).gt.abs(cffl)) THEN
                dltr=cffl
              ELSE IF (abs(dltl).gt.abs(cffr)) THEN
                dltl=cffr
              END IF
              qr(i,k)=qc(i,k)+dltr
              ql(i,k)=qc(i,k)-dltl
            END DO
          END DO
!
!  After this moment reconstruction is considered complete. The next
!  stage is to compute vertical advective fluxes, FC. It is expected
!  that sinking may occurs relatively fast, the algorithm is designed
!  to be free of CFL criterion, which is achieved by allowing
!  integration bounds for semi-Lagrangian advective flux to use as
!  many grid boxes in upstream direction as necessary.
!
!  In the two code segments below, WL is the z-coordinate of the
!  departure point for grid box interface z_w with the same indices;
!  FC is the finite volume flux; ksource(:,k) is index of vertical
!  grid box which contains the departure point (restricted by N(ng)).
!  During the search: also add in content of whole grid boxes
!  participating in FC.
!
          cff=dt(ng)*abs(wsed(ised,ng))
          DO k=1,n(ng)
            DO i=istr,iend
              fc(i,k-1)=0.0_r8
              wl(i,k)=z_w(i,j,k-1)+cff
              wr(i,k)=hz(i,j,k)*qc(i,k)
              ksource(i,k)=k
            END DO
          END DO
          DO k=1,n(ng)
            DO ks=k,n(ng)-1
              DO i=istr,iend
                IF (wl(i,k).gt.z_w(i,j,ks)) THEN
                  ksource(i,k)=ks+1
                  fc(i,k-1)=fc(i,k-1)+wr(i,ks)
                END IF
              END DO
            END DO
          END DO
!
!  Finalize computation of flux: add fractional part.
!
          DO k=1,n(ng)
            DO i=istr,iend
              ks=ksource(i,k)
              cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
              fc(i,k-1)=fc(i,k-1)+                                      &
     &                  hz(i,j,ks)*cu*                                  &
     &                  (ql(i,ks)+                                      &
     &                   cu*(0.5_r8*(qr(i,ks)-ql(i,ks))-                &
     &                       (1.5_r8-cu)*                               &
     &                       (qr(i,ks)+ql(i,ks)-2.0_r8*qc(i,ks))))
            END DO
          END DO
          DO i=istr,iend
            DO k=1,n(ng)
              t(i,j,k,nnew,indx)=qc(i,k)*hz(i,j,k)+(fc(i,k)-fc(i,k-1))
            END DO
            settling_flux(i,j,ised)=fc(i,0)
          END DO
        END DO sed_loop
      END DO j_loop
!
      RETURN

References mod_scalars::dt, mod_sediment::idsed, and mod_sediment::wsed.

Referenced by sed_settling().

Here is the caller graph for this function: