1 | #include "cppdefs.h"
|
---|
2 |
|
---|
3 | MODULE sed_bed_mod
|
---|
4 |
|
---|
5 | #if defined NONLINEAR && defined SEDIMENT && !defined COHESIVE_BED
|
---|
6 | !
|
---|
7 | !svn $Id: sed_bed.F 233 2008-09-18 02:05:28Z arango $
|
---|
8 | !==================================================== John C. Warner ===
|
---|
9 | ! Copyright (c) 2002-2008 The ROMS/TOMS Group Hernan G. Arango !
|
---|
10 | ! Licensed under a MIT/X style license !
|
---|
11 | ! See License_ROMS.txt !
|
---|
12 | !=======================================================================
|
---|
13 | ! !
|
---|
14 | ! This routine computes sediment bed layer stratigraphy. !
|
---|
15 | ! !
|
---|
16 | ! Warner, J.C., C.R. Sherwood, R.P. Signell, C.K. Harris, and H.G. !
|
---|
17 | ! Arango, 2008: Development of a three-dimensional, regional, !
|
---|
18 | ! coupled wave, current, and sediment-transport model, Computers !
|
---|
19 | ! & Geosciences, 34, 1284-1306. !
|
---|
20 | ! !
|
---|
21 | !=======================================================================
|
---|
22 | !
|
---|
23 | implicit none
|
---|
24 |
|
---|
25 | PRIVATE
|
---|
26 | PUBLIC :: sed_bed
|
---|
27 |
|
---|
28 | CONTAINS
|
---|
29 | !
|
---|
30 | !***********************************************************************
|
---|
31 | SUBROUTINE sed_bed (ng, tile)
|
---|
32 | !***********************************************************************
|
---|
33 | !
|
---|
34 | USE mod_param
|
---|
35 | USE mod_forces
|
---|
36 | USE mod_grid
|
---|
37 | USE mod_ocean
|
---|
38 | USE mod_stepping
|
---|
39 | # ifdef BBL_MODEL
|
---|
40 | USE mod_bbl
|
---|
41 | # endif
|
---|
42 | !
|
---|
43 | ! Imported variable declarations.
|
---|
44 | !
|
---|
45 | integer, intent(in) :: ng, tile
|
---|
46 | !
|
---|
47 | ! Local variable declarations.
|
---|
48 | !
|
---|
49 | # include "tile.h"
|
---|
50 | !
|
---|
51 | # ifdef PROFILE
|
---|
52 | CALL wclock_on (ng, iNLM, 16)
|
---|
53 | # endif
|
---|
54 | CALL sed_bed_tile (ng, tile, &
|
---|
55 | & LBi, UBi, LBj, UBj, &
|
---|
56 | & IminS, ImaxS, JminS, JmaxS, &
|
---|
57 | & nstp(ng), nnew(ng), &
|
---|
58 | # ifdef WET_DRY
|
---|
59 | & GRID(ng) % rmask_wet, &
|
---|
60 | # endif
|
---|
61 | # ifdef BBL_MODEL
|
---|
62 | & BBL(ng) % bustrc, &
|
---|
63 | & BBL(ng) % bvstrc, &
|
---|
64 | & BBL(ng) % bustrw, &
|
---|
65 | & BBL(ng) % bvstrw, &
|
---|
66 | & BBL(ng) % bustrcwmax, &
|
---|
67 | & BBL(ng) % bvstrcwmax, &
|
---|
68 | # else
|
---|
69 | & FORCES(ng) % bustr, &
|
---|
70 | & FORCES(ng) % bvstr, &
|
---|
71 | # endif
|
---|
72 | & OCEAN(ng) % t, &
|
---|
73 | # ifdef SUSPLOAD
|
---|
74 | & OCEAN(ng) % ero_flux, &
|
---|
75 | & OCEAN(ng) % settling_flux, &
|
---|
76 | # endif
|
---|
77 | # if defined SED_MORPH
|
---|
78 | & GRID(ng) % bed_thick, &
|
---|
79 | # endif
|
---|
80 | & OCEAN(ng) % bed, &
|
---|
81 | & OCEAN(ng) % bed_frac, &
|
---|
82 | & OCEAN(ng) % bed_mass, &
|
---|
83 | & OCEAN(ng) % bottom)
|
---|
84 | # ifdef PROFILE
|
---|
85 | CALL wclock_off (ng, iNLM, 16)
|
---|
86 | # endif
|
---|
87 | RETURN
|
---|
88 | END SUBROUTINE sed_bed
|
---|
89 | !
|
---|
90 | !***********************************************************************
|
---|
91 | SUBROUTINE sed_bed_tile (ng, tile, &
|
---|
92 | & LBi, UBi, LBj, UBj, &
|
---|
93 | & IminS, ImaxS, JminS, JmaxS, &
|
---|
94 | & nstp, nnew, &
|
---|
95 | # ifdef WET_DRY
|
---|
96 | & rmask_wet, &
|
---|
97 | # endif
|
---|
98 | # ifdef BBL_MODEL
|
---|
99 | & bustrc, bvstrc, &
|
---|
100 | & bustrw, bvstrw, &
|
---|
101 | & bustrcwmax, bvstrcwmax, &
|
---|
102 | # else
|
---|
103 | & bustr, bvstr, &
|
---|
104 | # endif
|
---|
105 | & t, &
|
---|
106 | # ifdef SUSPLOAD
|
---|
107 | & ero_flux, settling_flux, &
|
---|
108 | # endif
|
---|
109 | # if defined SED_MORPH
|
---|
110 | & bed_thick, &
|
---|
111 | # endif
|
---|
112 | & bed, bed_frac, bed_mass, &
|
---|
113 | & bottom)
|
---|
114 | !***********************************************************************
|
---|
115 | !
|
---|
116 | USE mod_param
|
---|
117 | USE mod_scalars
|
---|
118 | USE mod_sediment
|
---|
119 | !
|
---|
120 | USE bc_3d_mod, ONLY : bc_r3d_tile
|
---|
121 | # if defined EW_PERIODIC || defined NS_PERIODIC
|
---|
122 | USE exchange_2d_mod, ONLY : exchange_r2d_tile
|
---|
123 | # endif
|
---|
124 | # ifdef DISTRIBUTE
|
---|
125 | USE mp_exchange_mod, ONLY : mp_exchange3d, mp_exchange4d
|
---|
126 | # endif
|
---|
127 | !
|
---|
128 | ! Imported variable declarations.
|
---|
129 | !
|
---|
130 | integer, intent(in) :: ng, tile
|
---|
131 | integer, intent(in) :: LBi, UBi, LBj, UBj
|
---|
132 | integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
|
---|
133 | integer, intent(in) :: nstp, nnew
|
---|
134 | !
|
---|
135 | # ifdef ASSUMED_SHAPE
|
---|
136 | # ifdef WET_DRY
|
---|
137 | real(r8), intent(in) :: rmask_wet(LBi:,LBj:)
|
---|
138 | # endif
|
---|
139 | # ifdef BBL_MODEL
|
---|
140 | real(r8), intent(in) :: bustrc(LBi:,LBj:)
|
---|
141 | real(r8), intent(in) :: bvstrc(LBi:,LBj:)
|
---|
142 | real(r8), intent(in) :: bustrw(LBi:,LBj:)
|
---|
143 | real(r8), intent(in) :: bvstrw(LBi:,LBj:)
|
---|
144 | real(r8), intent(in) :: bustrcwmax(LBi:,LBj:)
|
---|
145 | real(r8), intent(in) :: bvstrcwmax(LBi:,LBj:)
|
---|
146 | # else
|
---|
147 | real(r8), intent(in) :: bustr(LBi:,LBj:)
|
---|
148 | real(r8), intent(in) :: bvstr(LBi:,LBj:)
|
---|
149 | # endif
|
---|
150 | # if defined SED_MORPH
|
---|
151 | real(r8), intent(inout):: bed_thick(LBi:,LBj:,:)
|
---|
152 | # endif
|
---|
153 | real(r8), intent(inout) :: t(LBi:,LBj:,:,:,:)
|
---|
154 | # ifdef SUSPLOAD
|
---|
155 | real(r8), intent(inout) :: ero_flux(LBi:,LBj:,:)
|
---|
156 | real(r8), intent(inout) :: settling_flux(LBi:,LBj:,:)
|
---|
157 | # endif
|
---|
158 | real(r8), intent(inout) :: bed(LBi:,LBj:,:,:)
|
---|
159 | real(r8), intent(inout) :: bed_frac(LBi:,LBj:,:,:)
|
---|
160 | real(r8), intent(inout) :: bed_mass(LBi:,LBj:,:,:,:)
|
---|
161 | real(r8), intent(inout) :: bottom(LBi:,LBj:,:)
|
---|
162 | # else
|
---|
163 | # ifdef WET_DRY
|
---|
164 | real(r8), intent(in) :: rmask_wet(LBi:UBi,LBj:UBj)
|
---|
165 | # endif
|
---|
166 | # ifdef BBL_MODEL
|
---|
167 | real(r8), intent(in) :: bustrc(LBi:UBi,LBj:UBj)
|
---|
168 | real(r8), intent(in) :: bvstrc(LBi:UBi,LBj:UBj)
|
---|
169 | real(r8), intent(in) :: bustrw(LBi:UBi,LBj:UBj)
|
---|
170 | real(r8), intent(in) :: bvstrw(LBi:UBi,LBj:UBj)
|
---|
171 | real(r8), intent(in) :: bustrcwmax(LBi:UBi,LBj:UBj)
|
---|
172 | real(r8), intent(in) :: bvstrcwmax(LBi:UBi,LBj:UBj)
|
---|
173 | # else
|
---|
174 | real(r8), intent(in) :: bustr(LBi:UBi,LBj:UBj)
|
---|
175 | real(r8), intent(in) :: bvstr(LBi:UBi,LBj:UBj)
|
---|
176 | # endif
|
---|
177 | # if defined SED_MORPH
|
---|
178 | real(r8), intent(inout):: bed_thick(LBi:UBi,LBj:UBj,2)
|
---|
179 | # endif
|
---|
180 | real(r8), intent(inout) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
|
---|
181 | # ifdef SUSPLOAD
|
---|
182 | real(r8), intent(inout) :: ero_flux(LBi:UBi,LBj:UBj,NST)
|
---|
183 | real(r8), intent(inout) :: settling_flux(LBi:UBi,LBj:UBj,NST)
|
---|
184 | # endif
|
---|
185 | real(r8), intent(inout) :: bed(LBi:UBi,LBj:UBj,Nbed,MBEDP)
|
---|
186 | real(r8), intent(inout) :: bed_frac(LBi:UBi,LBj:UBj,Nbed,NST)
|
---|
187 | real(r8), intent(inout) :: bed_mass(LBi:UBi,LBj:UBj,Nbed,1:2,NST)
|
---|
188 | real(r8), intent(inout) :: bottom(LBi:UBi,LBj:UBj,MBOTP)
|
---|
189 | # endif
|
---|
190 | !
|
---|
191 | ! Local variable declarations.
|
---|
192 | !
|
---|
193 | # ifdef DISTRIBUTE
|
---|
194 | # ifdef EW_PERIODIC
|
---|
195 | logical :: EWperiodic=.TRUE.
|
---|
196 | # else
|
---|
197 | logical :: EWperiodic=.FALSE.
|
---|
198 | # endif
|
---|
199 | # ifdef NS_PERIODIC
|
---|
200 | logical :: NSperiodic=.TRUE.
|
---|
201 | # else
|
---|
202 | logical :: NSperiodic=.FALSE.
|
---|
203 | # endif
|
---|
204 | # endif
|
---|
205 | integer :: Ksed, i, ised, j, k, ks
|
---|
206 | integer :: bnew
|
---|
207 |
|
---|
208 | real(r8), parameter :: eps = 1.0E-14_r8
|
---|
209 |
|
---|
210 | real(r8) :: cff, cff1, cff2, cff3
|
---|
211 | real(r8) :: thck_avail, thck_to_add
|
---|
212 |
|
---|
213 | real(r8), dimension(IminS:ImaxS,NST) :: dep_mass
|
---|
214 | real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tau_w
|
---|
215 |
|
---|
216 | # include "set_bounds.h"
|
---|
217 |
|
---|
218 | # ifdef BEDLOAD
|
---|
219 | bnew=nnew
|
---|
220 | # else
|
---|
221 | bnew=nstp
|
---|
222 | # endif
|
---|
223 | !
|
---|
224 | !-----------------------------------------------------------------------
|
---|
225 | ! Compute sediment bed layer stratigraphy.
|
---|
226 | !-----------------------------------------------------------------------
|
---|
227 | !
|
---|
228 | # if defined BEDLOAD_MPM || defined SUSPLOAD
|
---|
229 | # ifdef BBL_MODEL
|
---|
230 | DO j=Jstr-1,Jend+1
|
---|
231 | DO i=Istr-1,Iend+1
|
---|
232 | tau_w(i,j)=SQRT(bustrcwmax(i,j)*bustrcwmax(i,j)+ &
|
---|
233 | & bvstrcwmax(i,j)*bvstrcwmax(i,j))
|
---|
234 | # ifdef WET_DRY
|
---|
235 | tau_w(i,j)=tau_w(i,j)*rmask_wet(i,j)
|
---|
236 | # endif
|
---|
237 | END DO
|
---|
238 | END DO
|
---|
239 | # else
|
---|
240 | # ifdef EW_PERIODIC
|
---|
241 | # define I_RANGE Istr-1,Iend+1
|
---|
242 | # else
|
---|
243 | # define I_RANGE MAX(Istr-1,1),MIN(Iend+1,Lm(ng))
|
---|
244 | # endif
|
---|
245 | # ifdef NS_PERIODIC
|
---|
246 | # define J_RANGE Jstr-1,Jend+1
|
---|
247 | # else
|
---|
248 | # define J_RANGE MAX(Jstr-1,1),MIN(Jend+1,Mm(ng))
|
---|
249 | # endif
|
---|
250 | DO i=I_RANGE
|
---|
251 | DO j=J_RANGE
|
---|
252 | tau_w(i,j)=0.5_r8*SQRT((bustr(i,j)+bustr(i+1,j))* &
|
---|
253 | & (bustr(i,j)+bustr(i+1,j))+ &
|
---|
254 | & (bvstr(i,j)+bvstr(i,j+1))* &
|
---|
255 | & (bvstr(i,j)+bvstr(i,j+1)))
|
---|
256 | # ifdef WET_DRY
|
---|
257 | tau_w(i,j)=tau_w(i,j)*rmask_wet(i,j)
|
---|
258 | # endif
|
---|
259 | END DO
|
---|
260 | END DO
|
---|
261 | # undef I_RANGE
|
---|
262 | # undef J_RANGE
|
---|
263 | # endif
|
---|
264 | # endif
|
---|
265 | !
|
---|
266 | !-----------------------------------------------------------------------
|
---|
267 | ! Update bed properties according to ero_flux and dep_flux.
|
---|
268 | !-----------------------------------------------------------------------
|
---|
269 | !
|
---|
270 | # ifdef SUSPLOAD
|
---|
271 | J_LOOP : DO j=Jstr,Jend
|
---|
272 | SED_LOOP: DO ised=1,NST
|
---|
273 | !
|
---|
274 | ! The deposition and resuspension of sediment on the bottom "bed"
|
---|
275 | ! is due to precepitation flux FC(:,0), already computed, and the
|
---|
276 | ! resuspension (erosion, hence called ero_flux). The resuspension is
|
---|
277 | ! applied to the bottom-most grid box value qc(:,1) so the total mass
|
---|
278 | ! is conserved. Restrict "ero_flux" so that "bed" cannot go negative
|
---|
279 | ! after both fluxes are applied.
|
---|
280 | !
|
---|
281 | DO i=Istr,Iend
|
---|
282 | dep_mass(i,ised)=0.0_r8
|
---|
283 |
|
---|
284 | # ifdef SED_MORPH
|
---|
285 | !
|
---|
286 | ! Apply morphology factor.
|
---|
287 | !
|
---|
288 | ero_flux(i,j,ised)=ero_flux(i,j,ised)*morph_fac(ised,ng)
|
---|
289 | settling_flux(i,j,ised)=settling_flux(i,j,ised)* &
|
---|
290 | & morph_fac(ised,ng)
|
---|
291 | !
|
---|
292 | # endif
|
---|
293 | IF ((ero_flux(i,j,ised)-settling_flux(i,j,ised)).lt. &
|
---|
294 | & 0.0_r8) THEN
|
---|
295 | !
|
---|
296 | ! If first time step of deposit, then store deposit material in
|
---|
297 | ! temporary array, dep_mass.
|
---|
298 | !
|
---|
299 | IF ((time(ng).gt.(bed(i,j,1,iaged)+1.1_r8*dt(ng))).and. &
|
---|
300 | & (bed(i,j,1,ithck).gt.newlayer_thick(ng))) THEN
|
---|
301 | dep_mass(i,ised)=settling_flux(i,j,ised)- &
|
---|
302 | & ero_flux(i,j,ised)
|
---|
303 | END IF
|
---|
304 | bed(i,j,1,iaged)=time(ng)
|
---|
305 | END IF
|
---|
306 | !
|
---|
307 | ! Update bed mass arrays.
|
---|
308 | !
|
---|
309 | bed_mass(i,j,1,nnew,ised)=MAX(bed_mass(i,j,1,bnew,ised)- &
|
---|
310 | & (ero_flux(i,j,ised)- &
|
---|
311 | & settling_flux(i,j,ised)), &
|
---|
312 | & 0.0_r8)
|
---|
313 | DO k=2,Nbed
|
---|
314 | bed_mass(i,j,k,nnew,ised)=bed_mass(i,j,k,nstp,ised)
|
---|
315 | END DO
|
---|
316 | END DO
|
---|
317 | END DO SED_LOOP
|
---|
318 | !
|
---|
319 | ! If first time step of deposit, create new layer and combine bottom
|
---|
320 | ! two bed layers.
|
---|
321 | !
|
---|
322 | DO i=Istr,Iend
|
---|
323 | cff=0.0_r8
|
---|
324 | !
|
---|
325 | ! Determine if deposition ocurred here.
|
---|
326 | !
|
---|
327 | IF (Nbed.gt.1) THEN
|
---|
328 | DO ised=1,NST
|
---|
329 | cff=cff+dep_mass(i,ised)
|
---|
330 | END DO
|
---|
331 | IF (cff.gt.0.0_r8) THEN
|
---|
332 | !
|
---|
333 | ! Combine bottom layers.
|
---|
334 | !
|
---|
335 | bed(i,j,Nbed,iporo)=0.5_r8*(bed(i,j,Nbed-1,iporo)+ &
|
---|
336 | & bed(i,j,Nbed,iporo))
|
---|
337 | bed(i,j,Nbed,iaged)=0.5_r8*(bed(i,j,Nbed-1,iaged)+ &
|
---|
338 | & bed(i,j,Nbed,iaged))
|
---|
339 | DO ised=1,NST
|
---|
340 | bed_mass(i,j,Nbed,nnew,ised)= &
|
---|
341 | & bed_mass(i,j,Nbed-1,nnew,ised)+ &
|
---|
342 | & bed_mass(i,j,Nbed ,nnew,ised)
|
---|
343 | END DO
|
---|
344 | !
|
---|
345 | ! Push layers down.
|
---|
346 | !
|
---|
347 | DO k=Nbed-1,2,-1
|
---|
348 | bed(i,j,k,iporo)=bed(i,j,k-1,iporo)
|
---|
349 | bed(i,j,k,iaged)=bed(i,j,k-1,iaged)
|
---|
350 | DO ised =1,NST
|
---|
351 | bed_mass(i,j,k,nnew,ised)=bed_mass(i,j,k-1,nnew,ised)
|
---|
352 | END DO
|
---|
353 | END DO
|
---|
354 | !
|
---|
355 | ! Set new top layer parameters.
|
---|
356 | !
|
---|
357 | DO ised=1,NST
|
---|
358 | bed_mass(i,j,2,nnew,ised)=MAX(bed_mass(i,j,2,nnew,ised)- &
|
---|
359 | & dep_mass(i,ised),0.0_r8)
|
---|
360 | bed_mass(i,j,1,nnew,ised)=dep_mass(i,ised)
|
---|
361 | END DO
|
---|
362 | END IF
|
---|
363 | END IF !NBED=1
|
---|
364 | !
|
---|
365 | ! Recalculate thickness and fractions for all layers.
|
---|
366 | !
|
---|
367 | DO k=1,Nbed
|
---|
368 | cff3=0.0_r8
|
---|
369 | DO ised=1,NST
|
---|
370 | cff3=cff3+bed_mass(i,j,k,nnew,ised)
|
---|
371 | END DO
|
---|
372 | IF (cff3.eq.0.0_r8) THEN
|
---|
373 | cff3=eps
|
---|
374 | END IF
|
---|
375 | bed(i,j,k,ithck)=0.0_r8
|
---|
376 | DO ised=1,NST
|
---|
377 | bed_frac(i,j,k,ised)=bed_mass(i,j,k,nnew,ised)/cff3
|
---|
378 | bed(i,j,k,ithck)=MAX(bed(i,j,k,ithck)+ &
|
---|
379 | & bed_mass(i,j,k,nnew,ised)/ &
|
---|
380 | & (Srho(ised,ng)* &
|
---|
381 | & (1.0_r8-bed(i,j,k,iporo))),0.0_r8)
|
---|
382 | END DO
|
---|
383 | END DO
|
---|
384 | END DO
|
---|
385 | END DO J_LOOP
|
---|
386 | !
|
---|
387 | ! End of Suspended Sediment only section.
|
---|
388 | !
|
---|
389 | # endif
|
---|
390 | !
|
---|
391 | ! Ensure top bed layer thickness is greater or equal than active layer
|
---|
392 | ! thickness. If need to add sed to top layer, then entrain from lower
|
---|
393 | ! levels. Create new layers at bottom to maintain Nbed.
|
---|
394 | !
|
---|
395 | J_LOOP2 : DO j=Jstr,Jend
|
---|
396 | DO i=Istr,Iend
|
---|
397 | !
|
---|
398 | ! Calculate active layer thickness, bottom(i,j,iactv).
|
---|
399 | !
|
---|
400 | bottom(i,j,iactv)=MAX(0.0_r8, &
|
---|
401 | & 0.007_r8* &
|
---|
402 | & (tau_w(i,j)-bottom(i,j,itauc))*rho0)+ &
|
---|
403 | & 6.0_r8*bottom(i,j,isd50)
|
---|
404 | !
|
---|
405 | # ifdef SED_MORPH
|
---|
406 | !
|
---|
407 | ! Apply morphology factor.
|
---|
408 | !
|
---|
409 | bottom(i,j,iactv)=MAX(bottom(i,j,iactv)*morph_fac(1,ng), &
|
---|
410 | & bottom(i,j,iactv))
|
---|
411 | # endif
|
---|
412 | !
|
---|
413 | IF (bottom(i,j,iactv).gt.bed(i,j,1,ithck)) THEN
|
---|
414 | IF (Nbed.eq.1) THEN
|
---|
415 | bottom(i,j,iactv)=bed(i,j,1,ithck)
|
---|
416 | ELSE
|
---|
417 | thck_to_add=bottom(i,j,iactv)-bed(i,j,1,ithck)
|
---|
418 | thck_avail=0.0_r8
|
---|
419 | Ksed=1 ! initialize
|
---|
420 | DO k=2,Nbed
|
---|
421 | IF (thck_avail.lt.thck_to_add) THEN
|
---|
422 | thck_avail=thck_avail+bed(i,j,k,ithck)
|
---|
423 | Ksed=k
|
---|
424 | END IF
|
---|
425 | END DO
|
---|
426 | !
|
---|
427 | ! Catch here if there was not enough bed material.
|
---|
428 | !
|
---|
429 | IF (thck_avail.lt.thck_to_add) THEN
|
---|
430 | bottom(i,j,iactv)=bed(i,j,1,ithck)+thck_avail
|
---|
431 | thck_to_add=thck_avail
|
---|
432 | END IF
|
---|
433 | !
|
---|
434 | ! Update bed mass of top layer and fractional layer.
|
---|
435 | !
|
---|
436 | cff2=MAX(thck_avail-thck_to_add,0.0_r8)/ &
|
---|
437 | & MAX(bed(i,j,Ksed,ithck),eps)
|
---|
438 | DO ised=1,NST
|
---|
439 | cff1=0.0_r8
|
---|
440 | DO k=1,Ksed
|
---|
441 | cff1=cff1+bed_mass(i,j,k,nnew,ised)
|
---|
442 | END DO
|
---|
443 | cff3=cff2*bed_mass(i,j,Ksed,nnew,ised)
|
---|
444 | bed_mass(i,j,1 ,nnew,ised)=cff1-cff3
|
---|
445 | bed_mass(i,j,Ksed,nnew,ised)=cff3
|
---|
446 | END DO
|
---|
447 | !
|
---|
448 | ! Update thickness of fractional layer ksource_sed.
|
---|
449 | !
|
---|
450 | bed(i,j,Ksed,ithck)=MAX(thck_avail-thck_to_add,0.0_r8)
|
---|
451 | !
|
---|
452 | ! Upate bed fraction of top layer.
|
---|
453 | !
|
---|
454 | cff3=0.0_r8
|
---|
455 | DO ised=1,NST
|
---|
456 | cff3=cff3+bed_mass(i,j,1,nnew,ised)
|
---|
457 | END DO
|
---|
458 | IF (cff3.eq.0.0_r8) THEN
|
---|
459 | cff3=eps
|
---|
460 | END IF
|
---|
461 | DO ised=1,NST
|
---|
462 | bed_frac(i,j,1,ised)=bed_mass(i,j,1,nnew,ised)/cff3
|
---|
463 | END DO
|
---|
464 | !
|
---|
465 | ! Upate bed thickness of top layer.
|
---|
466 | !
|
---|
467 | bed(i,j,1,ithck)=bottom(i,j,iactv)
|
---|
468 | !
|
---|
469 | ! Pull all layers closer to the surface.
|
---|
470 | !
|
---|
471 | DO k=Ksed,Nbed
|
---|
472 | ks=Ksed-2
|
---|
473 | bed(i,j,k-ks,ithck)=bed(i,j,k,ithck)
|
---|
474 | bed(i,j,k-ks,iporo)=bed(i,j,k,iporo)
|
---|
475 | bed(i,j,k-ks,iaged)=bed(i,j,k,iaged)
|
---|
476 | DO ised=1,NST
|
---|
477 | bed_frac(i,j,k-ks,ised)=bed_frac(i,j,k,ised)
|
---|
478 | bed_mass(i,j,k-ks,nnew,ised)=bed_mass(i,j,k,nnew,ised)
|
---|
479 | END DO
|
---|
480 | END DO
|
---|
481 | !
|
---|
482 | ! Add new layers onto the bottom. Split what was in the bottom layer to
|
---|
483 | ! fill these new empty cells. ("ks" is the number of new layers).
|
---|
484 | !
|
---|
485 | ks=Ksed-2
|
---|
486 | cff=1.0_r8/REAL(ks+1,r8)
|
---|
487 | DO k=Nbed,Nbed-ks,-1
|
---|
488 | bed(i,j,k,ithck)=bed(i,j,Nbed-ks,ithck)*cff
|
---|
489 | bed(i,j,k,iaged)=bed(i,j,Nbed-ks,iaged)
|
---|
490 | DO ised=1,NST
|
---|
491 | bed_frac(i,j,k,ised)=bed_frac(i,j,Nbed-ks,ised)
|
---|
492 | bed_mass(i,j,k,nnew,ised)= &
|
---|
493 | & bed_mass(i,j,Nbed-ks,nnew,ised)*cff
|
---|
494 | END DO
|
---|
495 | END DO
|
---|
496 | END IF ! Nbed > 1
|
---|
497 | END IF ! increase top bed layer
|
---|
498 | END DO
|
---|
499 | END DO J_LOOP2
|
---|
500 | !
|
---|
501 | !-----------------------------------------------------------------------
|
---|
502 | ! Store old bed thickness.
|
---|
503 | !-----------------------------------------------------------------------
|
---|
504 | !
|
---|
505 | # if defined SED_MORPH
|
---|
506 | DO j=JstrR,JendR
|
---|
507 | DO i=IstrR,IendR
|
---|
508 | bed_thick(i,j,nnew)=0.0_r8
|
---|
509 | DO k=1,Nbed
|
---|
510 | bed_thick(i,j,nnew)=bed_thick(i,j,nnew)+ &
|
---|
511 | & bed(i,j,k,ithck)
|
---|
512 | END DO
|
---|
513 | END DO
|
---|
514 | END DO
|
---|
515 | # if defined EW_PERIODIC || defined NS_PERIODIC
|
---|
516 | CALL exchange_r2d_tile (ng, tile, &
|
---|
517 | & LBi, UBi, LBj, UBj, &
|
---|
518 | & bed_thick(:,:,nnew))
|
---|
519 | # endif
|
---|
520 | # endif
|
---|
521 | !
|
---|
522 | !-----------------------------------------------------------------------
|
---|
523 | ! Apply periodic or gradient boundary conditions to property arrays.
|
---|
524 | !-----------------------------------------------------------------------
|
---|
525 | !
|
---|
526 | DO ised=1,NST
|
---|
527 | CALL bc_r3d_tile (ng, tile, &
|
---|
528 | & LBi, UBi, LBj, UBj, 1, Nbed, &
|
---|
529 | & bed_frac(:,:,:,ised))
|
---|
530 | CALL bc_r3d_tile (ng, tile, &
|
---|
531 | & LBi, UBi, LBj, UBj, 1, Nbed, &
|
---|
532 | & bed_mass(:,:,:,nnew,ised))
|
---|
533 | END DO
|
---|
534 | # ifdef DISTRIBUTE
|
---|
535 | CALL mp_exchange4d (ng, tile, iNLM, 2, &
|
---|
536 | & LBi, UBi, LBj, UBj, 1, Nbed, 1, NST, &
|
---|
537 | & NghostPoints, EWperiodic, NSperiodic, &
|
---|
538 | & bed_frac, &
|
---|
539 | & bed_mass(:,:,:,nnew,:))
|
---|
540 | # endif
|
---|
541 |
|
---|
542 | DO i=1,MBEDP
|
---|
543 | CALL bc_r3d_tile (ng, tile, &
|
---|
544 | & LBi, UBi, LBj, UBj, 1, Nbed, &
|
---|
545 | & bed(:,:,:,i))
|
---|
546 | END DO
|
---|
547 | # ifdef DISTRIBUTE
|
---|
548 | CALL mp_exchange4d (ng, tile, iNLM, 1, &
|
---|
549 | & LBi, UBi, LBj, UBj, 1, Nbed, 1, MBEDP, &
|
---|
550 | & NghostPoints, EWperiodic, NSperiodic, &
|
---|
551 | & bed)
|
---|
552 | # endif
|
---|
553 |
|
---|
554 | RETURN
|
---|
555 | END SUBROUTINE sed_bed_tile
|
---|
556 | #endif
|
---|
557 | END MODULE sed_bed_mod
|
---|