ROMS
Loading...
Searching...
No Matches
sed_settling.F
Go to the documentation of this file.
1#include "cppdefs.h"
2#undef NEUMANN
3#undef LINEAR_CONTINUATION
4
5 MODULE sed_settling_mod
6
7#if defined NONLINEAR && defined SEDIMENT && defined SUSPLOAD
8!
9!git $Id$
10!=======================================================================
11! Copyright (c) 2002-2025 The ROMS Group !
12! Licensed under a MIT/X style license Hernan G. Arango !
13! See License_ROMS.md Alexander F. Shchepetkin !
14!==================================================== John C. Warner ===
15! !
16! This routine computes the vertical settling (sinking) of suspended !
17! sediment via a semi-Lagrangian advective flux algorithm. It uses a !
18! parabolic, vertical reconstructuion of the suspended sediment in !
19! the water column with PPT/WENO constraints to avoid oscillations. !
20! !
21! References: !
22! !
23! Colella, P. and P. Woodward, 1984: The piecewise parabolic method !
24! (PPM) for gas-dynamical simulations, J. Comp. Phys., 54, 174-201. !
25! !
26! Liu, X.D., S. Osher, and T. Chan, 1994: Weighted essentially !
27! nonoscillatory shemes, J. Comp. Phys., 115, 200-212. !
28! !
29! Warner, J.C., C.R. Sherwood, R.P. Signell, C.K. Harris, and H.G. !
30! Arango, 2008: Development of a three-dimensional, regional, !
31! coupled wave, current, and sediment-transport model, Computers !
32! & Geosciences, 34, 1284-1306. !
33! !
34!=======================================================================
35!
36 implicit none
37!
38 PRIVATE
39 PUBLIC :: sed_settling
40!
41 CONTAINS
42!
43!***********************************************************************
44 SUBROUTINE sed_settling (ng, tile)
45!***********************************************************************
46!
47 USE mod_param
48 USE mod_forces
49 USE mod_grid
50 USE mod_ocean
51 USE mod_sedbed
52 USE mod_stepping
53!
54! Imported variable declarations.
55!
56 integer, intent(in) :: ng, tile
57!
58! Local variable declarations.
59!
60 character (len=*), parameter :: myfile = &
61 & __FILE__
62!
63# include "tile.h"
64!
65# ifdef PROFILE
66 CALL wclock_on (ng, inlm, 16, __line__, myfile)
67# endif
68 CALL sed_settling_tile (ng, tile, &
69 & lbi, ubi, lbj, ubj, &
70 & imins, imaxs, jmins, jmaxs, &
71 & nstp(ng), nnew(ng), &
72 & grid(ng) % Hz, &
73 & grid(ng) % z_w, &
74 & sedbed(ng) % settling_flux, &
75 & ocean(ng) % t)
76# ifdef PROFILE
77 CALL wclock_off (ng, inlm, 16, __line__, myfile)
78# endif
79!
80 RETURN
81 END SUBROUTINE sed_settling
82!
83!***********************************************************************
84 SUBROUTINE sed_settling_tile (ng, tile, &
85 & LBi, UBi, LBj, UBj, &
86 & IminS, ImaxS, JminS, JmaxS, &
87 & nstp, nnew, &
88 & Hz, z_w, &
89 & settling_flux, &
90 & t)
91!***********************************************************************
92!
93 USE mod_param
94 USE mod_scalars
95 USE mod_sediment
96!
97!
98! Imported variable declarations.
99!
100 integer, intent(in) :: ng, tile
101 integer, intent(in) :: LBi, UBi, LBj, UBj
102 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
103 integer, intent(in) :: nstp, nnew
104!
105# ifdef ASSUMED_SHAPE
106 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
107 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
108 real(r8), intent(inout) :: settling_flux(LBi:,LBj:,:)
109 real(r8), intent(inout) :: t(LBi:,LBj:,:,:,:)
110# else
111 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
112 real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
113 real(r8), intent(inout) :: settling_flux(LBi:UBi,LBj:UBj,NST)
114 real(r8), intent(inout) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
115# endif
116!
117! Local variable declarations.
118!
119 integer :: i, indx, ised, j, k, ks
120
121 real(r8) :: cff, cu, cffL, cffR, dltL, dltR
122
123 integer, dimension(IminS:ImaxS,N(ng)) :: ksource
124
125 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FC
126
127 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv
128 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv2
129 real(r8), dimension(IminS:ImaxS,N(ng)) :: Hz_inv3
130 real(r8), dimension(IminS:ImaxS,N(ng)) :: qc
131 real(r8), dimension(IminS:ImaxS,N(ng)) :: qR
132 real(r8), dimension(IminS:ImaxS,N(ng)) :: qL
133 real(r8), dimension(IminS:ImaxS,N(ng)) :: WR
134 real(r8), dimension(IminS:ImaxS,N(ng)) :: WL
135
136# include "set_bounds.h"
137!
138!-----------------------------------------------------------------------
139! Add sediment vertical sinking (settling) term.
140!-----------------------------------------------------------------------
141!
142! Compute inverse thicknesses to avoid repeated divisions.
143!
144 j_loop : DO j=jstr,jend
145 DO k=1,n(ng)
146 DO i=istr,iend
147 hz_inv(i,k)=1.0_r8/hz(i,j,k)
148 END DO
149 END DO
150 DO k=1,n(ng)-1
151 DO i=istr,iend
152 hz_inv2(i,k)=1.0_r8/(hz(i,j,k)+hz(i,j,k+1))
153 END DO
154 END DO
155 DO k=2,n(ng)-1
156 DO i=istr,iend
157 hz_inv3(i,k)=1.0_r8/(hz(i,j,k-1)+hz(i,j,k)+hz(i,j,k+1))
158 END DO
159 END DO
160!
161! Copy concentration of suspended sediment into scratch array "qc"
162! (q-central, restrict it to be positive) which is hereafter
163! interpreted as a set of grid-box averaged values for sediment
164! concentration.
165!
166 sed_loop: DO ised=1,nst
167 indx=idsed(ised)
168 DO k=1,n(ng)
169 DO i=istr,iend
170 qc(i,k)=t(i,j,k,nnew,indx)*hz_inv(i,k)
171 END DO
172 END DO
173!
174!-----------------------------------------------------------------------
175! Vertical sinking of suspended sediment.
176!-----------------------------------------------------------------------
177!
178! Reconstruct vertical profile of suspended sediment "qc" in terms
179! of a set of parabolic segments within each grid box. Then, compute
180! semi-Lagrangian flux due to sinking.
181!
182 DO k=n(ng)-1,1,-1
183 DO i=istr,iend
184 fc(i,k)=(qc(i,k+1)-qc(i,k))*hz_inv2(i,k)
185 END DO
186 END DO
187 DO k=2,n(ng)-1
188 DO i=istr,iend
189 dltr=hz(i,j,k)*fc(i,k)
190 dltl=hz(i,j,k)*fc(i,k-1)
191 cff=hz(i,j,k-1)+2.0_r8*hz(i,j,k)+hz(i,j,k+1)
192 cffr=cff*fc(i,k)
193 cffl=cff*fc(i,k-1)
194!
195! Apply PPM monotonicity constraint to prevent oscillations within the
196! grid box.
197!
198 IF ((dltr*dltl).le.0.0_r8) THEN
199 dltr=0.0_r8
200 dltl=0.0_r8
201 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
202 dltr=cffl
203 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
204 dltl=cffr
205 END IF
206!
207! Compute right and left side values (qR,qL) of parabolic segments
208! within grid box Hz(k); (WR,WL) are measures of quadratic variations.
209!
210! NOTE: Although each parabolic segment is monotonic within its grid
211! box, monotonicity of the whole profile is not guaranteed,
212! because qL(k+1)-qR(k) may still have different sign than
213! qc(k+1)-qc(k). This possibility is excluded, after qL and qR
214! are reconciled using WENO procedure.
215!
216 cff=(dltr-dltl)*hz_inv3(i,k)
217 dltr=dltr-cff*hz(i,j,k+1)
218 dltl=dltl+cff*hz(i,j,k-1)
219 qr(i,k)=qc(i,k)+dltr
220 ql(i,k)=qc(i,k)-dltl
221 wr(i,k)=(2.0_r8*dltr-dltl)**2
222 wl(i,k)=(dltr-2.0_r8*dltl)**2
223 END DO
224 END DO
225 cff=1.0e-14_r8
226 DO k=2,n(ng)-2
227 DO i=istr,iend
228 dltl=max(cff,wl(i,k ))
229 dltr=max(cff,wr(i,k+1))
230 qr(i,k)=(dltr*qr(i,k)+dltl*ql(i,k+1))/(dltr+dltl)
231 ql(i,k+1)=qr(i,k)
232 END DO
233 END DO
234 DO i=istr,iend
235 fc(i,n(ng))=0.0_r8 ! no-flux boundary condition
236# if defined LINEAR_CONTINUATION
237 ql(i,n(ng))=qr(i,n(ng)-1)
238 qr(i,n(ng))=2.0_r8*qc(i,n(ng))-ql(i,n(ng))
239# elif defined NEUMANN
240 ql(i,n(ng))=qr(i,n(ng)-1)
241 qr(i,n(ng))=1.5_r8*qc(i,n(ng))-0.5_r8*ql(i,n(ng))
242# else
243 qr(i,n(ng))=qc(i,n(ng)) ! default strictly monotonic
244 ql(i,n(ng))=qc(i,n(ng)) ! conditions
245 qr(i,n(ng)-1)=qc(i,n(ng))
246# endif
247# if defined LINEAR_CONTINUATION
248 qr(i,1)=ql(i,2)
249 ql(i,1)=2.0_r8*qc(i,1)-qr(i,1)
250# elif defined NEUMANN
251 qr(i,1)=ql(i,2)
252 ql(i,1)=1.5_r8*qc(i,1)-0.5_r8*qr(i,1)
253# else
254 ql(i,2)=qc(i,1) ! bottom grid boxes are
255 qr(i,1)=qc(i,1) ! re-assumed to be
256 ql(i,1)=qc(i,1) ! piecewise constant.
257# endif
258 END DO
259!
260! Apply monotonicity constraint again, since the reconciled interfacial
261! values may cause a non-monotonic behavior of the parabolic segments
262! inside the grid box.
263!
264 DO k=1,n(ng)
265 DO i=istr,iend
266 dltr=qr(i,k)-qc(i,k)
267 dltl=qc(i,k)-ql(i,k)
268 cffr=2.0_r8*dltr
269 cffl=2.0_r8*dltl
270 IF ((dltr*dltl).lt.0.0_r8) THEN
271 dltr=0.0_r8
272 dltl=0.0_r8
273 ELSE IF (abs(dltr).gt.abs(cffl)) THEN
274 dltr=cffl
275 ELSE IF (abs(dltl).gt.abs(cffr)) THEN
276 dltl=cffr
277 END IF
278 qr(i,k)=qc(i,k)+dltr
279 ql(i,k)=qc(i,k)-dltl
280 END DO
281 END DO
282!
283! After this moment reconstruction is considered complete. The next
284! stage is to compute vertical advective fluxes, FC. It is expected
285! that sinking may occurs relatively fast, the algorithm is designed
286! to be free of CFL criterion, which is achieved by allowing
287! integration bounds for semi-Lagrangian advective flux to use as
288! many grid boxes in upstream direction as necessary.
289!
290! In the two code segments below, WL is the z-coordinate of the
291! departure point for grid box interface z_w with the same indices;
292! FC is the finite volume flux; ksource(:,k) is index of vertical
293! grid box which contains the departure point (restricted by N(ng)).
294! During the search: also add in content of whole grid boxes
295! participating in FC.
296!
297 cff=dt(ng)*abs(wsed(ised,ng))
298 DO k=1,n(ng)
299 DO i=istr,iend
300 fc(i,k-1)=0.0_r8
301 wl(i,k)=z_w(i,j,k-1)+cff
302 wr(i,k)=hz(i,j,k)*qc(i,k)
303 ksource(i,k)=k
304 END DO
305 END DO
306 DO k=1,n(ng)
307 DO ks=k,n(ng)-1
308 DO i=istr,iend
309 IF (wl(i,k).gt.z_w(i,j,ks)) THEN
310 ksource(i,k)=ks+1
311 fc(i,k-1)=fc(i,k-1)+wr(i,ks)
312 END IF
313 END DO
314 END DO
315 END DO
316!
317! Finalize computation of flux: add fractional part.
318!
319 DO k=1,n(ng)
320 DO i=istr,iend
321 ks=ksource(i,k)
322 cu=min(1.0_r8,(wl(i,k)-z_w(i,j,ks-1))*hz_inv(i,ks))
323 fc(i,k-1)=fc(i,k-1)+ &
324 & hz(i,j,ks)*cu* &
325 & (ql(i,ks)+ &
326 & cu*(0.5_r8*(qr(i,ks)-ql(i,ks))- &
327 & (1.5_r8-cu)* &
328 & (qr(i,ks)+ql(i,ks)-2.0_r8*qc(i,ks))))
329 END DO
330 END DO
331 DO i=istr,iend
332 DO k=1,n(ng)
333 t(i,j,k,nnew,indx)=qc(i,k)*hz(i,j,k)+(fc(i,k)-fc(i,k-1))
334 END DO
335 settling_flux(i,j,ised)=fc(i,0)
336 END DO
337 END DO sed_loop
338 END DO j_loop
339!
340 RETURN
341 END SUBROUTINE sed_settling_tile
342#endif
343 END MODULE sed_settling_mod
mod_forces
Definition mod_forces.F:2
mod_grid
Definition mod_grid.F:2
mod_grid::grid
type(t_grid), dimension(:), allocatable grid
Definition mod_grid.F:365
mod_ocean
Definition mod_ocean.F:2
mod_ocean::ocean
type(t_ocean), dimension(:), allocatable ocean
Definition mod_ocean.F:351
mod_param
Definition mod_param.F:2
mod_param::inlm
integer, parameter inlm
Definition mod_param.F:662
mod_scalars
Definition mod_scalars.F:2
mod_scalars::dt
real(dp), dimension(:), allocatable dt
Definition mod_scalars.F:269
mod_sedbed
Definition sedbed_mod.h:1
mod_sedbed::sedbed
type(t_sedbed), dimension(:), allocatable sedbed
Definition sedbed_mod.h:157
mod_sediment
Definition sediment_mod.h:1
mod_sediment::idsed
integer, dimension(:), allocatable idsed
Definition sediment_mod.h:88
mod_sediment::wsed
real(r8), dimension(:,:), allocatable wsed
Definition sediment_mod.h:175
mod_stepping
Definition mod_stepping.F:2
mod_stepping::nnew
integer, dimension(:), allocatable nnew
Definition mod_stepping.F:69
mod_stepping::nstp
integer, dimension(:), allocatable nstp
Definition mod_stepping.F:71
sed_settling_mod
Definition sed_settling.F:5
sed_settling_mod::sed_settling
subroutine, public sed_settling(ng, tile)
Definition sed_settling.F:45
sed_settling_mod::sed_settling_tile
subroutine sed_settling_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, nstp, nnew, hz, z_w, settling_flux, t)
Definition sed_settling.F:91
wclock_off
recursive subroutine wclock_off(ng, model, region, line, routine)
Definition timers.F:148
wclock_on
recursive subroutine wclock_on(ng, model, region, line, routine)
Definition timers.F:3