163
164
169
171# ifdef BEDLOAD
173# endif
174# ifdef DISTRIBUTE
176# endif
177
178
179
180 integer, intent(in) :: ng, tile
181 integer, intent(in) :: LBi, UBi, LBj, UBj
182 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
183 integer, intent(in) :: nstp, nnew
184
185# ifdef ASSUMED_SHAPE
186 real(r8), intent(in) :: pm(LBi:,LBj:)
187 real(r8), intent(in) :: pn(LBi:,LBj:)
188# ifdef MASKING
189 real(r8), intent(in) :: rmask(LBi:,LBj:)
190 real(r8), intent(in) :: umask(LBi:,LBj:)
191 real(r8), intent(in) :: vmask(LBi:,LBj:)
192# endif
193# ifdef WET_DRY
194 real(r8), intent(in) :: rmask_wet(LBi:,LBj:)
195# endif
196 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
197# ifdef BBL_MODEL
198 real(r8), intent(in) :: bustrc(LBi:,LBj:)
199 real(r8), intent(in) :: bvstrc(LBi:,LBj:)
200 real(r8), intent(in) :: bustrw(LBi:,LBj:)
201 real(r8), intent(in) :: bvstrw(LBi:,LBj:)
202 real(r8), intent(in) :: bustrcwmax(LBi:,LBj:)
203 real(r8), intent(in) :: bvstrcwmax(LBi:,LBj:)
204 real(r8), intent(in) :: Dwave(LBi:,LBj:)
205 real(r8), intent(in) :: Pwave_bot(LBi:,LBj:)
206# endif
207 real(r8), intent(in) :: bustr(LBi:,LBj:)
208 real(r8), intent(in) :: bvstr(LBi:,LBj:)
209# if defined BEDLOAD_SOULSBY
210 real(r8), intent(in) :: Hwave(LBi:,LBj:)
211 real(r8), intent(in) :: Lwave(LBi:,LBj:)
212 real(r8), intent(in) :: angler(LBi:,LBj:)
213# endif
214# if defined SED_MORPH
215 real(r8), intent(inout):: bed_thick(LBi:,LBj:,:)
216# endif
217# if defined BEDLOAD_MPM || defined BEDLOAD_SOULSBY
218 real(r8), intent(in) :: h(LBi:,LBj:)
219 real(r8), intent(in) :: om_r(LBi:,LBj:)
220 real(r8), intent(in) :: om_u(LBi:,LBj:)
221 real(r8), intent(in) :: on_r(LBi:,LBj:)
222 real(r8), intent(in) :: on_v(LBi:,LBj:)
223 real(r8), intent(inout) :: bedldu(LBi:,LBj:,:)
224 real(r8), intent(inout) :: bedldv(LBi:,LBj:,:)
225# endif
226 real(r8), intent(inout) :: t(LBi:,LBj:,:,:,:)
227 real(r8), intent(inout) :: bed(LBi:,LBj:,:,:)
228 real(r8), intent(inout) :: bed_frac(LBi:,LBj:,:,:)
229 real(r8), intent(inout) :: bed_mass(LBi:,LBj:,:,:,:)
230 real(r8), intent(inout) :: bottom(LBi:,LBj:,:)
231# else
232 real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
233 real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
234# ifdef MASKING
235 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
236 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
237 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
238# endif
239# ifdef WET_DRY
240 real(r8), intent(in) :: rmask_wet(LBi:UBi,LBj:UBj)
241# endif
242 real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
243# ifdef BBL_MODEL
244 real(r8), intent(in) :: bustrc(LBi:UBi,LBj:UBj)
245 real(r8), intent(in) :: bvstrc(LBi:UBi,LBj:UBj)
246 real(r8), intent(in) :: bustrw(LBi:UBi,LBj:UBj)
247 real(r8), intent(in) :: bvstrw(LBi:UBi,LBj:UBj)
248 real(r8), intent(in) :: bustrcwmax(LBi:UBi,LBj:UBj)
249 real(r8), intent(in) :: bvstrcwmax(LBi:UBi,LBj:UBj)
250 real(r8), intent(in) :: Dwave(LBi:UBi,LBj:UBj)
251 real(r8), intent(in) :: Pwave_bot(LBi:UBi,LBj:UBj)
252# endif
253 real(r8), intent(in) :: bustr(LBi:UBi,LBj:UBj)
254 real(r8), intent(in) :: bvstr(LBi:UBi,LBj:UBj)
255# if defined BEDLOAD_SOULSBY
256 real(r8), intent(in) :: Hwave(LBi:UBi,LBj:UBj)
257 real(r8), intent(in) :: Lwave(LBi:UBi,LBj:UBj)
258 real(r8), intent(in) :: angler(LBi:UBi,LBj:UBj)
259# endif
260# if defined SED_MORPH
261 real(r8), intent(inout):: bed_thick(LBi:UBi,LBj:UBj,3)
262# endif
263# if defined BEDLOAD_MPM || defined BEDLOAD_SOULSBY
264 real(r8), intent(in) :: h(LBi:UBi,LBj:UBj)
265 real(r8), intent(in) :: om_r(LBi:UBi,LBj:UBj)
266 real(r8), intent(in) :: om_u(LBi:UBi,LBj:UBj)
267 real(r8), intent(in) :: on_r(LBi:UBi,LBj:UBj)
268 real(r8), intent(in) :: on_v(LBi:UBi,LBj:UBj)
269 real(r8), intent(inout) :: bedldu(LBi:UBi,LBj:UBj,NST)
270 real(r8), intent(inout) :: bedldv(LBi:UBi,LBj:UBj,NST)
271# endif
272 real(r8), intent(inout) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
273 real(r8), intent(inout) :: bed(LBi:UBi,LBj:UBj,Nbed,MBEDP)
274 real(r8), intent(inout) :: bed_frac(LBi:UBi,LBj:UBj,Nbed,NST)
275 real(r8), intent(inout) :: bed_mass(LBi:UBi,LBj:UBj,Nbed,1:2,NST)
276 real(r8), intent(inout) :: bottom(LBi:UBi,LBj:UBj,MBOTP)
277# endif
278
279
280
281 integer :: i, ised, j, k
282
283 real(r8), parameter :: eps = 1.0e-14_r8
284
285 real(r8) :: cff, cff1, cff2, cff3, cff4, cff5
286
287 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tau_w
288# ifdef BSTRESS_UPWIND
289 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tau_wX
290 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tau_wE
291# endif
292# ifdef BEDLOAD
293 real(r8) :: a_slopex, a_slopey, sed_angle
294 real(r8) :: bedld, bedld_mass, dzdx, dzdy
295 real(r8) :: smgd, smgdr, osmgd, Umag
296 real(r8) :: rhs_bed, Ua, Ra, phi, Clim
297
298 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FX
299 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FE
300 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FX_r
301 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FE_r
302# endif
303# if defined BEDLOAD_MPM
304 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: angleu
305 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: anglev
306# endif
307# if defined BEDLOAD_SOULSBY
308 real(r8) :: theta_mean, theta_wav, w_asym
309 real(r8) :: theta_max, theta_max1, theta_max2
310 real(r8) :: phi_x1, phi_x2, phi_x, phi_y, Dstp
311 real(r8) :: bedld_x, bedld_y, tau_cur, waven, wavec
312
313 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: phic
314 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: phicw
315 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tau_wav
316 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tau_mean
317 real(r8), parameter :: kdmax = 100.0_r8
318# endif
319
320# include "set_bounds.h"
321
322
323
324
325
326# if defined BEDLOAD_MPM || defined SUSPLOAD
327# ifdef BBL_MODEL
328 DO j=jstr-1,jend+1
329 DO i=istr-1,iend+1
330 tau_w(i,j)=sqrt(bustrcwmax(i,j)*bustrcwmax(i,j)+ &
331 & bvstrcwmax(i,j)*bvstrcwmax(i,j))
332# ifdef WET_DRY
333 tau_w(i,j)=tau_w(i,j)*rmask_wet(i,j)
334# endif
335 END DO
336 END DO
337# else
338 DO j=jstrm1,jendp1
339 DO i=istrm1,iendp1
340# ifdef BSTRESS_UPWIND
341 cff1=0.5_r8*(1.0_r8+sign(1.0_r8,bustr(i+1,j)))
342 cff2=0.5_r8*(1.0_r8-sign(1.0_r8,bustr(i+1,j)))
343 cff3=0.5_r8*(1.0_r8+sign(1.0_r8,bustr(i ,j)))
344 cff4=0.5_r8*(1.0_r8-sign(1.0_r8,bustr(i ,j)))
345 tau_wx(i,j)=cff3*(cff1*bustr(i,j)+ &
346 & cff2*0.5_r8*(bustr(i,j)+bustr(i+1,j)))+ &
347 & cff4*(cff2*bustr(i+1,j)+ &
348 & cff1*0.5_r8*(bustr(i,j)+bustr(i+1,j)))
349 cff1=0.5_r8*(1.0_r8+sign(1.0_r8,bvstr(i,j+1)))
350 cff2=0.5_r8*(1.0_r8-sign(1.0_r8,bvstr(i,j+1)))
351 cff3=0.5_r8*(1.0_r8+sign(1.0_r8,bvstr(i,j)))
352 cff4=0.5_r8*(1.0_r8-sign(1.0_r8,bvstr(i,j)))
353 tau_we(i,j)=cff3*(cff1*bvstr(i,j)+ &
354 & cff2*0.5_r8*(bvstr(i,j)+bvstr(i,j+1)))+ &
355 & cff4*(cff2*bvstr(i,j+1)+ &
356 & cff1*0.5_r8*(bvstr(i,j)+bvstr(i,j+1)))
357# endif
358 tau_w(i,j)=0.5_r8*sqrt((bustr(i,j)+bustr(i+1,j))* &
359 & (bustr(i,j)+bustr(i+1,j))+ &
360 & (bvstr(i,j)+bvstr(i,j+1))* &
361 & (bvstr(i,j)+bvstr(i,j+1)))
362# ifdef WET_DRY
363 tau_w(i,j)=tau_w(i,j)*rmask_wet(i,j)
364# endif
365 END DO
366 END DO
367# endif
368# endif
369
370# ifdef BEDLOAD
371
372
373
374
375
376
377
378 sed_angle=dtan(33.0_r8*
pi/180.0_r8)
379
380
381
382 DO j=jstrm1,jendp1
383 DO i=istrm1,iendp1
384# if defined BEDLOAD_SOULSBY
385
386
387
388
389 IF (bustrc(i,j).eq.0.0_r8) THEN
390 phic(i,j)=0.5_r8*
pi*sign(1.0_r8,bvstrc(i,j))
391 ELSE
392 phic(i,j)=atan2(bvstrc(i,j),bustrc(i,j))
393 ENDIF
394 phicw(i,j)=1.5_r8*
pi-dwave(i,j)-phic(i,j)-angler(i,j)
395
396
397
398 tau_cur=sqrt(bustrc(i,j)*bustrc(i,j)+ &
399 & bvstrc(i,j)*bvstrc(i,j))
400 tau_wav(i,j)=sqrt(bustrw(i,j)*bustrw(i,j)+ &
401 & bvstrw(i,j)*bvstrw(i,j))
402 tau_mean(i,j)=tau_cur*(1.0_r8+1.2_r8*((tau_wav(i,j)/ &
403 & (tau_cur+tau_wav(i,j)+eps))**3.2_r8))
404
405# elif defined BEDLOAD_MPM
406 cff1=0.5_r8*(bustr(i,j)+bustr(i+1,j))
407 cff2=0.5_r8*(bvstr(i,j)+bvstr(i,j+1))
408 umag=sqrt(cff1*cff1+cff2*cff2)+eps
409 angleu(i,j)=cff1/umag
410 anglev(i,j)=cff2/umag
411# endif
412 END DO
413 END DO
414
417 osmgd=1.0_r8/smgd
418 smgdr=sqrt(smgd)*
sd50(ised,ng)*
srho(ised,ng)
419
420 DO j=jstrm1,jendp1
421 DO i=istrm1,iendp1
422# ifdef BEDLOAD_SOULSBY
423
424
425
426 dstp=z_w(i,j,
n(ng))+h(i,j)
427 waven=2.0_r8*
pi/(lwave(i,j)+eps)
428 wavec=sqrt(
g/waven*tanh(waven*dstp))
429 cff4=min(waven*dstp,kdmax)
430 cff1=-0.1875_r8*wavec*(waven*dstp)**2/(sinh(cff4))**4
431 cff2=0.125_r8*
g*hwave(i,j)**2/(wavec*dstp+eps)
432 cff3=
pi*hwave(i,j)/(pwave_bot(i,j)*sinh(cff4)+eps)
433 w_asym=max(min((cff1-cff2)/cff3,0.2_r8),0.0_r8)
434 w_asym=0.0_r8
435
436
437
438 theta_wav=tau_wav(i,j)*osmgd+eps
439 theta_mean=tau_mean(i,j)*osmgd
440
441 cff1=theta_wav*(1.0_r8+w_asym)
442 cff2=theta_wav*(1.0_r8-w_asym)
443 theta_max1=sqrt((theta_mean+ &
444 & cff1*cos(phicw(i,j)))**2+ &
445 & (cff1*sin(phicw(i,j)))**2)
446 theta_max2=sqrt((theta_mean+ &
447 & cff2*cos(phicw(i,j)+
pi))**2+ &
448 & (cff2*sin(phicw(i,j)+
pi))**2)
449 theta_max=max(theta_max1,theta_max2)
450
451
452
453 cff3=0.5_r8*(1.0_r8+sign(1.0_r8, &
454 & theta_max/
tau_ce(ised,ng)-1.0_r8))
455
456
457
458
459 phi_x1=12.0_r8*sqrt(theta_mean)* &
460 & max((theta_mean-
tau_ce(ised,ng)),0.0_r8)
461 phi_x2=12.0_r8*(0.9534_r8+0.1907*cos(2.0_r8*phicw(i,j)))* &
462 & sqrt(theta_wav)*theta_mean+ &
463 & 12.0_r8*(0.229_r8*w_asym*theta_wav**1.5_r8* &
464 & cos(phicw(i,j)))
465
466 IF (abs(phi_x2).gt.phi_x1) THEN
467 phi_x=phi_x2
468 ELSE
469 phi_x=phi_x1
470 END IF
471 bedld_x=phi_x*smgdr*cff3
472
473 cff5=theta_wav**1.5_r8+1.5_r8*(theta_mean**1.5_r8)
474 phi_y=12.0_r8*0.1907_r8*theta_wav*theta_wav* &
475 & (theta_mean*sin(2.0_r8*phicw(i,j))+1.2_r8*w_asym* &
476 & theta_wav*sin(phicw(i,j)))/cff5*cff3
477 bedld_y=phi_y*smgdr
478
479
480
481
482 fx_r(i,j)=(bedld_x*cos(phic(i,j))-bedld_y*sin(phic(i,j)))* &
484 fe_r(i,j)=(bedld_x*sin(phic(i,j))+bedld_y*cos(phic(i,j)))* &
486# elif defined BEDLOAD_MPM
487# ifdef BSTRESS_UPWIND
488
489
490
491
492
493
494 bedld=8.0_r8*(max((abs(tau_wx(i,j))*osmgd-0.047_r8), &
495 & 0.0_r8)**1.5_r8)*smgdr* &
496 & sign(1.0_r8,tau_wx(i,j))
497 fx_r(i,j)=bedld*on_r(i,j)*
dt(ng)
498 bedld=8.0_r8*(max((abs(tau_we(i,j))*osmgd-0.047_r8), &
499 & 0.0_r8)**1.5_r8)*smgdr* &
500 & sign(1.0_r8,tau_we(i,j))
501 fe_r(i,j)=bedld*om_r(i,j)*
dt(ng)
502# else
503
504
505
506
507 bedld=8.0_r8*(max((tau_w(i,j)*osmgd-0.047_r8), &
508 & 0.0_r8)**1.5_r8)*smgdr
509
510
511
512
513 fx_r(i,j)=angleu(i,j)*bedld*on_r(i,j)*
dt(ng)
514 fe_r(i,j)=anglev(i,j)*bedld*om_r(i,j)*
dt(ng)
515# endif
516# endif
517
518
519
520 cff1=0.5_r8*(1.0_r8+sign(1.0_r8,fx_r(i,j)))
521 cff2=0.5_r8*(1.0_r8-sign(1.0_r8,fx_r(i,j)))
522# if defined SLOPE_NEMETH
523 dzdx=(h(i+1,j)-h(i ,j))/om_u(i+1,j)*cff1+ &
524 & (h(i-1,j)-h(i ,j))/om_u(i ,j)*cff2
525 dzdy=(h(i,j+1)-h(i,j ))/on_v(i,j+1)*cff1+ &
526 & (h(i,j-1)-h(i,j ))/on_v(i ,j)*cff2
527# ifdef BEDLOAD_MPM
528 cff=abs(tau_w(i,j))
529# else
530 cff=abs(tau_mean(i,j))
531# endif
532 a_slopex=0.3_r8*cff**0.5_r8*0.002_r8*dzdx+ &
533 & 0.3_r8*cff**1.5_r8*3.330_r8*dzdx
534 a_slopey=0.3_r8*cff**0.5_r8*0.002_r8*dzdy+ &
535 & 0.3_r8*cff**1.5_r8*3.330_r8*dzdy
536
537
538
539 fx_r(i,j)=fx_r(i,j)+a_slopex
540 fe_r(i,j)=fe_r(i,j)+a_slopey
541# elif defined SLOPE_LESSER
542 dzdx=min(((h(i+1,j)-h(i ,j))/om_u(i+1,j)*cff1+ &
543 & (h(i ,j)-h(i-1,j))/om_u(i ,j)*cff2),0.52_r8)* &
544 & sign(1.0_r8,fx_r(i,j))
545 dzdy=min(((h(i,j+1)-h(i,j ))/on_v(i,j+1)*cff1+ &
546 & (h(i,j )-h(i,j-1))/on_v(i ,j)*cff2),0.52_r8)* &
547 & sign(1.0_r8,fe_r(i,j))
548 cff=datan(dzdx)
549 a_slopex=sed_angle/(cos(cff)*(sed_angle-dzdx))
550 cff=datan(dzdy)
551 a_slopey=sed_angle/(cos(cff)*(sed_angle-dzdy))
552
553
554
555 fx_r(i,j)=fx_r(i,j)*a_slopex
556 fe_r(i,j)=fe_r(i,j)*a_slopey
557# endif
558
559
560# ifdef SED_MORPH
561
562
563
566
567# endif
568
569
570
571
574
575
576
577 bedld_mass=abs(fx_r(i,j))+abs(fe_r(i,j))+eps
578 fx_r(i,j)=min(abs(fx_r(i,j)), &
579 & bed_mass(i,j,1,nstp,ised)* &
580 & om_r(i,j)*on_r(i,j)*abs(fx_r(i,j))/ &
581 & bedld_mass)* &
582 & sign(1.0_r8,fx_r(i,j))
583 fe_r(i,j)=min(abs(fe_r(i,j)), &
584 & bed_mass(i,j,1,nstp,ised)* &
585 & om_r(i,j)*on_r(i,j)*abs(fe_r(i,j))/ &
586 & bedld_mass)* &
587 & sign(1.0_r8,fe_r(i,j))
588 END DO
589 END DO
590
591
592
594 IF (
domain(ng)%Western_Edge(tile))
THEN
595 DO j=jstrm1,jendp1
596 fx_r(istr-1,j)=fx_r(istr,j)
597 fe_r(istr-1,j)=fe_r(istr,j)
598 END DO
599 END IF
600 END IF
602 IF (
domain(ng)%Eastern_Edge(tile))
THEN
603 DO j=jstrm1,jendp1
604 fx_r(iend+1,j)=fx_r(iend,j)
605 fe_r(iend+1,j)=fe_r(iend,j)
606 END DO
607 END IF
608 END IF
609
611 IF (
domain(ng)%Southern_Edge(tile))
THEN
612 DO i=istrm1,iendp1
613 fx_r(i,jstr-1)=fx_r(i,jstr)
614 fe_r(i,jstr-1)=fe_r(i,jstr)
615 END DO
616 END IF
617 END IF
619 IF (
domain(ng)%Northern_Edge(tile))
THEN
620 DO i=istrm1,iendp1
621 fx_r(i,jend+1)=fx_r(i,jend)
622 fe_r(i,jend+1)=fe_r(i,jend)
623 END DO
624 END IF
625 END IF
626
629 IF (
domain(ng)%SouthWest_Corner(tile))
THEN
630 fx_r(istr-1,jstr-1)=0.5_r8*(fx_r(istr ,jstr-1)+ &
631 & fx_r(istr-1,jstr ))
632 fe_r(istr-1,jstr-1)=0.5_r8*(fe_r(istr ,jstr-1)+ &
633 & fe_r(istr-1,jstr ))
634 END IF
635 END IF
636
639 IF (
domain(ng)%SouthEast_Corner(tile))
THEN
640 fx_r(iend+1,jstr-1)=0.5_r8*(fx_r(iend ,jstr-1)+ &
641 & fx_r(iend+1,jstr ))
642 fe_r(iend+1,jstr-1)=0.5_r8*(fe_r(iend ,jstr-1)+ &
643 & fe_r(iend+1,jstr ))
644 END IF
645 END IF
646
649 IF (
domain(ng)%NorthWest_Corner(tile))
THEN
650 fx_r(istr-1,jend+1)=0.5_r8*(fx_r(istr-1,jend )+ &
651 & fx_r(istr ,jend+1))
652 fe_r(istr-1,jend+1)=0.5_r8*(fe_r(istr-1,jend )+ &
653 & fe_r(istr ,jend+1))
654 END IF
655 END IF
656
659 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
660 fx_r(iend+1,jend+1)=0.5_r8*(fx_r(iend+1,jend )+ &
661 & fx_r(iend ,jend+1))
662 fe_r(iend+1,jend+1)=0.5_r8*(fe_r(iend+1,jend )+ &
663 & fe_r(iend ,jend+1))
664 END IF
665 END IF
666
667
668
669 DO j=jstr-1,jend+1
670 DO i=istr,iend+1
671 cff1=0.5_r8*(1.0_r8+sign(1.0_r8,fx_r(i,j)))
672 cff2=0.5_r8*(1.0_r8-sign(1.0_r8,fx_r(i,j)))
673 fx(i,j)=0.5_r8*(1.0_r8+sign(1.0_r8,fx_r(i-1,j)))* &
674 & (cff1*fx_r(i-1,j)+ &
675 & cff2*0.5_r8*(fx_r(i-1,j)+fx_r(i,j)))+ &
676 & 0.5_r8*(1.0_r8-sign(1.0_r8,fx_r(i-1,j)))* &
677 & (cff2*fx_r(i ,j)+ &
678 & cff1*0.5_r8*(fx_r(i-1,j)+fx_r(i,j)))
679# ifdef MASKING
680 fx(i,j)=fx(i,j)*umask(i,j)
681# endif
682 END DO
683 END DO
684 DO j=jstr,jend+1
685 DO i=istr-1,iend+1
686 cff1=0.5_r8*(1.0_r8+sign(1.0_r8,fe_r(i,j)))
687 cff2=0.5_r8*(1.0_r8-sign(1.0_r8,fe_r(i,j)))
688 fe(i,j)=0.5_r8*(1.0_r8+sign(1.0_r8,fe_r(i,j-1)))* &
689 & (cff1*fe_r(i,j-1)+ &
690 & cff2*0.5_r8*(fe_r(i,j-1)+fe_r(i,j)))+ &
691 & 0.5_r8*(1.0_r8-sign(1.0_r8,fe_r(i,j-1)))* &
692 & (cff2*fe_r(i ,j)+ &
693 & cff1*0.5_r8*(fe_r(i,j-1)+fe_r(i,j)))
694# ifdef MASKING
695 fe(i,j)=fe(i,j)*vmask(i,j)
696# endif
697 END DO
698 END DO
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
723 IF (
domain(ng)%Western_Edge(tile))
THEN
725 DO j=jstr-1,jend+1
726 fx(istr,j)=0.0_r8
727 END DO
728 END IF
729 END IF
730 END IF
732 IF (
domain(ng)%Eastern_Edge(tile))
THEN
734 DO j=jstr-1,jend+1
735 fx(iend+1,j)=0.0_r8
736 END DO
737 END IF
738 END IF
739 END IF
740
742 IF (
domain(ng)%Southern_Edge(tile))
THEN
744 DO i=istr-1,iend+1
745 fe(i,jstr)=0.0_r8
746 END DO
747 END IF
748 END IF
749 END IF
751 IF (
domain(ng)%Northern_Edge(tile))
THEN
753 DO i=istr-1,iend+1
754 fe(i,jend+1)=0.0_r8
755 END DO
756 END IF
757 END IF
758 END IF
759
760
761
762 DO j=jstr,jend
763 DO i=istr,iend
764 cff=(fx(i+1,j)-fx(i,j)+ &
765 & fe(i,j+1)-fe(i,j))*pm(i,j)*pn(i,j)
766 bed_mass(i,j,1,nnew,ised)=max(bed_mass(i,j,1,nstp,ised)- &
767 & cff,0.0_r8)
768# if !defined SUSPLOAD
770 bed_mass(i,j,k,nnew,ised)=bed_mass(i,j,k,nstp,ised)
771 END DO
772# endif
774 & cff/(
srho(ised,ng)* &
775 & (1.0_r8-bed(i,j,1,
iporo))), &
776 & 0.0_r8)
777# ifdef MASKING
779# endif
780 END DO
781 END DO
782
783
784
785
786
788 DO j=jstrr,jendr
789 DO i=istr,iendr
790 bedldu(i,j,ised)=fx(i,j)*(pn(i-1,j)+pn(i,j))*cff
791 END DO
792 END DO
793 DO j=jstr,jendr
794 DO i=istrr,iendr
795 bedldv(i,j,ised)=fe(i,j)*(pm(i,j-1)+pm(i,j))*cff
796 END DO
797 END DO
798 END DO
799
800
801
802
803
804 DO j=jstr,jend
805 DO i=istr,iend
806 cff3=0.0_r8
808 cff3=cff3+bed_mass(i,j,1,nnew,ised)
809 END DO
810 IF (cff3.eq.0.0_r8) THEN
811 cff3=eps
812 END IF
814 bed_frac(i,j,1,ised)=bed_mass(i,j,1,nnew,ised)/cff3
815 END DO
816
817 cff1=1.0_r8
818 cff2=1.0_r8
819 cff3=1.0_r8
820 cff4=1.0_r8
822 cff1=cff1*
tau_ce(ised,ng)**bed_frac(i,j,1,ised)
823 cff2=cff2*
sd50(ised,ng)**bed_frac(i,j,1,ised)
824 cff3=cff3*(
wsed(ised,ng)+eps)**bed_frac(i,j,1,ised)
825 cff4=cff4*
srho(ised,ng)**bed_frac(i,j,1,ised)
826 END DO
827 bottom(i,j,
itauc)=cff1
829 bottom(i,j,
iwsed)=cff3
830 bottom(i,j,
idens)=max(cff4,1050.0_r8)
831 END DO
832 END DO
833# endif
834
835
836
837
838
841 & lbi, ubi, lbj, ubj, 1,
nbed, &
842 & bed_frac(:,:,:,ised))
844 & lbi, ubi, lbj, ubj, 1,
nbed, &
845 & bed_mass(:,:,:,nnew,ised))
846# ifdef BEDLOAD
849 & lbi, ubi, lbj, ubj, &
850 & bedldu(:,:,ised))
852 & lbi, ubi, lbj, ubj, &
853 & bedldv(:,:,ised))
854 END IF
855# endif
856 END DO
857# ifdef DISTRIBUTE
859 & lbi, ubi, lbj, ubj, 1,
nbed, 1,
nst, &
862 & bed_frac, &
863 & bed_mass(:,:,:,nnew,:))
864# ifdef BEDLOAD
867 & lbi, ubi, lbj, ubj, 1,
nst, &
870 & bedldu, bedldv)
871 END IF
872# endif
873# endif
874
877 & lbi, ubi, lbj, ubj, 1,
nbed, &
878 & bed(:,:,:,i))
879 END DO
880# ifdef DISTRIBUTE
882 & lbi, ubi, lbj, ubj, 1,
nbed, 1,
mbedp, &
885 & bed)
886# endif
887
889 & lbi, ubi, lbj, ubj, 1,
mbotp, &
890 & bottom)
891# ifdef DISTRIBUTE
893 & lbi, ubi, lbj, ubj, 1,
mbotp, &
896 & bottom)
897# endif
898
899 RETURN
subroutine bc_r3d_tile(ng, tile, lbi, ubi, lbj, ubj, lbk, ubk, a)
subroutine exchange_u2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
subroutine exchange_v2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
integer, dimension(:), allocatable istvar
integer, dimension(:), allocatable n
type(t_lbc), dimension(:,:,:), allocatable lbc
type(t_domain), dimension(:), allocatable domain
real(dp), dimension(:), allocatable dt
logical, dimension(:), allocatable ewperiodic
logical, dimension(:), allocatable nsperiodic
logical, dimension(:,:), allocatable compositegrid
integer, parameter isouth
real(r8), dimension(:), allocatable zob
integer, parameter inorth
real(r8), dimension(:), allocatable bedload_coeff
real(r8), dimension(:,:), allocatable srho
integer, dimension(:), allocatable idsed
real(r8), dimension(:,:), allocatable morph_fac
real(r8), dimension(:,:), allocatable sd50
real(r8), dimension(:,:), allocatable wsed
real(r8), dimension(:,:), allocatable tau_ce
subroutine mp_exchange4d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, lbk, ubk, lbt, ubt, nghost, ew_periodic, ns_periodic, a, b, c)
subroutine mp_exchange3d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, lbk, ubk, nghost, ew_periodic, ns_periodic, a, b, c, d)