131 & LBi, UBi, LBj, UBj, &
132 & IminS, ImaxS, JminS, JmaxS, &
136 & rmask, umask, vmask, &
145 & bustrcwmax, bvstrcwmax, &
146 & Dwave, Pwave_bot, &
150# if defined BEDLOAD_SOULSBY
154# if defined SED_MORPH
157# if defined BEDLOAD_MPM || defined BEDLOAD_SOULSBY
158 & h, om_r, om_u, on_r, on_v, &
161 & bed, bed_frac, bed_mass, &
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
186 real(r8),
intent(in) :: pm(LBi:,LBj:)
187 real(r8),
intent(in) :: pn(LBi:,LBj:)
189 real(r8),
intent(in) :: rmask(LBi:,LBj:)
190 real(r8),
intent(in) :: umask(LBi:,LBj:)
191 real(r8),
intent(in) :: vmask(LBi:,LBj:)
194 real(r8),
intent(in) :: rmask_wet(LBi:,LBj:)
196 real(r8),
intent(in) :: z_w(LBi:,LBj:,0:)
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:)
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:)
214# if defined SED_MORPH
215 real(r8),
intent(inout):: bed_thick(LBi:,LBj:,:)
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:,:)
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:,:)
232 real(r8),
intent(in) :: pm(LBi:UBi,LBj:UBj)
233 real(r8),
intent(in) :: pn(LBi:UBi,LBj:UBj)
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)
240 real(r8),
intent(in) :: rmask_wet(LBi:UBi,LBj:UBj)
242 real(r8),
intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
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)
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)
260# if defined SED_MORPH
261 real(r8),
intent(inout):: bed_thick(LBi:UBi,LBj:UBj,3)
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)
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)
281 integer :: i, ised, j, k
283 real(r8),
parameter :: eps = 1.0e-14_r8
285 real(r8) :: cff, cff1, cff2, cff3, cff4, cff5
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
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
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
303# if defined BEDLOAD_MPM
304 real(r8),
dimension(IminS:ImaxS,JminS:JmaxS) :: angleu
305 real(r8),
dimension(IminS:ImaxS,JminS:JmaxS) :: anglev
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
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
320# include "set_bounds.h"
326# if defined BEDLOAD_MPM || defined SUSPLOAD
330 tau_w(i,j)=sqrt(bustrcwmax(i,j)*bustrcwmax(i,j)+ &
331 & bvstrcwmax(i,j)*bvstrcwmax(i,j))
333 tau_w(i,j)=tau_w(i,j)*rmask_wet(i,j)
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)))
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)))
363 tau_w(i,j)=tau_w(i,j)*rmask_wet(i,j)
378 sed_angle=dtan(33.0_r8*
pi/180.0_r8)
384# if defined BEDLOAD_SOULSBY
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))
392 phic(i,j)=atan2(bvstrc(i,j),bustrc(i,j))
394 phicw(i,j)=1.5_r8*
pi-dwave(i,j)-phic(i,j)-angler(i,j)
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))
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
418 smgdr=sqrt(smgd)*
sd50(ised,ng)*
srho(ised,ng)
422# ifdef BEDLOAD_SOULSBY
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)
438 theta_wav=tau_wav(i,j)*osmgd+eps
439 theta_mean=tau_mean(i,j)*osmgd
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)
453 cff3=0.5_r8*(1.0_r8+sign(1.0_r8, &
454 & theta_max/
tau_ce(ised,ng)-1.0_r8))
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* &
466 IF (abs(phi_x2).gt.phi_x1)
THEN
471 bedld_x=phi_x*smgdr*cff3
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
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
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)
507 bedld=8.0_r8*(max((tau_w(i,j)*osmgd-0.047_r8), &
508 & 0.0_r8)**1.5_r8)*smgdr
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)
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
530 cff=abs(tau_mean(i,j))
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
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))
549 a_slopex=sed_angle/(cos(cff)*(sed_angle-dzdx))
551 a_slopey=sed_angle/(cos(cff)*(sed_angle-dzdy))
555 fx_r(i,j)=fx_r(i,j)*a_slopex
556 fe_r(i,j)=fe_r(i,j)*a_slopey
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))/ &
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))/ &
587 & sign(1.0_r8,fe_r(i,j))
594 IF (
domain(ng)%Western_Edge(tile))
THEN
596 fx_r(istr-1,j)=fx_r(istr,j)
597 fe_r(istr-1,j)=fe_r(istr,j)
602 IF (
domain(ng)%Eastern_Edge(tile))
THEN
604 fx_r(iend+1,j)=fx_r(iend,j)
605 fe_r(iend+1,j)=fe_r(iend,j)
611 IF (
domain(ng)%Southern_Edge(tile))
THEN
613 fx_r(i,jstr-1)=fx_r(i,jstr)
614 fe_r(i,jstr-1)=fe_r(i,jstr)
619 IF (
domain(ng)%Northern_Edge(tile))
THEN
621 fx_r(i,jend+1)=fx_r(i,jend)
622 fe_r(i,jend+1)=fe_r(i,jend)
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 ))
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 ))
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))
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))
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)))
680 fx(i,j)=fx(i,j)*umask(i,j)
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)))
695 fe(i,j)=fe(i,j)*vmask(i,j)
723 IF (
domain(ng)%Western_Edge(tile))
THEN
732 IF (
domain(ng)%Eastern_Edge(tile))
THEN
742 IF (
domain(ng)%Southern_Edge(tile))
THEN
751 IF (
domain(ng)%Northern_Edge(tile))
THEN
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)- &
768# if !defined SUSPLOAD
770 bed_mass(i,j,k,nnew,ised)=bed_mass(i,j,k,nstp,ised)
774 & cff/(
srho(ised,ng)* &
775 & (1.0_r8-bed(i,j,1,
iporo))), &
790 bedldu(i,j,ised)=fx(i,j)*(pn(i-1,j)+pn(i,j))*cff
795 bedldv(i,j,ised)=fe(i,j)*(pm(i,j-1)+pm(i,j))*cff
808 cff3=cff3+bed_mass(i,j,1,nnew,ised)
810 IF (cff3.eq.0.0_r8)
THEN
814 bed_frac(i,j,1,ised)=bed_mass(i,j,1,nnew,ised)/cff3
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)
827 bottom(i,j,
itauc)=cff1
829 bottom(i,j,
iwsed)=cff3
830 bottom(i,j,
idens)=max(cff4,1050.0_r8)
841 & lbi, ubi, lbj, ubj, 1, nbed, &
842 & bed_frac(:,:,:,ised))
844 & lbi, ubi, lbj, ubj, 1, nbed, &
845 & bed_mass(:,:,:,nnew,ised))
849 & lbi, ubi, lbj, ubj, &
852 & lbi, ubi, lbj, ubj, &
859 & lbi, ubi, lbj, ubj, 1, nbed, 1, nst, &
863 & bed_mass(:,:,:,nnew,:))
867 & lbi, ubi, lbj, ubj, 1, nst, &
877 & lbi, ubi, lbj, ubj, 1, nbed, &
882 & lbi, ubi, lbj, ubj, 1, nbed, 1, mbedp, &
889 & lbi, ubi, lbj, ubj, 1, mbotp, &
893 & lbi, ubi, lbj, ubj, 1, mbotp, &