ROMS
Loading...
Searching...
No Matches
ssw_bbl.h
Go to the documentation of this file.
1#define M94WC
2#undef SGWC
3#define N92_RIPRUF
4#define CRS_FIX
5
6 MODULE bbl_mod
7!
8!git $Id$
9!================================================== Hernan G. Arango ===
10! Copyright (c) 2002-2025 The ROMS Group Chris Sherwood !
11! Licensed under a MIT/X style license Rich Signell !
12! See License_ROMS.md John C. Warner !
13!=======================================================================
14! !
15! This routine compute bottom stresses for the case when the wave !
16! solution in the wave boundary layer is based on a 2-layer eddy !
17! viscosity that is linear increasing above Zo and constant above !
18! Z1. !
19! !
20! Reference: !
21! !
22! Styles, R. and S.M. glenn, 2000: Modeling stratified wave and !
23! current bottom boundary layers in the continental shelf, JGR, !
24! 105, 24119-24139. !
25! !
26!=======================================================================
27!
28 implicit none
29!
30 PRIVATE
31 PUBLIC :: bblm
32!
33 CONTAINS
34!
35!***********************************************************************
36 SUBROUTINE bblm (ng, tile)
37!***********************************************************************
38!
39 USE mod_param
40 USE mod_bbl
41 USE mod_forces
42 USE mod_grid
43 USE mod_ocean
44 USE mod_sedbed
45 USE mod_stepping
46!
47! Imported variable declarations.
48!
49 integer, intent(in) :: ng, tile
50!
51! Local variable declarations.
52!
53 character (len=*), parameter :: MyFile = &
54 & __FILE__
55!
56#include "tile.h"
57!
58#ifdef PROFILE
59 CALL wclock_on (ng, inlm, 37, __line__, myfile)
60#endif
61 CALL ssw_bbl_tile (ng, tile, &
62 & lbi, ubi, lbj, ubj, &
63 & imins, imaxs, jmins, jmaxs, &
64 & nrhs(ng), &
65 & grid(ng) % h, &
66 & grid(ng) % z_r, &
67 & grid(ng) % z_w, &
68 & grid(ng) % angler, &
69 & grid(ng) % ZoBot, &
70#if defined SSW_CALC_UB
71 & forces(ng) % Hwave, &
72#else
73 & forces(ng) % Uwave_rms, &
74#endif
75 & forces(ng) % Dwave, &
76 & forces(ng) % Pwave_bot, &
77#ifdef BEDLOAD
78 & sedbed(ng) % bedldu, &
79 & sedbed(ng) % bedldv, &
80#endif
81 & sedbed(ng) % bottom, &
82 & ocean(ng) % rho, &
83 & ocean(ng) % u, &
84 & ocean(ng) % v, &
85#if defined SSW_LOGINT_STOKES
86 & ocean(ng) % u_stokes, &
87 & ocean(ng) % v_stokes, &
88#endif
89#if defined SSW_CALC_UB
90 & ocean(ng) % zeta, &
91#endif
92 & bbl(ng) % Iconv, &
93 & bbl(ng) % Ubot, &
94 & bbl(ng) % Vbot, &
95 & bbl(ng) % Ur, &
96 & bbl(ng) % Vr, &
97 & bbl(ng) % bustrc, &
98 & bbl(ng) % bvstrc, &
99 & bbl(ng) % bustrw, &
100 & bbl(ng) % bvstrw, &
101 & bbl(ng) % bustrcwmax, &
102 & bbl(ng) % bvstrcwmax, &
103 & forces(ng) % bustr, &
104 & forces(ng) % bvstr)
105#ifdef PROFILE
106 CALL wclock_off (ng, inlm, 37, __line__, myfile)
107#endif
108!
109 RETURN
110 END SUBROUTINE bblm
111!
112!***********************************************************************
113 SUBROUTINE ssw_bbl_tile (ng, tile, &
114 & LBi, UBi, LBj, UBj, &
115 & IminS, ImaxS, JminS, JmaxS, &
116 & nrhs, &
117 & h, z_r, z_w, angler, ZoBot, &
118#if defined SSW_CALC_UB
119 & Hwave, &
120#else
121 & Uwave_rms, &
122#endif
123 & Dwave, Pwave_bot, &
124#ifdef BEDLOAD
125 & bedldu, bedldv, &
126#endif
127 & bottom, rho, &
128 & u, v, &
129#if defined SSW_LOGINT_STOKES
130 & u_stokes, v_stokes, &
131#endif
132#if defined SSW_CALC_UB
133 & zeta, &
134#endif
135 & Iconv, &
136 & Ubot, Vbot, Ur, Vr, &
137 & bustrc, bvstrc, &
138 & bustrw, bvstrw, &
139 & bustrcwmax, bvstrcwmax, &
140 & bustr, bvstr)
141!***********************************************************************
142!
143 USE mod_param
144 USE mod_parallel
145 USE mod_iounits
146 USE mod_scalars
147 USE mod_sediment
148!
149 USE bc_2d_mod
150#ifdef DISTRIBUTE
151 USE mp_exchange_mod, ONLY : mp_exchange2d
152#endif
153
154!
155! Imported variable declarations.
156!
157 integer, intent(in) :: ng, tile
158 integer, intent(in) :: LBi, UBi, LBj, UBj
159 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
160 integer, intent(in) :: nrhs
161
162#ifdef ASSUMED_SHAPE
163 integer, intent(inout) :: Iconv(LBi:,LBj:)
164
165 real(r8), intent(in) :: h(LBi:,LBj:)
166 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
167 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
168 real(r8), intent(in) :: angler(LBi:,LBj:)
169 real(r8), intent(in) :: ZoBot(LBi:,LBj:)
170# if defined SSW_CALC_UB
171 real(r8), intent(in) :: Hwave(LBi:,LBj:)
172# else
173 real(r8), intent(in) :: Uwave_rms(LBi:,LBj:)
174# endif
175 real(r8), intent(in) :: Dwave(LBi:,LBj:)
176 real(r8), intent(in) :: Pwave_bot(LBi:,LBj:)
177# ifdef BEDLOAD
178 real(r8), intent(in) :: bedldu(LBi:,LBj:,:)
179 real(r8), intent(in) :: bedldv(LBi:,LBj:,:)
180# endif
181 real(r8), intent(inout) :: bottom(LBi:,LBj:,:)
182 real(r8), intent(in) :: rho(LBi:,LBj:,:)
183 real(r8), intent(in) :: u(LBi:,LBj:,:,:)
184 real(r8), intent(in) :: v(LBi:,LBj:,:,:)
185# if defined SSW_LOGINT_STOKES
186 real(r8), intent(in) :: u_stokes(LBi:,LBj:,:)
187 real(r8), intent(in) :: v_stokes(LBi:,LBj:,:)
188# endif
189# if defined SSW_CALC_UB
190 real(r8), intent(in) :: zeta(LBi:,LBj:,:)
191# endif
192 real(r8), intent(out) :: Ubot(LBi:,LBj:)
193 real(r8), intent(out) :: Vbot(LBi:,LBj:)
194 real(r8), intent(out) :: Ur(LBi:,LBj:)
195 real(r8), intent(out) :: Vr(LBi:,LBj:)
196 real(r8), intent(out) :: bustrc(LBi:,LBj:)
197 real(r8), intent(out) :: bvstrc(LBi:,LBj:)
198 real(r8), intent(out) :: bustrw(LBi:,LBj:)
199 real(r8), intent(out) :: bvstrw(LBi:,LBj:)
200 real(r8), intent(out) :: bustrcwmax(LBi:,LBj:)
201 real(r8), intent(out) :: bvstrcwmax(LBi:,LBj:)
202 real(r8), intent(out) :: bustr(LBi:,LBj:)
203 real(r8), intent(out) :: bvstr(LBi:,LBj:)
204#else
205 integer, intent(inout) :: Iconv(LBi:UBi,LBj:UBj)
206
207 real(r8), intent(in) :: h(LBi:UBi,LBj:UBj)
208 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
209 real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
210 real(r8), intent(in) :: angler(LBi:UBi,LBj:UBj)
211 real(r8), intent(in) :: ZoBot(LBi:UBi,LBj:UBj)
212# if defined SSW_CALC_UB
213 real(r8), intent(in) :: Hwave(LBi:UBi,LBj:UBj)
214# else
215 real(r8), intent(in) :: Uwave_rms(LBi:UBi,LBj:UBj)
216# endif
217 real(r8), intent(in) :: Dwave(LBi:UBi,LBj:UBj)
218 real(r8), intent(in) :: Pwave_bot(LBi:UBi,LBj:UBj)
219# ifdef BEDLOAD
220 real(r8), intent(in) :: bedldu(LBi:UBi,LBj:UBj,1:NST)
221 real(r8), intent(in) :: bedldv(LBi:UBi,LBj:UBj,1:NST)
222# endif
223 real(r8), intent(inout) :: bottom(LBi:UBi,LBj:UBj,MBOTP)
224 real(r8), intent(in) :: rho(LBi:UBi,LBj:UBj,N(ng))
225 real(r8), intent(in) :: u(LBi:UBi,LBj:UBj,N(ng),2)
226 real(r8), intent(in) :: v(LBi:UBi,LBj:UBj,N(ng),2)
227# if defined SSW_LOGINT_STOKES
228 real(r8), intent(in) :: u_stokes(LBi:UBi,LBj:UBj,N(ng))
229 real(r8), intent(in) :: v_stokes(LBi:UBi,LBj:UBj,N(ng))
230# endif
231# if defined SSW_CALC_UB
232 real(r8), intent(in) :: zeta(LBi:UBi,LBj:UBj,3)
233# endif
234 real(r8), intent(out) :: Ubot(LBi:UBi,LBj:UBj)
235 real(r8), intent(out) :: Vbot(LBi:UBi,LBj:UBj)
236 real(r8), intent(out) :: Ur(LBi:UBi,LBj:UBj)
237 real(r8), intent(out) :: Vr(LBi:UBi,LBj:UBj)
238 real(r8), intent(out) :: bustrc(LBi:UBi,LBj:UBj)
239 real(r8), intent(out) :: bvstrc(LBi:UBi,LBj:UBj)
240 real(r8), intent(out) :: bustrw(LBi:UBi,LBj:UBj)
241 real(r8), intent(out) :: bvstrw(LBi:UBi,LBj:UBj)
242 real(r8), intent(out) :: bustrcwmax(LBi:UBi,LBj:UBj)
243 real(r8), intent(out) :: bvstrcwmax(LBi:UBi,LBj:UBj)
244 real(r8), intent(out) :: bustr(LBi:UBi,LBj:UBj)
245 real(r8), intent(out) :: bvstr(LBi:UBi,LBj:UBj)
246#endif
247!
248! Local variable declarations.
249!
250 logical :: ITERATE
251
252 integer :: Iter, i, j, k
253
254 real(r8), parameter :: eps = 1.0e-10_r8
255
256 real(r8) :: Kbh, Kbh2, Kdh
257 real(r8) :: taucr, wsedr, tstar, coef_st
258 real(r8) :: coef_b1, coef_b2, coef_b3, d0
259 real(r8) :: dolam, dolam1, doeta1, doeta2, fdo_etaano
260 real(r8) :: lamorb, lamanorb
261 real(r8) :: m_ubr, m_wr, m_ucr, m_zr, m_phicw, m_kb
262 real(r8) :: m_ustrc, m_ustrwm, m_ustrr, m_fwc, m_zoa, m_dwc
263 real(r8) :: zo, Dstp
264 real(r8) :: Kb, Kdelta, Ustr
265 real(r8) :: anglec, anglew
266 real(r8) :: cff, cff1, cff2, cff3, og, fac, fac1, fac2
267 real(r8) :: sg_ab, sg_abokb, sg_a1, sg_b1, sg_chi, sg_c1, d50
268 real(r8) :: sg_epsilon, ssw_eta, sg_fofa, sg_fofb, sg_fofc, sg_fwm
269 real(r8) :: sg_kbs, ssw_lambda, sg_mu, sg_phicw, sg_ro, sg_row
270 real(r8) :: sg_shdnrm, sg_shld, sg_shldcr, sg_scf, rhos, sg_star
271 real(r8) :: sg_ub, sg_ubokur, sg_ubouc, sg_ubouwm, sg_ur
272 real(r8) :: sg_ustarc, sg_ustarcw, sg_ustarwm, sg_znot, sg_znotp
273 real(r8) :: sg_zr, sg_zrozn, sg_z1, sg_z1ozn, sg_z2, z1, z2
274 real(r8) :: zoMIN, zoMAX
275 real(r8) :: coef_fd
276
277 real(r8), parameter :: twopi=2.0_r8*pi
278!
279 real(r8), parameter :: absolute_zoMIN = 5.0d-5 ! in Harris-Wiberg
280!! real(r8), parameter :: absolute_zoMIN = 5.0d-8 ! in Harris-Wiberg
281 real(r8), parameter :: Cd_fd = 0.5_r8
282
283 real(r8), parameter :: K1 = 0.6666666666_r8 ! Coefficients for
284 real(r8), parameter :: K2 = 0.3555555555_r8 ! explicit
285 real(r8), parameter :: K3 = 0.1608465608_r8 ! wavenumber
286 real(r8), parameter :: K4 = 0.0632098765_r8 ! calculation
287 real(r8), parameter :: K5 = 0.0217540484_r8 ! (Dean and
288 real(r8), parameter :: K6 = 0.0065407983_r8 ! Dalrymple, 1991)
289
290 real(r8), parameter :: coef_a1=0.095_r8 ! Coefficients for
291 real(r8), parameter :: coef_a2=0.442_r8 ! ripple predictor
292 real(r8), parameter :: coef_a3=2.280_r8 ! (Wiberg-Harris)
293
294#if defined GM82_RIPRUF
295 real(r8), parameter :: ar = 27.7_r8/30.0_r8 ! Grant-Madsen (1982)
296#elif defined N92_RIPRUF
297 real(r8), parameter :: ar = 0.267_r8 ! Nielsen (1992)
298#elif defined R88_RIPRUF
299 real(r8), parameter :: ar = 0.533_r8 ! Raudkivi (1988)
300#else
301 no ripple roughness coeff. chosen
302#endif
303
304 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Ab
305 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Fwave_bot
306 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Tauc
307 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Tauw
308 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Taucwmax
309 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Ur_sg
310 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Vr_sg
311 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Ub
312 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Ucur
313 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Umag
314 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Vcur
315 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Zr
316 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: phic
317 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: phicw
318 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: rheight
319 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: rlength
320 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: u100
321 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: znot
322 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: znotc
323 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: zoN
324 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: zoST
325 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: zoBF
326 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: zoDEF
327 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: zoBIO
328 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: thck_wbl
329#if defined BEDLOAD_VANDERA_MADSEN_UDELTA
330 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: ksd_wbl
331 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: ustrc_wbl
332#endif
333#if defined BEDLOAD_VANDERA_MADSEN_UDELTA || \
334 defined bedload_vandera_direct_udelta
335 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: udelta_wbl
336 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Zr_wbl
337 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: phic_sgwbl
338#endif
339!
340 real(r8), dimension(1:N(ng)) :: Urz, Vrz
341#if defined BEDLOAD_VANDERA_MADSEN_UDELTA || \
342 defined bedload_vandera_direct_udelta
343 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Ur_sgwbl
344 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Vr_sgwbl
345 real(r8) :: Ucur_sgwbl, Vcur_sgwbl
346#endif
347!
348#include "set_bounds.h"
349!
350!-----------------------------------------------------------------------
351! Set currents above the bed.
352!-----------------------------------------------------------------------
353!
354 DO j=jstrv-1,jend
355 DO i=istru-1,iend
356!
357! Calculate bottom cell thickness
358!
359 zr(i,j)=z_r(i,j,1)-z_w(i,j,0)
360!
361#if defined SSW_LOGINT
362!
363 dstp=z_r(i,j,n(ng))-z_w(i,j,0)
364!
365# if defined SSW_LOGINT_DIRECT
366!
367! Cap minimum Zr 0.9*depth and user input.
368!
369 cff1=min(0.9_r8*dstp, max(zr(i,j), sg_zwbl(ng)))
370# elif defined SSW_LOGINT_WBL
371!
372! Cap minimum Zr 0.9*depth and computed wave bound layer
373!
374 cff1=min(0.98_r8*dstp,max(zr(i,j), bottom(i,j,idtbl)*1.1_r8))
375# else
376!
377! Use the original coded ssw_logint formulation
378!
379 cff1=sg_z1min
380# endif
381!
382! If using the logarithmic interpolation
383!
384 DO k=1,n(ng)
385 urz(k)=0.5_r8*(u(i,j,k,nrhs)+u(i+1,j,k,nrhs))
386 vrz(k)=0.5_r8*(v(i,j,k,nrhs)+v(i,j+1,k,nrhs))
387# ifdef SSW_LOGINT_STOKES
388 urz(k)=urz(k)+0.5_r8*(u_stokes(i,j,k)+u_stokes(i+1,j,k))
389 vrz(k)=vrz(k)+0.5_r8*(v_stokes(i,j,k)+v_stokes(i,j+1,k))
390# endif
391 END DO
392 CALL log_interp( n(ng), dstp, cff1, &
393 & urz, vrz, &
394 & z_r(i,j,:), z_w(i,j,:), &
395 & bottom(i,j,isd50), bottom(i,j,izapp), &
396 & zr(i,j), &
397 & ur_sg(i,j), vr_sg(i,j) )
398#else
399!
400! Regular method to get reference velocity for Madsen
401! without using log interp by the bottom cell.
402!
403 ur_sg(i,j)=0.5_r8*(u(i,j,1,nrhs)+u(i+1,j,1,nrhs))
404 vr_sg(i,j)=0.5_r8*(v(i,j,1,nrhs)+v(i,j+1,1,nrhs))
405# ifdef SSW_LOGINT_STOKES
406 ur_sg(i,j)=ur_sg(i,j)+0.5_r8*(u_stokes(i,j,1)+u_stokes(i+1,j,1))
407 vr_sg(i,j)=vr_sg(i,j)+0.5_r8*(v_stokes(i,j,1)+v_stokes(i,j+1,1))
408# endif
409#endif
410 END DO
411 END DO
412!
413!-----------------------------------------------------------------------
414! Compute bottom stresses.
415!-----------------------------------------------------------------------
416!
417 DO j=jstrv-1,jend
418 DO i=istru-1,iend
419!
420! Set bed wave orbital velocity and excursion amplitude. Use data
421! from wave models (SWAN) or use Dean and Dalrymple (1991) 6th-degree
422! polynomial to approximate wave number on shoaling water.
423
424 fwave_bot(i,j)=twopi/max(pwave_bot(i,j),1.0_r8)
425#ifdef SSW_CALC_UB
426 kdh=(h(i,j)+zeta(i,j,nrhs))*fwave_bot(i,j)**2/g
427 kbh2=kdh*kdh+ &
428 & kdh/(1.0_r8+kdh*(k1+kdh*(k2+kdh*(k3+kdh*(k4+ &
429 & kdh*(k5+k6*kdh))))))
430 kbh=sqrt(kbh2)
431 ab(i,j)=0.5_r8*hwave(i,j)/sinh(kbh)+eps
432 ub(i,j)=fwave_bot(i,j)*ab(i,j)+eps
433#else
434 ub(i,j)=max(uwave_rms(i,j),0.0_r8)+eps
435 ab(i,j)=ub(i,j)/fwave_bot(i,j)+eps
436#endif
437!
438! Compute bottom current magnitude at RHO-points.
439!
440 ucur(i,j)=ur_sg(i,j)
441 vcur(i,j)=vr_sg(i,j)
442!
443 umag(i,j)=sqrt(ucur(i,j)*ucur(i,j)+vcur(i,j)*vcur(i,j)+eps)
444!
445! Compute angle between currents and waves (radians)
446!
447 IF (ucur(i,j).eq.0.0_r8) THEN
448 phic(i,j)=0.5_r8*pi*sign(1.0_r8,vcur(i,j))
449 ELSE
450 phic(i,j)=atan2(vcur(i,j),ucur(i,j))
451 ENDIF
452 phicw(i,j)=1.5_r8*pi-dwave(i,j)-phic(i,j)-angler(i,j)
453 END DO
454 END DO
455!
456! Loop over RHO points.
457!
458 DO j=jstrv-1,jend
459 DO i=istru-1,iend
460!
461! Load sediment properties, stresses, and roughness from previous time
462! step (stresses are in m2 s-2).
463!
464 d50=bottom(i,j,isd50)
465 rhos=bottom(i,j,idens)/(rho(i,j,1)+1000.0_r8)
466 wsedr=bottom(i,j,iwsed)
467 taucr=bottom(i,j,itauc)
468 tauc(i,j)=sqrt(bustrc(i,j)**2+bvstrc(i,j)**2)
469 tauw(i,j)=sqrt(bustrw(i,j)**2+bvstrw(i,j)**2)
470 taucwmax(i,j)=sqrt( bustrcwmax(i,j)**2+bvstrcwmax(i,j)**2)
471!
472 rheight(i,j)=bottom(i,j,irhgt)
473 rlength(i,j)=bottom(i,j,irlen)
474 zomax=0.9_r8*zr(i,j)
475 zomin=max(absolute_zomin,2.5_r8*d50/30.0_r8)
476!
477! Initialize arrays.
478!
479 zon(i,j)=min(max(2.5_r8*d50/30.0_r8, zomin ),zomax)
480 zost(i,j)=0.0_r8
481 zobf(i,j)=0.0_r8
482 zobio(i,j)=0.0_r8
483!! zoDEF(i,j)=bottom(i,j,izdef)
484 zodef(i,j)=zobot(i,j)
485
486#ifdef SSW_CALC_ZNOT
487!
488! Calculate components of roughness and sum: zo = zoN + zoST + zoBF
489! Determine whether sediment is in motion. Use Shields criterion to
490! determine if sediment is mobile.
491!
492 tstar=taucwmax(i,j)/(taucr+eps)
493 IF (tstar.lt.1.0_r8) THEN ! no motion
494 zost(i,j)=0.0_r8
495 zobf(i,j)=ar*rheight(i,j)**2/(rlength(i,j)+eps)
496 ELSE
497!
498! Threshold of motion exceeded - calculate new zoST and zoBF
499! Calculate saltation roughness according to Wiberg & Rubin (1989)
500! (Eqn. 11 in Harris & Wiberg, 2001)
501! (d50 is in m, but this formula needs cm)
502!
503 coef_st=0.0204_r8*log(100.0_r8*d50+eps)**2+ &
504 & 0.0220_r8*log(100.0_r8*d50+eps)+0.0709_r8
505 zost(i,j)=0.056_r8*d50*0.68_r8*tstar/ &
506 & (1.0_r8+coef_st*tstar)
507 IF (zost(i,j).lt.0.0_r8) THEN
508 IF (master) THEN
509 WRITE (stdout,'(/,a)') &
510 & 'Warning zoST < 0: tstar, d50, coef_st:'
511 WRITE (stdout,*) tstar, d50, coef_st
512 END IF
513 END IF
514!
515! Calculate ripple height and wavelength.
516! Use Malarkey and Davies (2003) explict version of Wiberg and Harris.
517!
518 coef_b1=1.0_r8/coef_a1
519 coef_b2=0.5_r8*(1.0_r8 + coef_a2)*coef_b1
520 coef_b3=coef_b2**2-coef_a3*coef_b1
521 d0=2.0_r8*ab(i,j)
522 IF ((d0/d50).gt.13000.0_r8) THEN ! sheet flow
523 rheight(i,j)=0.0_r8
524 rlength(i,j)=535.0_r8*d50 ! does not matter since
525 ELSE ! rheight=0
526 dolam1=d0/(535.0_r8*d50)
527 doeta1=exp(coef_b2-sqrt(coef_b3-coef_b1*log(dolam1)))
528 lamorb=0.62_r8*d0
529 lamanorb=535.0_r8*d50
530 IF (doeta1.lt.20.0_r8) THEN
531 dolam=1.0_r8/0.62_r8
532 ELSE IF (doeta1.gt.100.0_r8) THEN
533 dolam=dolam1
534 ELSE
535 fdo_etaano=-log(lamorb/lamanorb)* &
536 & log(0.01_r8*doeta1)/log(5.0_r8)
537 dolam=dolam1*exp(-fdo_etaano)
538 END IF
539 doeta2=exp(coef_b2-sqrt(coef_b3-coef_b1*log(dolam)))
540 rheight(i,j)=d0/doeta2
541 rlength(i,j)=d0/dolam
542 END IF
543!
544! Value of ar can range from 0.3 to 3 (Soulsby, 1997, p. 124)
545!
546 zobf(i,j)=ar*rheight(i,j)**2/rlength(i,j)
547 END IF
548 zo=zon(i,j)
549# ifdef SSW_ZOBL
550 zo=zo+zost(i,j)
551# endif
552# ifdef SSW_ZORIP
553 zo=zo+zobf(i,j)
554# endif
555# ifdef SSW_ZOBIO
556 zo=zo+zobio(i,j)
557# endif
558#else
559 IF (zodef(i,j).lt.absolute_zomin) THEN
560 zodef(i,j)=absolute_zomin
561 IF (master) THEN
562 WRITE (stdout,*) ' Warning: default zo < 0.05 mm,', &
563 & ' replaced with: ', zodef
564 END IF
565 END IF
566 zo=zodef(i,j)
567#endif
568!
569! Compute stresses.
570!
571! Default stress calcs for pure currents
572!
573 zo=min(max(zo,zomin),zomax)
574!
575 cff1=vonkar/log(zr(i,j)/zo)
576 cff2=min(cdb_max,max(cdb_min,cff1*cff1))
577 tauc(i,j)=cff2*umag(i,j)*umag(i,j)
578 tauw(i,j)=0.0_r8
579 taucwmax(i,j)=tauc(i,j)
580 znot(i,j)=zo
581 znotc(i,j)=zo
582#if defined BEDLOAD_VANDERA_MADSEN_UDELTA
583 ksd_wbl(i,j)=zo
584 ustrc_wbl(i,j)=sqrt(tauc(i,j)+eps)
585#endif
586!
587 IF ((umag(i,j).le.eps).and.(ub(i,j).ge.eps)) THEN
588!
589! Pure waves - use wave friction factor approach from Madsen
590! (1994, eqns 32-33).
591!
592 sg_abokb=ab(i,j)/(30.0_r8*zo)
593 sg_fwm=0.3_r8
594 IF ((sg_abokb.gt.0.2_r8).and.(sg_abokb.le.100.0_r8)) THEN
595 sg_fwm=exp(-8.82_r8+7.02_r8*sg_abokb**(-0.078_r8))
596 ELSE IF (sg_abokb.gt.100.0_r8)THEN
597 sg_fwm=exp(-7.30_r8+5.61_r8*sg_abokb**(-0.109_r8))
598 END IF
599 tauc(i,j)= 0.0_r8
600 tauw(i,j)= 0.5_r8*sg_fwm*ub(i,j)*ub(i,j)
601 taucwmax(i,j)=tauw(i,j)
602 znot(i,j)=zo
603 znotc(i,j)=zo
604 ELSE IF ((umag(i,j).gt.eps).and.(ub(i,j).ge.eps).and. &
605 & ((zr(i,j)/zo).le.1.0_r8)) THEN
606!
607! Waves and currents, but zr <= zo.
608!
609 IF (master) THEN
610 WRITE (stdout,*) ' Warning: w-c calcs ignored because', &
611 & ' zr <= zo'
612 END IF
613 ELSE IF ((umag(i,j).gt.eps).and.(ub(i,j).ge.eps).and. &
614 & ((zr(i,j)/zo).gt.1.0_r8)) THEN
615!
616! Waves and currents, zr > zo.
617!
618#if defined SGWC
619 sg_zrozn=zr(i,j)/zo
620 sg_ubokur=ub(i,j)/(sg_kappa*umag(i,j))
621 sg_row=ab(i,j)/zo
622 sg_a1=1.0d-6
623 sg_phicw=phicw(i,j)
624 CALL sg_bstress (sg_row, sg_zrozn, sg_phicw, sg_ubokur, &
625 & sg_a1, sg_mu, sg_epsilon, sg_ro, sg_fofa)
626 sg_abokb=ab(i,j)/(30.0_r8*zo)
627 IF (sg_abokb.le.100.0_r8) THEN
628 sg_fwm=exp(-8.82_r8+7.02_r8*sg_abokb**(-0.078_r8))
629 ELSE
630 sg_fwm=exp(-7.30_r8+5.61_r8*sg_abokb**(-0.109_r8))
631 END IF
632 sg_ubouwm=sqrt(2.0_r8/sg_fwm)
633!
634! Determine the maximum ratio of wave over combined shear stresses,
635! sg_ubouwm (ub/ustarwm).
636!
637 CALL sg_purewave (sg_row, sg_ubouwm, sg_znotp, sg_ro)
638!
639! Set initial guess of the ratio of wave over shear stress, sg_c1
640! (ub/ustarc).
641!
642 sg_b1=sg_ubouwm
643 sg_fofb=-sg_fofa
644 sg_c1=0.5_r8*(sg_a1+sg_b1)
645 CALL sg_bstress (sg_row, sg_zrozn, sg_phicw, sg_ubokur, &
646 & sg_c1, sg_mu, sg_epsilon, sg_ro, sg_fofc)
647!
648! Solve PDE via bi-section method.
649!
650 iterate=.true.
651 DO iter=1,sg_n
652 IF (iterate) THEN
653 IF ((sg_fofb*sg_fofc).lt.0.0_r8) THEN
654 sg_a1=sg_c1
655 ELSE
656 sg_b1=sg_c1
657 END IF
658 sg_c1=0.5_r8*(sg_a1+sg_b1)
659 CALL sg_bstress (sg_row, sg_zrozn, sg_phicw, sg_ubokur, &
660 & sg_c1, sg_mu, sg_epsilon, sg_ro, &
661 & sg_fofc)
662 iterate=(sg_b1-sg_c1) .ge. sg_tol
663 IF (iterate) iconv(i,j)=iter
664 END IF
665 END DO
666 sg_ubouc=sg_c1
667!
668! Compute bottom shear stress magnitude (m/s).
669!
670 sg_ustarcw=ub(i,j)/sg_ubouc
671 sg_ustarwm=sg_mu*sg_ustarcw
672!! sg_ustarc=MIN(sg_ustarcdef,sg_epsilon*sg_ustarcw) !original
673 sg_ustarc=max(sqrt(tauc(i,j)),sg_epsilon*sg_ustarcw)
674 tauc(i,j)=sg_ustarc*sg_ustarc
675 tauw(i,j)=sg_ustarwm*sg_ustarwm
676 taucwmax(i,j)=sqrt((tauc(i,j)+ &
677 & tauw(i,j)*cos(phicw(i,j)))**2+ &
678 & (tauw(i,j)*sin(phicw(i,j)))**2)
679!
680! Compute apparent hydraulic roughness (m).
681!
682 IF (sg_epsilon.gt.0.0_r8) THEN
683 sg_z1=sg_alpha*sg_kappa*ab(i,j)/sg_ubouc
684 sg_z2=sg_z1/sg_epsilon
685 sg_z1ozn=sg_z1/zo
686 znotc(i,j)=sg_z2* &
687 & exp(-(1.0_r8-sg_epsilon+ &
688 & sg_epsilon*log(sg_z1ozn)))
689!
690! Compute mean (m/s) current at 100 cm above the bottom.
691!
692 IF (sg_z100.gt.sg_z2) THEN
693 u100(i,j)=sg_ustarc* &
694 & (log(sg_z100/sg_z2)+1.0_r8-sg_epsilon+ &
695 & sg_epsilon*log(sg_z1ozn))/sg_kappa
696 ELSE IF ((sg_z100.le.sg_z2).and.(zr(i,j).gt.sg_z1)) THEN
697 u100(i,j)=sg_ustarc*sg_epsilon* &
698 & (sg_z100/sg_z1-1.0_r8+log(sg_z1ozn))/sg_kappa
699 ELSE
700 u100(i,j)=sg_ustarc*sg_epsilon* &
701 & log(sg_z100/zo)/sg_kappa
702 END IF
703 END IF
704#elif defined M94WC
705 m_ubr=ub(i,j)
706 m_wr=fwave_bot(i,j)
707 m_ucr=umag(i,j)
708 m_zr=zr(i,j)
709 m_phicw=phicw(i,j)
710 m_kb=30.0_r8*zo
711 dstp=z_r(i,j,n(ng))-z_w(i,j,0)
712 CALL madsen94 (m_ubr, m_wr, m_ucr, &
713 & m_zr, m_phicw, m_kb, dstp, &
714 & m_ustrc, m_ustrwm, m_ustrr, m_fwc, m_zoa, &
715 & m_dwc)
716 tauc(i,j)=m_ustrc*m_ustrc
717 tauw(i,j)=m_ustrwm*m_ustrwm
718 taucwmax(i,j)=m_ustrr*m_ustrr
719 znotc(i,j)=min( m_zoa, zomax )
720 u100(i,j)=(m_ustrc/vonkar)*log(1.0_r8/m_zoa)
721 thck_wbl(i,j)=m_dwc
722#endif
723#if defined SSW_FORM_DRAG_COR
724 IF (rheight(i,j).gt.(zon(i,j)+zost(i,j))) THEN
725 coef_fd=0.5_r8*cd_fd*(rheight(i,j)/rlength(i,j))* &
726 & (1.0_r8/(vonkar*vonkar))* &
727 & (log(rheight(i,j)/ &
728 & (zon(i,j)+zost(i,j)))-1.0_r8)**2
729 taucwmax(i,j)=taucwmax(i,j)/(1.0_r8+coef_fd)
730 taucwmax(i,j)=taucwmax(i,j)*(1.0_r8+8.0_r8* &
731 & rheight(i,j)/rlength(i,j))
732 END IF
733#endif
734#if defined BEDLOAD_VANDERA_MADSEN_UDELTA
735 ksd_wbl(i,j)=m_zoa
736 ustrc_wbl(i,j)=m_ustrc
737#endif
738 END IF
739 END DO
740 END DO
741!
742#if defined BEDLOAD_VANDERA_MADSEN_UDELTA
743!
744! Find the near-bottom current velocity(udelta) at a given elevation.
745! Use the Madsen output of current shear stress, apparent roughness
746! to get udelta.
747! Find the angle at that near-bottom current velocity.
748!
749 DO j=jstr,jend
750 DO i=istr,iend
751!
752 dstp=z_r(i,j,n(ng))-z_w(i,j,0)
753!
754! Use wave boundary layer (wbl) thickness based on Madsen to get
755! near bottom current velocity.
756!
757 cff=min( 0.98_r8*dstp, thck_wbl(i,j) )
758!
759! Make sure that wbl is under total depth and greater than
760! apparent roughness.
761!
762 cff1=max(cff, 1.1_r8*ksd_wbl(i,j))
763 cff2=log(cff1/ksd_wbl(i,j))
764 udelta_wbl(i,j)=(ustrc_wbl(i,j)/vonkar)*cff2
765!
766 DO k=1,n(ng)
767 urz(k)=0.5_r8*(u(i,j,k,nrhs)+u(i+1,j,k,nrhs))
768 vrz(k)=0.5_r8*(v(i,j,k,nrhs)+v(i,j+1,k,nrhs))
769# ifdef SSW_LOGINT_STOKES
770 urz(k)=urz(k)+0.5_r8*(u_stokes(i,j,k)+u_stokes(i+1,j,k))
771 vrz(k)=vrz(k)+0.5_r8*(v_stokes(i,j,k)+v_stokes(i,j+1,k))
772# endif
773 END DO
774 CALL log_interp( n(ng), dstp, cff1, &
775 & urz, vrz, &
776 & z_r(i,j,:), z_w(i,j,:), &
777 & bottom(i,j,isd50), bottom(i,j,izapp), &
778 & zr_wbl(i,j), &
779 & ur_sgwbl(i,j), vr_sgwbl(i,j) )
780!
781! Compute bottom current magnitude at RHO-points.
782!
783 ucur_sgwbl=ur_sgwbl(i,j)
784 vcur_sgwbl=vr_sgwbl(i,j)
785!
786 IF (ucur_sgwbl.eq.0.0_r8) THEN
787 phic_sgwbl(i,j)=0.5_r8*pi*sign(1.0_r8,vcur_sgwbl)
788 ELSE
789 phic_sgwbl(i,j)=atan2(vcur_sgwbl,ucur_sgwbl)
790 ENDIF
791 END DO
792 END DO
793#endif
794#if defined BEDLOAD_VANDERA_DIRECT_UDELTA
795!
796! Find the near-bottom current velocity directly at a given
797! elevation (doesnot require Madsen output)
798! Find the angle at that near-bottom current velocity.
799!
800 DO j=jstr,jend
801 DO i=istr,iend
802 dstp=z_r(i,j,n(ng))-z_w(i,j,0)
803 cff=min( 0.98_r8*dstp, sg_zwbl(ng) )
804 DO k=1,n(ng)
805 urz(k)=0.5_r8*(u(i,j,k,nrhs)+u(i+1,j,k,nrhs))
806 vrz(k)=0.5_r8*(v(i,j,k,nrhs)+v(i,j+1,k,nrhs))
807# ifdef SSW_LOGINT_STOKES
808 urz(k)=urz(k)+0.5_r8*(u_stokes(i,j,k)+u_stokes(i+1,j,k))
809 vrz(k)=vrz(k)+0.5_r8*(v_stokes(i,j,k)+v_stokes(i,j+1,k))
810# endif
811 END DO
812 CALL log_interp( n(ng), dstp, cff, &
813 & urz, vrz, &
814 & z_r(i,j,:), z_w(i,j,:), &
815 & bottom(i,j,isd50), bottom(i,j,izapp), &
816 & zr_wbl(i,j), &
817 & ur_sgwbl(i,j), vr_sgwbl(i,j) )
818!
819! Compute bottom current magnitude at RHO-points.
820!
821 ucur_sgwbl=ur_sgwbl(i,j)
822 vcur_sgwbl=vr_sgwbl(i,j)
823 udelta_wbl(i,j)=sqrt(ur_sgwbl(i,j)*ur_sgwbl(i,j)+ &
824 & vr_sgwbl(i,j)*vr_sgwbl(i,j)+eps)
825!
826 IF (ucur_sgwbl.eq.0.0_r8) THEN
827 phic_sgwbl(i,j)=0.5_r8*pi*sign(1.0_r8,vcur_sgwbl)
828 ELSE
829 phic_sgwbl(i,j)=atan2(vcur_sgwbl,ucur_sgwbl)
830 ENDIF
831 END DO
832 END DO
833!
834#endif
835!
836!-----------------------------------------------------------------------
837! Compute kinematic bottom stress components due current and wind-
838! induced waves.
839!-----------------------------------------------------------------------
840!
841 DO j=jstr,jend
842 DO i=istru,iend
843 anglec=0.5_r8*(ur_sg(i,j)+ur_sg(i-1,j))/ &
844 & (0.5_r8*(umag(i-1,j)+umag(i,j)))
845 bustr(i,j)=0.5_r8*(tauc(i-1,j)+tauc(i,j))*anglec
846#ifdef WET_DRY
847 cff2=0.75_r8*0.5_r8*(z_w(i-1,j,1)+z_w(i,j,1)- &
848 & z_w(i-1,j,0)-z_w(i,j,0))
849 bustr(i,j)=sign(1.0_r8,bustr(i,j))*min(abs(bustr(i,j)), &
850 & abs(u(i,j,1,nrhs))*cff2/dt(ng))
851#endif
852 END DO
853 END DO
854 DO j=jstrv,jend
855 DO i=istr,iend
856 anglec=0.5_r8*(vr_sg(i,j)+vr_sg(i,j-1))/ &
857 & (0.5_r8*(umag(i,j-1)+umag(i,j)))
858 bvstr(i,j)=0.5_r8*(tauc(i,j-1)+tauc(i,j))*anglec
859#ifdef WET_DRY
860 cff2=0.75_r8*0.5_r8*(z_w(i,j-1,1)+z_w(i,j,1)- &
861 & z_w(i,j-1,0)-z_w(i,j,0))
862 bvstr(i,j)=sign(1.0_r8,bvstr(i,j))*min(abs(bvstr(i,j)), &
863 & abs(v(i,j,1,nrhs))*cff2/dt(ng))
864#endif
865 END DO
866 END DO
867 DO j=jstr,jend
868 DO i=istr,iend
869 anglec=ucur(i,j)/umag(i,j)
870 anglew=cos(1.5_r8*pi-dwave(i,j)-angler(i,j))
871 bustrc(i,j)=tauc(i,j)*anglec
872 bustrw(i,j)=tauw(i,j)*anglew
873 bustrcwmax(i,j)=taucwmax(i,j)*anglew
874 ubot(i,j)=ub(i,j)*anglew
875 ur(i,j)=ucur(i,j)
876!
877 anglec=vcur(i,j)/umag(i,j)
878 anglew=sin(1.5_r8*pi-dwave(i,j)-angler(i,j))
879 bvstrc(i,j)=tauc(i,j)*anglec
880 bvstrw(i,j)=tauw(i,j)*anglew
881 bvstrcwmax(i,j)=taucwmax(i,j)*anglew
882 vbot(i,j)=ub(i,j)*anglew
883 vr(i,j)=vcur(i,j)
884!
885 bottom(i,j,irlen)=rlength(i,j)
886 bottom(i,j,irhgt)=rheight(i,j)
887 bottom(i,j,ibwav)=ab(i,j)
888 bottom(i,j,izdef)=zodef(i,j)
889 bottom(i,j,izapp)=znotc(i,j)
890 bottom(i,j,iznik)=zon(i,j)
891 bottom(i,j,izbio)=zobio(i,j)
892 bottom(i,j,izbfm)=zobf(i,j)
893 bottom(i,j,izbld)=zost(i,j)
894 bottom(i,j,izwbl)=znot(i,j)
895 bottom(i,j,idtbl)=thck_wbl(i,j)
896#if defined BEDLOAD_VANDERA_MADSEN_UDELTA
897 bottom(i,j,idksd)=ksd_wbl(i,j)
898 bottom(i,j,idusc)=ustrc_wbl(i,j)
899#endif
900#if defined BEDLOAD_VANDERA_MADSEN_UDELTA || \
901 defined bedload_vandera_direct_udelta
902 bottom(i,j,idubl)=udelta_wbl(i,j)
903 bottom(i,j,idzrw)=zr_wbl(i,j)
904 bottom(i,j,idpcx)=phic_sgwbl(i,j)
905#else
906 bottom(i,j,idpcx)=phic(i,j)
907 bottom(i,j,idpwc)=phicw(i,j)
908#endif
909 END DO
910 END DO
911!
912! Apply periodic or gradient boundary conditions for output
913! purposes only.
914!
915 CALL bc_u2d_tile (ng, tile, &
916 & lbi, ubi, lbj, ubj, &
917 & bustr)
918 CALL bc_v2d_tile (ng, tile, &
919 & lbi, ubi, lbj, ubj, &
920 & bvstr)
921 CALL bc_r2d_tile (ng, tile, &
922 & lbi, ubi, lbj, ubj, &
923 & bustrc)
924 CALL bc_r2d_tile (ng, tile, &
925 & lbi, ubi, lbj, ubj, &
926 & bvstrc)
927 CALL bc_r2d_tile (ng, tile, &
928 & lbi, ubi, lbj, ubj, &
929 & bustrw)
930 CALL bc_r2d_tile (ng, tile, &
931 & lbi, ubi, lbj, ubj, &
932 & bvstrw)
933 CALL bc_r2d_tile (ng, tile, &
934 & lbi, ubi, lbj, ubj, &
935 & bustrcwmax)
936 CALL bc_r2d_tile (ng, tile, &
937 & lbi, ubi, lbj, ubj, &
938 & bvstrcwmax)
939 CALL bc_r2d_tile (ng, tile, &
940 & lbi, ubi, lbj, ubj, &
941 & ubot)
942 CALL bc_r2d_tile (ng, tile, &
943 & lbi, ubi, lbj, ubj, &
944 & vbot)
945 CALL bc_r2d_tile (ng, tile, &
946 & lbi, ubi, lbj, ubj, &
947 & ur)
948 CALL bc_r2d_tile (ng, tile, &
949 & lbi, ubi, lbj, ubj, &
950 & vr)
951 CALL bc_r2d_tile (ng, tile, &
952 & lbi, ubi, lbj, ubj, &
953 & bottom(:,:,irlen))
954 CALL bc_r2d_tile (ng, tile, &
955 & lbi, ubi, lbj, ubj, &
956 & bottom(:,:,irhgt))
957 CALL bc_r2d_tile (ng, tile, &
958 & lbi, ubi, lbj, ubj, &
959 & bottom(:,:,ibwav))
960 CALL bc_r2d_tile (ng, tile, &
961 & lbi, ubi, lbj, ubj, &
962 & bottom(:,:,izdef))
963 CALL bc_r2d_tile (ng, tile, &
964 & lbi, ubi, lbj, ubj, &
965 & bottom(:,:,izapp))
966 CALL bc_r2d_tile (ng, tile, &
967 & lbi, ubi, lbj, ubj, &
968 & bottom(:,:,iznik))
969 CALL bc_r2d_tile (ng, tile, &
970 & lbi, ubi, lbj, ubj, &
971 & bottom(:,:,izbio))
972 CALL bc_r2d_tile (ng, tile, &
973 & lbi, ubi, lbj, ubj, &
974 & bottom(:,:,izbfm))
975 CALL bc_r2d_tile (ng, tile, &
976 & lbi, ubi, lbj, ubj, &
977 & bottom(:,:,izbld))
978 CALL bc_r2d_tile (ng, tile, &
979 & lbi, ubi, lbj, ubj, &
980 & bottom(:,:,izwbl))
981#if defined BEDLOAD_VANDERA_MADSEN_UDELTA
982 CALL bc_r2d_tile (ng, tile, &
983 & lbi, ubi, lbj, ubj, &
984 & bottom(:,:,idtbl))
985 CALL bc_r2d_tile (ng, tile, &
986 & lbi, ubi, lbj, ubj, &
987 & bottom(:,:,idksd))
988 CALL bc_r2d_tile (ng, tile, &
989 & lbi, ubi, lbj, ubj, &
990 & bottom(:,:,idusc))
991#endif
992#if defined BEDLOAD_VANDERA_MADSEN_UDELTA || \
993 defined bedload_vandera_direct_udelta
994 CALL bc_r2d_tile (ng, tile, &
995 & lbi, ubi, lbj, ubj, &
996 & bottom(:,:,idubl))
997 CALL bc_r2d_tile (ng, tile, &
998 & lbi, ubi, lbj, ubj, &
999 & bottom(:,:,idzrw))
1000#else
1001 CALL bc_r2d_tile (ng, tile, &
1002 & lbi, ubi, lbj, ubj, &
1003 & bottom(:,:,idpwc))
1004#endif
1005 CALL bc_r2d_tile (ng, tile, &
1006 & lbi, ubi, lbj, ubj, &
1007 & bottom(:,:,idpcx))
1008#ifdef DISTRIBUTE
1009 CALL mp_exchange2d (ng, tile, inlm, 4, &
1010 & lbi, ubi, lbj, ubj, &
1011 & nghostpoints, &
1012 & ewperiodic(ng), nsperiodic(ng), &
1013 & bustr, bvstr, bustrc, bvstrc)
1014 CALL mp_exchange2d (ng, tile, inlm, 4, &
1015 & lbi, ubi, lbj, ubj, &
1016 & nghostpoints, &
1017 & ewperiodic(ng), nsperiodic(ng), &
1018 & bustrw, bvstrw, bustrcwmax, bvstrcwmax)
1019 CALL mp_exchange2d (ng, tile, inlm, 4, &
1020 & lbi, ubi, lbj, ubj, &
1021 & nghostpoints, &
1022 & ewperiodic(ng), nsperiodic(ng), &
1023 & ubot, vbot, ur, vr)
1024 CALL mp_exchange2d (ng, tile, inlm, 3, &
1025 & lbi, ubi, lbj, ubj, &
1026 & nghostpoints, &
1027 & ewperiodic(ng), nsperiodic(ng), &
1028 & bottom(:,:,irlen), &
1029 & bottom(:,:,irhgt), &
1030 & bottom(:,:,ibwav))
1031 CALL mp_exchange2d (ng, tile, inlm, 3, &
1032 & lbi, ubi, lbj, ubj, &
1033 & nghostpoints, &
1034 & ewperiodic(ng), nsperiodic(ng), &
1035 & bottom(:,:,izdef), &
1036 & bottom(:,:,izapp), &
1037 & bottom(:,:,iznik))
1038 CALL mp_exchange2d (ng, tile, inlm, 4, &
1039 & lbi, ubi, lbj, ubj, &
1040 & nghostpoints, &
1041 & ewperiodic(ng), nsperiodic(ng), &
1042 & bottom(:,:,izbio), &
1043 & bottom(:,:,izbfm), &
1044 & bottom(:,:,izbld), &
1045 & bottom(:,:,izwbl))
1046# if defined BEDLOAD_VANDERA_MADSEN_UDELTA
1047 CALL mp_exchange2d (ng, tile, inlm, 3, &
1048 & lbi, ubi, lbj, ubj, &
1049 & nghostpoints, &
1050 & ewperiodic(ng), nsperiodic(ng), &
1051 & bottom(:,:,idtbl), &
1052 & bottom(:,:,idksd), &
1053 & bottom(:,:,idusc))
1054# endif
1055# if defined BEDLOAD_VANDERA_MADSEN_UDELTA || \
1056 defined bedload_vandera_direct_udelta
1057 CALL mp_exchange2d (ng, tile, inlm, 2, &
1058 & lbi, ubi, lbj, ubj, &
1059 & nghostpoints, &
1060 & ewperiodic(ng), nsperiodic(ng), &
1061 & bottom(:,:,idzrw), &
1062 & bottom(:,:,idubl))
1063# else
1064 CALL mp_exchange2d (ng, tile, inlm, 1, &
1065 & lbi, ubi, lbj, ubj, &
1066 & nghostpoints, &
1067 & ewperiodic(ng), nsperiodic(ng), &
1068 & bottom(:,:,idpwc))
1069# endif
1070 CALL mp_exchange2d (ng, tile, inlm, 1, &
1071 & lbi, ubi, lbj, ubj, &
1072 & nghostpoints, &
1073 & ewperiodic(ng), nsperiodic(ng), &
1074 & bottom(:,:,idpcx))
1075#endif
1076!
1077 RETURN
1078 END SUBROUTINE ssw_bbl_tile
1079
1080#ifdef SGWC
1081!
1082 SUBROUTINE sg_bstress (sg_row, sg_zrozn, sg_phicw, sg_ubokur, &
1083 & sg_ubouc, sg_mu, sg_epsilon, sg_ro, &
1084 & sg_fofx)
1085!
1086!=======================================================================
1087! !
1088! This routine computes bottom stresses via bottom boundary layer !
1089! formulation of Styles and Glenn (1999). !
1090! !
1091! On Input: !
1092! !
1093! sg_row Ratio of wave excursion amplitude over roughness. !
1094! sg_zrozn Ratio of height of current over roughness. !
1095! sg_phiwc Angle between wave and currents (radians). !
1096! sg_ubokur Ratio of wave over current velocity: !
1097! ub/(vonKar*ur) !
1098! sg_ubouc Ratio of bed wave orbital over bottom shear stress !
1099! (ub/ustarc), first guess. !
1100! !
1101! On Output: !
1102! !
1103! sg_ubouc Ratio of bed wave orbital over bottom shear stress !
1104! (ub/ustarc), iterated value. !
1105! sg_mu Ratio between wave and current bottom shear !
1106! stresses (ustarwm/ustarc). !
1107! sg_epsilon Ratio between combined (wave and current) and !
1108! current bottom shear stresses (ustarc/ustarcw). !
1109! sg_ro Internal friction Rossby number: !
1110! ustarc/(omega*znot) !
1111! sg_fofx Root of PDE used for convergence. !
1112! !
1113!=======================================================================
1114!
1115 USE mod_param
1116 USE mod_scalars
1117!
1118! Imported variable declarations.
1119!
1120 real(r8), intent(in) :: sg_row, sg_zrozn, sg_phicw, sg_ubokur
1121
1122 real(r8), intent(inout) :: sg_ubouc
1123
1124 real(r8), intent(out) :: sg_mu, sg_epsilon, sg_ro, sg_fofx
1125!
1126! Local variable declarations.
1127!
1128 logical :: ITERATE
1129
1130 integer :: Iter
1131
1132 real(r8) :: cff, sg_bei, sg_beip, sg_ber, sg_berp, sg_cosphi
1133 real(r8) :: sg_eps2, sg_kei, sg_keip, sg_ker, sg_kerp, sg_mu2
1134 real(r8) :: sg_phi, sg_ror, sg_x, sg_z2p, sg_znotp, sg_zroz1
1135 real(r8) :: sg_zroz2, sg_z1ozn, sg_z2ozn
1136
1137 complex(c8) :: sg_argi, sg_bnot, sg_bnotp, sg_b1, sg_b1p
1138 complex(c8) :: sg_gammai, sg_knot, sg_knotp, sg_k1, sg_k1p
1139 complex(c8) :: sg_ll, sg_nn
1140!
1141!-----------------------------------------------------------------------
1142! Compute bottom stresses.
1143!-----------------------------------------------------------------------
1144!
1145! Compute nondimensional bottom wave shear, phi. Iterate to make
1146! sure that there is an upper limit in "ubouc". It usually requires
1147! only one pass.
1148!
1149 iterate=.true.
1150 DO iter=1,sg_n
1151 IF (iterate) THEN
1152 sg_ro=sg_row/sg_ubouc
1153 sg_znotp=1.0_r8/(sg_kappa*sg_ro)
1154 IF ((sg_z1p/sg_znotp).gt.1.0_r8) THEN
1155 sg_x=2.0_r8*sqrt(sg_znotp)
1156 IF (sg_x.le.8.0_r8) THEN
1157 CALL sg_kelvin8m (sg_x, sg_ber, sg_bei, sg_ker, sg_kei, &
1158 & sg_berp, sg_beip, sg_kerp, sg_keip)
1159 ELSE
1160 CALL sg_kelvin8p (sg_x, sg_ker, sg_kei, sg_ber, sg_bei, &
1161 & sg_kerp, sg_keip, sg_berp, sg_beip)
1162 END IF
1163 cff=1.0_r8/sqrt(sg_znotp)
1164 sg_bnot =cmplx(sg_ber,sg_bei,8)
1165 sg_knot =cmplx(sg_ker,sg_kei,8)
1166 sg_bnotp=cmplx(sg_berp,sg_beip,8)*cff
1167 sg_knotp=cmplx(sg_kerp,sg_keip,8)*cff
1168!
1169 sg_x=2.0_r8*sqrt(sg_z1p)
1170 IF (sg_x.le.8.0_r8) THEN
1171 CALL sg_kelvin8m (sg_x, sg_ber, sg_bei, sg_ker, sg_kei, &
1172 & sg_berp, sg_beip, sg_kerp, sg_keip)
1173 ELSE
1174 CALL sg_kelvin8p (sg_x, sg_ker, sg_kei, sg_ber, sg_bei, &
1175 & sg_kerp, sg_keip, sg_berp, sg_beip)
1176 END IF
1177 cff=1.0_r8/sqrt(sg_z1p)
1178 sg_b1 =cmplx(sg_ber,sg_bei,8)
1179 sg_k1 =cmplx(sg_ker,sg_kei,8)
1180 sg_b1p=cmplx(sg_berp,sg_beip,8)*cff
1181 sg_k1p=cmplx(sg_kerp,sg_keip,8)*cff
1182!
1183 sg_ll=sg_mp*sg_b1+sg_b1p
1184 sg_nn=sg_mp*sg_k1+sg_k1p
1185 sg_argi=sg_bnotp*sg_nn/(sg_bnot*sg_nn-sg_knot*sg_ll)+ &
1186 & sg_knotp*sg_ll/(sg_knot*sg_ll-sg_bnot*sg_nn)
1187 sg_gammai=-sg_kappa*sg_znotp*sg_argi
1188 sg_phi=cabs(sg_gammai)
1189 ELSE
1190 sg_gammai=-sg_kappa*sg_z1p*sg_mp
1191 sg_phi=cabs(sg_gammai)
1192 END IF
1193!
1194 IF (sg_ubouc.gt.(1.0_r8/sg_phi)) THEN
1195 sg_ubouc=1.0_r8/sg_phi
1196 ELSE
1197 iterate=.false.
1198 END IF
1199 END IF
1200 END DO
1201!
1202! Compute ratio of wave over current bottom shear stresses.
1203!
1204 sg_mu=sqrt(sg_ubouc*sg_phi)
1205!
1206! Compute ratio of current over combined bottom shear stresses.
1207!
1208 IF (sg_mu.eq.1.0_r8) THEN
1209 sg_epsilon=0.0_r8
1210 ELSE
1211 sg_mu2=sg_mu*sg_mu
1212 sg_cosphi=abs(cos(sg_phicw))
1213 sg_eps2=-sg_mu2*sg_cosphi+ &
1214 & sqrt(1.0_r8+sg_mu2*sg_mu2*(sg_cosphi*sg_cosphi-1.0_r8))
1215 sg_epsilon=sqrt(sg_eps2)
1216 END IF
1217!
1218! Determine root of PDE used for convergence.
1219!
1220 IF (sg_epsilon.ne.0.0_r8) THEN
1221 sg_z2p=sg_z1p/sg_epsilon
1222 sg_ror=sg_ro/sg_zrozn
1223 sg_zroz1=1.0_r8/(sg_alpha*sg_kappa*sg_ror)
1224 sg_zroz2=sg_epsilon*sg_zroz1
1225 sg_z1ozn=sg_alpha*sg_kappa*sg_ro
1226 sg_z2ozn=sg_z1ozn/sg_epsilon
1227!
1228 IF ((sg_zroz2.gt.1.0_r8).and.(sg_z1ozn.gt.1.0_r8)) THEN
1229 sg_fofx=-sg_ubouc+sg_ubokur*sg_epsilon* &
1230 & (log(sg_zroz2)+1.0_r8-sg_epsilon+ &
1231 & sg_epsilon*log(sg_z1ozn))
1232 ELSE IF ((sg_zroz2.le.1.0_r8).and.(sg_zroz1.gt.1.0_r8).and. &
1233 & (sg_z1ozn.gt.1.0_r8)) THEN
1234 sg_fofx=-sg_ubouc+sg_ubokur*sg_epsilon*sg_epsilon* &
1235 & (sg_zroz1-1.0_r8+log(sg_z1ozn))
1236 ELSE IF ((sg_zroz1.le.1.0_r8).and.(sg_z1ozn.gt.1.0_r8)) THEN
1237 sg_fofx=-sg_ubouc+sg_ubokur*sg_epsilon*sg_epsilon* &
1238 & log(sg_zrozn)
1239 ELSE IF ((sg_zroz2.gt.1.0_r8).and.(sg_z1ozn.le.1.0_r8).and. &
1240 & (sg_z2ozn.gt.1.0_r8)) THEN
1241 sg_fofx=-sg_ubouc+sg_ubokur*sg_epsilon* &
1242 & (log(sg_zroz2)+1.0_r8-1.0_r8/sg_z2ozn)
1243 ELSE IF ((sg_zroz2.le.1.0_r8).and.(sg_zroz1.gt.1.0_r8).and. &
1244 & (sg_z1ozn.le.1.0_r8).and.(sg_z2ozn.gt.1.0_r8)) THEN
1245 sg_fofx=-sg_ubouc+sg_ubokur*sg_epsilon*sg_epsilon* &
1246 & (sg_zroz1-1.0_r8/sg_z1ozn)
1247 ELSE IF ((sg_zroz2.gt.1.0_r8).and.(sg_z2ozn.le.1.0_r8)) THEN
1248 sg_fofx=-sg_ubouc+sg_ubokur*sg_epsilon*log(sg_zrozn)
1249 END IF
1250 END IF
1251!
1252 RETURN
1253 END SUBROUTINE sg_bstress
1254!
1255 SUBROUTINE sg_purewave (sg_row, sg_ubouwm, sg_znotp, sg_ro)
1256!
1257!=======================================================================
1258! !
1259! This routine determines the maximum ratio of waves over combined !
1260! bottom shear stress. !
1261! !
1262! On Input: !
1263! !
1264! sg_row Ratio of wave excursion amplitude over roughness. !
1265! sg_ubouwm Maximum ratio of waves over combined bottom shear !
1266! stress. !
1267! !
1268! On Output: !
1269! !
1270! sg_ubouwm Maximum ratio of waves over combined bottom shear !
1271! stress. !
1272! sg_znotp Ratio of hydraulic roughness over scaled height !
1273! of bottom boundary layer. !
1274! sg_ro Internal friction Rossby number: !
1275! ustarc/(omega*znot) !
1276! !
1277!=======================================================================
1278!
1279 USE mod_param
1280 USE mod_scalars
1281!
1282! Imported variable declarations.
1283!
1284 real(r8), intent(in) :: sg_row
1285
1286 real(r8), intent(inout) :: sg_ubouwm
1287
1288 real(r8), intent(out) :: sg_znotp, sg_ro
1289!
1290! Local variable declarations.
1291!
1292 integer :: Iter
1293
1294 real(r8) :: cff, sg_bei, sg_beip, sg_ber, sg_berp, sg_kei
1295 real(r8) :: sg_keip, sg_ker, sg_kerp, sg_phi, sg_ubouwmn, sg_x
1296
1297 complex(c8) :: sg_argi, sg_bnot, sg_bnotp, sg_b1, sg_b1p
1298 complex(c8) :: sg_gammai, sg_knot, sg_knotp, sg_k1, sg_k1p
1299 complex(c8) :: sg_ll, sg_nn
1300!
1301!-----------------------------------------------------------------------
1302! Compute wind-induced wave stress.
1303!-----------------------------------------------------------------------
1304!
1305 DO iter=1,sg_n
1306 sg_ro=sg_row/sg_ubouwm
1307 sg_znotp=1.0_r8/(sg_kappa*sg_ro)
1308 IF (sg_z1p/sg_znotp.gt.1.0_r8) THEN
1309 sg_x=2.0_r8*sqrt(sg_znotp)
1310 IF (sg_x.le.8.0_r8) THEN
1311 CALL sg_kelvin8m (sg_x, sg_ber, sg_bei, sg_ker, sg_kei, &
1312 & sg_berp, sg_beip, sg_kerp, sg_keip)
1313 ELSE
1314 CALL sg_kelvin8p (sg_x, sg_ker, sg_kei, sg_ber, sg_bei, &
1315 & sg_kerp, sg_keip, sg_berp, sg_beip)
1316 END IF
1317 cff=1.0_r8/sqrt(sg_znotp)
1318 sg_bnot =cmplx(sg_ber,sg_bei,8)
1319 sg_knot =cmplx(sg_ker,sg_kei,8)
1320 sg_bnotp=cmplx(sg_berp,sg_beip,8)*cff
1321 sg_knotp=cmplx(sg_kerp,sg_keip,8)*cff
1322!
1323 sg_x=2.0*sqrt(sg_z1p)
1324 IF (sg_x.le.8.0_r8) THEN
1325 CALL sg_kelvin8m (sg_x, sg_ber, sg_bei, sg_ker, sg_kei, &
1326 & sg_berp, sg_beip, sg_kerp, sg_keip)
1327 ELSE
1328 CALL sg_kelvin8p (sg_x, sg_ker, sg_kei, sg_ber, sg_bei, &
1329 & sg_kerp, sg_keip, sg_berp, sg_beip)
1330 END IF
1331 cff=1.0_r8/sqrt(sg_z1p)
1332 sg_b1 =cmplx(sg_ber,sg_bei,8)
1333 sg_k1 =cmplx(sg_ker,sg_kei,8)
1334 sg_b1p=cmplx(sg_berp,sg_beip,8)*cff
1335 sg_k1p=cmplx(sg_kerp,sg_keip,8)*cff
1336!
1337 sg_ll=sg_mp*sg_b1+sg_b1p
1338 sg_nn=sg_mp*sg_k1+sg_k1p
1339 sg_argi=sg_bnotp*sg_nn/(sg_bnot*sg_nn-sg_knot*sg_ll)+ &
1340 & sg_knotp*sg_ll/(sg_knot*sg_ll-sg_bnot*sg_nn)
1341 sg_gammai=-sg_kappa*sg_znotp*sg_argi
1342 sg_phi=cabs(sg_gammai)
1343 ELSE
1344 sg_gammai=-sg_kappa*sg_z1p*sg_mp
1345 sg_phi=cabs(sg_gammai)
1346 END IF
1347!
1348 sg_ubouwmn=1.0_r8/sg_phi
1349 IF (abs((sg_ubouwmn-sg_ubouwm)/sg_ubouwmn).le.sg_tol) THEN
1350 sg_ubouwm=sg_ubouwmn
1351 RETURN
1352 ELSE
1353 sg_ubouwm=sg_ubouwmn
1354 END IF
1355 END DO
1356!
1357 RETURN
1358 END SUBROUTINE sg_purewave
1359!
1360 SUBROUTINE sg_kelvin8m (x, ber, bei, ker, kei, berp, beip, &
1361 & kerp, keip)
1362!
1363!=======================================================================
1364! !
1365! This rotuine computes the Kelvin functions for arguments less !
1366! than eight (p 384 Abram and Stegun). !
1367! !
1368!=======================================================================
1369!
1370 USE mod_scalars
1371!
1372! Imported variable declarations.
1373!
1374 real(r8), intent(in) :: x
1375 real(r8), intent(out) :: ber, bei, ker, kei
1376 real(r8), intent(out) :: berp, beip, kerp, keip
1377!
1378! Local variable declarations.
1379!
1380 integer :: i
1381
1382 real(r8) :: cff, xhalf
1383
1384 real(r8), dimension(28) :: xp
1385!
1386!-----------------------------------------------------------------------
1387! Compute Kelvin functions.
1388!-----------------------------------------------------------------------
1389!
1390 cff=0.125_r8*x
1391 xp(1)=cff
1392 DO i=2,28
1393 xp(i)=xp(i-1)*cff
1394 END DO
1395 xhalf=0.5_r8*x
1396!
1397 ber=1.0_r8- &
1398 & 64.0_r8*xp(4)+113.77777774_r8*xp(8)- &
1399 & 32.36345652_r8*xp(12)+2.64191397_r8*xp(16)- &
1400 & 0.08349609_r8*xp(20)+0.00122552_r8*xp(24)- &
1401 & 0.00000901_r8*xp(28)
1402 bei=16.0_r8*xp(2)-113.77777774_r8*xp(6)+ &
1403 & 72.81777742*xp(10)-10.56765779_r8*xp(14)+ &
1404 & 0.52185615_r8*xp(18)-0.01103667_r8*xp(22)+ &
1405 & 0.00011346*xp(26)
1406!
1407 ker=-ber*log(xhalf)+0.25_r8*pi*bei- &
1408 & 0.57721566_r8-59.05819744*xp(4)+ &
1409 & 171.36272133_r8*xp(8)-60.60977451_r8*xp(12)+ &
1410 & 5.65539121_r8*xp(16)-0.19636347_r8*xp(20)+ &
1411 & 0.00309699_r8*xp(24)-0.00002458_r8*xp(28)
1412 kei=-bei*log(xhalf)-0.25_r8*pi*ber+ &
1413 & 6.76454936_r8*xp(2)-142.91827687_r8*xp(6)+ &
1414 & 124.23569650_r8*xp(10)-21.30060904_r8*xp(14)+ &
1415 & 1.17509064_r8*xp(18)-0.02695875_r8*xp(22)+ &
1416 & 0.00029532_r8*xp(26)
1417!
1418 berp=x*(-4.0_r8*xp(2)+14.22222222_r8*xp(6)- &
1419 & 6.06814810_r8*xp(10)+0.66047849_r8*xp(14)- &
1420 & 0.02609253_r8*xp(18)+0.00045957_r8*xp(22)- &
1421 & 0.00000394_r8*xp(26))
1422 beip=x*(0.5_r8-10.66666666_r8*xp(4)+11.37777772_r8*xp(8)- &
1423 & 2.31167514_r8*xp(12)+0.14677204_r8*xp(16)- &
1424 & 0.00379386_r8*xp(20)+0.00004609_r8*xp(24))
1425!
1426 kerp=-berp*log(xhalf)-ber/x+0.25*pi*beip+ &
1427 & x*(-3.69113734_r8*xp(2)+21.42034017_r8*xp(6)- &
1428 & 11.36433272_r8*xp(10)+1.41384780_r8*xp(14)- &
1429 & 0.06136358_r8*xp(18)+0.00116137_r8*xp(22)- &
1430 & 0.00001075*xp(26))
1431 keip=-beip*log(xhalf)-bei/x-0.25_r8*pi*berp+ &
1432 & x*(0.21139217_r8-13.39858846_r8*xp(4)+ &
1433 & 19.41182758_r8*xp(8)-4.65950823_r8*xp(12)+ &
1434 & 0.33049424_r8*xp(16)-0.00926707_r8*xp(20)+ &
1435 & 0.00011997_r8*xp(24))
1436!
1437 RETURN
1438 END SUBROUTINE sg_kelvin8m
1439!
1440 SUBROUTINE sg_kelvin8p (x, ker, kei, ber, bei, kerp, keip, &
1441 & berp, beip)
1442!
1443!=======================================================================
1444! !
1445! This rotuine computes the Kelvin functions for arguments greater !
1446! than eight. !
1447! !
1448!=======================================================================
1449!
1450 USE mod_scalars
1451!
1452! Imported variable declarations.
1453!
1454 real(r8), intent(in) :: x
1455 real(r8), intent(out) :: ker, kei, ber, bei
1456 real(r8), intent(out) :: kerp, keip, berp, beip
1457!
1458! Local variable declarations.
1459!
1460 integer :: i
1461
1462 real(r8) :: cff, xhalf
1463
1464 real(r8), dimension(6) :: xm, xp
1465
1466 complex(c8) :: argm, argp, fofx, gofx, phim, phip, thetam, thetap
1467!
1468!-----------------------------------------------------------------------
1469! Compute Kelvin functions.
1470!-----------------------------------------------------------------------
1471!
1472 cff=8.0_r8/x
1473 xp(1)=cff
1474 xm(1)=-cff
1475 DO i=2,6
1476 xp(i)=xp(i-1)*cff
1477 xm(i)=-xm(i-1)*cff
1478 END DO
1479!
1480 thetap=cmplx(0.0_r8,-0.3926991_r8,8)+ &
1481 & cmplx(0.0110486_r8,-0.0110485_r8,8)*xp(1)+ &
1482 & cmplx(0.0_r8,-0.0009765_r8,8)*xp(2)+ &
1483 & cmplx(-0.0000906_r8,-0.0000901_r8,8)*xp(3)+ &
1484 & cmplx(-0.0000252_r8,0.0_r8,8)*xp(4)+ &
1485 & cmplx(-0.0000034_r8,0.0000051_r8,8)*xp(5)+ &
1486 & cmplx(0.0000006,0.0000019,8)*xp(6)
1487 thetam=cmplx(0.0_r8,-0.3926991_r8,8)+ &
1488 & cmplx(0.0110486_r8,-0.0110485_r8,8)*xm(1)+ &
1489 & cmplx(0.0_r8,-0.0009765_r8,8)*xm(2)+ &
1490 & cmplx(-0.0000906_r8,-0.0000901_r8,8)*xm(3)+ &
1491 & cmplx(-0.0000252_r8,0.0_r8,8)*xm(4)+ &
1492 & cmplx(-0.0000034_r8,0.0000051_r8,8)*xm(5)+ &
1493 & cmplx(0.0000006_r8,0.0000019_r8,8)*xm(6)
1494!
1495 phip=cmplx(0.7071068_r8,0.7071068_r8,8)+ &
1496 & cmplx(-0.0625001_r8,-0.0000001_r8,8)*xp(1)+ &
1497 & cmplx(-0.0013813_r8,0.0013811_r8,8)*xp(2)+ &
1498 & cmplx(0.0000005_r8,0.0002452_r8,8)*xp(3)+ &
1499 & cmplx(0.0000346_r8,0.0000338_r8,8)*xp(4)+ &
1500 & cmplx(0.0000117_r8,-0.0000024_r8,8)*xp(5)+ &
1501 & cmplx(0.0000016_r8,-0.0000032_r8,8)*xp(6)
1502 phim=cmplx(0.7071068_r8,0.7071068_r8,8)+ &
1503 & cmplx(-0.0625001_r8,-0.0000001_r8,8)*xm(1)+ &
1504 & cmplx(-0.0013813_r8,0.0013811_r8,8)*xm(2)+ &
1505 & cmplx(0.0000005_r8,0.0002452_r8,8)*xm(3)+ &
1506 & cmplx(0.0000346_r8,0.0000338_r8,8)*xm(4)+ &
1507 & cmplx(0.0000117_r8,-0.0000024_r8,8)*xm(5)+ &
1508 & cmplx(0.0000016_r8,-0.0000032_r8,8)*xm(6)
1509!
1510 cff=x/sqrt(2.0_r8)
1511 argm=-cff*cmplx(1.0_r8,1.0_r8,8)+thetam
1512 fofx=sqrt(pi/(2.0_r8*x))*cexp(argm)
1513 ker=real(fofx)
1514 kei=aimag(fofx)
1515!
1516 argp=cff*cmplx(1.0_r8,1.0_r8,8)+thetap
1517 gofx=1.0_r8/sqrt(2.0_r8*pi*x)*cexp(argp)
1518 ber=real(gofx)-kei/pi
1519 bei=aimag(gofx)+ker/pi
1520!
1521 kerp=real(-fofx*phim)
1522 keip=aimag(-fofx*phim)
1523!
1524 berp=real(gofx*phip)-keip/pi
1525 beip=aimag(gofx*phip)+kerp/pi
1526!
1527 RETURN
1528 END SUBROUTINE sg_kelvin8p
1529#endif
1530
1531#ifdef M94WC
1532 SUBROUTINE madsen94 (ubr, wr, ucr, zr, phiwc, kN, Dstp, &
1533 & ustrc, ustrwm, ustrr, fwc, zoa, dwc_va)
1534!
1535!=======================================================================
1536! !
1537! Grant-Madsen model from Madsen (1994). !
1538! !
1539! On Input: !
1540! !
1541! ubr Rep. wave-orbital velocity amplitude outside WBL (m/s). !
1542! wr Rep. angular wave frequency, 2* pi/T (rad/s). !
1543! ucr Current velocity at height zr (m/s). !
1544! zr Reference height for current velocity (m). !
1545! phiwc Angle between currents and waves at zr (radians). !
1546! kN Bottom roughness height, like Nikuradse k, (m). !
1547! Dstp Total water depth (m). !
1548! !
1549! On Output: !
1550! !
1551! ustrc Current friction velocity, u*c (m/s). !
1552! ustrwm Wave maximum friction velocity, u*wm (m/s). !
1553! ustrr Wave-current combined friction velocity, u*r (m/s). !
1554! fwc Wave friction factor (nondimensional). !
1555! zoa Apparent bottom roughness (m). !
1556! dwc_va Wave-boundary layer thickness (m). !
1557! !
1558!=======================================================================
1559!
1560 USE mod_param
1561 USE mod_scalars
1562!
1563! Imported variable declarations.
1564!
1565 real(r8), intent(in) :: ubr, wr, ucr, zr, phiwc, Dstp
1566 real(r8), intent(inout) :: kN
1567 real(r8), intent(out) :: ustrc, ustrwm, ustrr, fwc, zoa
1568 real(r8), intent(out) :: dwc_va
1569!
1570! Local variable declarations.
1571!
1572 integer, parameter :: MAXIT = 20
1573 integer :: i, nit
1574 integer :: iverbose = 1
1575
1576 real(r8) :: bigsqr, cosphiwc, cukw, diff, lndw, lnln, lnzr
1577 real(r8) :: phicwc, zo
1578 real(r8) :: dval = 99.99_r8
1579
1580 real(r8), dimension(MAXIT) :: Cmu
1581 real(r8), dimension(MAXIT) :: dwc
1582 real(r8), dimension(MAXIT) :: fwci
1583 real(r8), dimension(MAXIT) :: rmu
1584 real(r8), dimension(MAXIT) :: ustrci
1585 real(r8), dimension(MAXIT) :: ustrr2
1586 real(r8), dimension(MAXIT) :: ustrwm2
1587!
1588!-----------------------------------------------------------------------
1589! Compute bottom friction velocities and roughness.
1590!-----------------------------------------------------------------------
1591!
1592! Set special default values.
1593!
1594 ustrc=dval
1595 ustrwm=dval
1596 ustrr=dval
1597# ifdef CRS_FIX
1598 kn=min(kn,0.9_r8*zr)
1599# endif
1600 zoa=kn/30.0_r8
1601 phicwc=phiwc
1602
1603 zo = kn/30.0_r8
1604
1605 IF (ubr.le.0.01_r8) THEN
1606 dwc_va=kn
1607 zoa=dwc_va
1608 fwc=0.0_r8
1609 IF (ucr.le. 0.01_r8) THEN ! no waves or currents
1610 ustrc=0.0_r8
1611 ustrwm=0.0_r8
1612 ustrr=0.0_r8
1613 RETURN
1614 END IF
1615 ustrc=ucr*vonkar/log(zr/zo) ! no waves
1616 ustrwm=0.0_r8
1617 ustrr=ustrc
1618 RETURN
1619 END IF
1620!
1621! Iterate to compute friction velocities, roughness, and wave friction
1622! factor. Notice that the computation of the wave friction factor
1623! has been inlined for efficiency.
1624!
1625 cosphiwc=abs(cos(phiwc))
1626 rmu(1)=0.0_r8
1627 cmu(1)=1.0_r8
1628
1629 cukw=cmu(1)*ubr/(kn*wr)
1630# if defined CRS_FIX
1631!
1632! New fwc CRS calculation
1633!
1634 fwci(1)=cmu(1)*0.3_r8
1635 IF ((cukw.gt.0.352_r8).and.(cukw.le.100.0_r8)) THEN ! Eq 32/33
1636 fwci(1)=cmu(1)*exp(7.02_r8*cukw**(-0.078_r8)-8.82_r8)
1637 ELSE IF (cukw.gt.100.0_r8) THEN
1638 fwci(1)=cmu(1)*exp(5.61_r8*cukw**(-0.109_r8)-7.30_r8)
1639 END IF
1640# else
1641!
1642! Original method fwc calculation
1643!
1644 IF ((cukw.gt.0.2_r8).and.(cukw.le.100.0_r8)) THEN ! Eq 32/33
1645 fwci(1)=cmu(1)*exp(7.02_r8*cukw**(-0.078_r8)-8.82_r8)
1646 ELSE IF ((cukw.gt.100.).and.(cukw.le.10000.0_r8)) THEN
1647 fwci(1)=cmu(1)*exp(5.61_r8*cukw**(-0.109_r8)-7.30_r8)
1648 ELSE IF (cukw.gt.10000.0_r8 ) THEN
1649 fwci(1)=cmu(1)*exp(5.61_r8*10000.0_r8**(-0.109_r8)-7.30_r8)
1650 ELSE
1651 fwci(1)=cmu(1)*0.43_r8
1652 END IF
1653# endif
1654!
1655 ustrwm2(1)=0.5_r8*fwci(1)*ubr*ubr ! Eq 29
1656 ustrr2(1)=cmu(1)*ustrwm2(1) ! Eq 26
1657 ustrr=sqrt(ustrr2(1))
1658 IF (cukw.ge.8.0_r8) THEN
1659 dwc(1)=2.0_r8*vonkar*ustrr/wr
1660# if defined CRS_FIX
1661 dwc(1)=min( 0.9_r8*zr, dwc(1) )
1662# endif
1663 ELSE
1664 dwc(1)=kn
1665 END IF
1666 lnzr=log(zr/dwc(1))
1667 lndw=log(dwc(1)/zo)
1668 lnln=lnzr/lndw
1669 bigsqr=-1.0_r8+sqrt(1.0_r8+((4.0_r8*vonkar*lndw)/ &
1670 & (lnzr*lnzr))*ucr/ustrr)
1671 ustrci(1)=0.5_r8*ustrr*lnln*bigsqr
1672!
1673 i=1
1674 diff=1.0_r8
1675 DO WHILE ((i.lt.maxit).and.(diff.gt.0.000005_r8))
1676 i=i+1
1677 rmu(i)=ustrci(i-1)*ustrci(i-1)/ustrwm2(i-1)
1678 cmu(i)=sqrt(1.0_r8+ &
1679 & 2.0_r8*rmu(i)*cosphiwc+rmu(i)*rmu(i)) ! Eq 27
1680 cukw=cmu(i)*ubr/(kn*wr)
1681# ifdef CRS_FIX
1682!
1683! New fwc CRS calculation
1684!
1685 fwci(i)=cmu(i)*0.3_r8
1686 IF ((cukw.gt.0.352_r8).and.(cukw.le.100.0_r8)) THEN ! Eq 32/33
1687 fwci(i)=cmu(i)*exp(7.02_r8*cukw**(-0.078_r8)-8.82_r8)
1688 ELSE IF (cukw.gt.100.0_r8) THEN
1689 fwci(i)=cmu(i)*exp(5.61_r8*cukw**(-0.109_r8)-7.30_r8)
1690 END IF
1691# else
1692!
1693! Original method fwc calculation
1694!
1695 IF ((cukw.gt.0.2_r8).and.(cukw.le.100.0_r8)) THEN ! Eq 32/33
1696 fwci(i)=cmu(i)*exp(7.02_r8*cukw**(-0.078_r8)-8.82_r8)
1697 ELSE IF ((cukw.gt.100.).and.(cukw.le.10000.0_r8)) THEN
1698 fwci(i)=cmu(i)*exp(5.61_r8*cukw**(-0.109_r8)-7.30_r8)
1699 ELSE IF (cukw.gt.10000.0_r8 ) THEN
1700 fwci(i)=cmu(i)*exp(5.61_r8*10000.0_r8**(-0.109_r8)-7.30_r8)
1701 ELSE
1702 fwci(i)=cmu(i)*0.43_r8
1703 END IF
1704# endif
1705!
1706 ustrwm2(i)=0.5_r8*fwci(i)*ubr*ubr ! Eq 29
1707 ustrr2(i)=cmu(i)*ustrwm2(i) ! Eq 26
1708 ustrr=sqrt(ustrr2(i))
1709!! IF ((Cmu(1)*ubr/(kN*wr)).ge.8.0_r8) THEN ! HGA Why 1?
1710 IF (cukw.ge.8.0_r8) THEN
1711 dwc(i)=2.0_r8*vonkar*ustrr/wr ! Eq 36
1712# if defined CRS_FIX
1713 dwc(i)=min( 0.9_r8*zr, dwc(i) )
1714# endif
1715 ELSE
1716 dwc(i)=kn
1717 END IF
1718 lnzr=log(zr/dwc(i))
1719 lndw=log(dwc(i)/zo)
1720 lnln=lnzr/lndw
1721 bigsqr=-1.0_r8+sqrt(1.0_r8+((4.0_r8*vonkar*lndw)/ &
1722 & (lnzr*lnzr))*ucr/ustrr)
1723 ustrci(i)=0.5_r8*ustrr*lnln*bigsqr ! Eq 38
1724 diff=abs((fwci(i)-fwci(i-1))/fwci(i))
1725 END DO
1726 ustrwm=sqrt(ustrwm2(i))
1727 ustrc=ustrci(i)
1728 ustrr=sqrt(ustrr2(i))
1729 phicwc=phiwc
1730 zoa=exp(log(dwc(i))-(ustrc/ustrr)*log(dwc(i)/zo)) ! Eq 11
1731 dwc_va=dwc(i)
1732 fwc=fwci(i)
1733!
1734 RETURN
1735 END SUBROUTINE madsen94
1736#endif
1737
1738#if defined SSW_LOGINT || defined BEDLOAD_VANDERA_MADSEN_UDELTA
1739!
1740 SUBROUTINE log_interp (kmax, Dstp, sg_loc, u_1d, v_1d, &
1741 & z_r_1d, z_w_1d, &
1742 & d50, zapp_loc, &
1743 & Zr_sg, &
1744 & Ur_sg, Vr_sg)
1745!
1746!=======================================================================
1747! !
1748! Find the near-bottom current velocity in x- and y- directions !
1749! that corresponds to user input elevation to get near-bottom !
1750! current velocity (m) !
1751! !
1752!=======================================================================
1753!
1754 USE mod_param
1755 USE mod_sediment
1756 USE mod_scalars
1757!
1758 implicit none
1759!
1760! Imported variable declarations.
1761!
1762 integer, intent(in) :: kmax
1763 real(r8), intent(in) :: Dstp, sg_loc
1764 real(r8), intent(in) :: u_1d(1:kmax), v_1d(1:kmax)
1765 real(r8), intent(in) :: z_r_1d(1:kmax), z_w_1d(0:kmax)
1766 real(r8), intent(in) :: d50, zapp_loc
1767 real(r8), intent(out) :: Zr_sg
1768 real(r8), intent(out) :: Ur_sg, Vr_sg
1769!
1770! Local variables.
1771!
1772 integer :: k
1773 real(r8) :: z1, z2, Zr
1774 real(r8) :: fac, fac1, fac2
1775!
1776 zr=z_r_1d(1)-z_w_1d(0)
1777!
1778 IF ( sg_loc.ge.zr ) THEN
1779!
1780! If chosen height to get near bottom-current velocity lies
1781! within any vertical level, perform logarithmic interpolation.
1782!
1783 DO k=2,kmax
1784 z1=z_r_1d(k-1)-z_w_1d(0)
1785 z2=z_r_1d(k )-z_w_1d(0)
1786 IF ( ( z1.le.sg_loc ).and.( sg_loc.lt.z2 )) THEN
1787 fac=1.0_r8/log(z2/z1)
1788 fac1=fac*log(z2/sg_loc)
1789 fac2=fac*log(sg_loc/z1)
1790!
1791 ur_sg=fac1*u_1d(k-1)+fac2*u_1d(k)
1792 vr_sg=fac1*v_1d(k-1)+fac2*v_1d(k)
1793 zr_sg=sg_loc
1794 ENDIF
1795 IF ((k.eq.kmax).and.(sg_loc.ge.z2)) THEN
1796 ur_sg=u_1d(k)
1797 vr_sg=v_1d(k)
1798 zr_sg=z2
1799 END IF
1800 END DO
1801 ELSE !IF ( sg_loc.lt.Zr ) THEN
1802 z1=max( 2.5_r8*d50/30.0_r8, zapp_loc, 1.0e-10_r8 )
1803 z2=zr
1804!
1805 IF ( sg_loc.lt.z1 ) THEN
1806!
1807! If chosen height is less than the bottom roughness
1808! perform linear interpolation.
1809!
1810 z1=sg_loc
1811 fac=z1/z2
1812!
1813 ur_sg=fac*u_1d(1)
1814 vr_sg=fac*v_1d(1)
1815 zr_sg=sg_loc
1816!
1817 ELSE !IF ( sg_loc.gt.z1 ) THEN
1818!
1819! If chosen height is less than the bottom cell thickness
1820! perform logarithmic interpolation with bottom roughness.
1821!
1822 fac=1.0_r8/log(z2/z1)
1823 fac2=fac*log(sg_loc/z1)
1824!
1825 ur_sg=fac2*u_1d(1)
1826 vr_sg=fac2*v_1d(1)
1827 zr_sg=sg_loc
1828 END IF
1829 END IF
1830!
1831 RETURN
1832 END SUBROUTINE log_interp
1833#endif
1834 END MODULE bbl_mod
bbl_mod
Definition mb_bbl.h:1
bbl_mod::sg_purewave
subroutine sg_purewave(sg_row, sg_ubouwm, sg_znotp, sg_ro)
Definition sg_bbl.h:875
bbl_mod::log_interp
subroutine log_interp(kmax, dstp, sg_loc, u_1d, v_1d, z_r_1d, z_w_1d, d50, zapp_loc, zr_sg, ur_sg, vr_sg)
Definition ssw_bbl.h:1745
bbl_mod::bblm
subroutine, public bblm(ng, tile)
Definition mb_bbl.h:55
bbl_mod::sg_kelvin8p
subroutine sg_kelvin8p(x, ker, kei, ber, bei, kerp, keip, berp, beip)
Definition sg_bbl.h:1061
bbl_mod::ssw_bbl_tile
subroutine ssw_bbl_tile(ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, nrhs, h, z_r, z_w, angler, zobot, hwave, uwave_rms, dwave, pwave_bot, bedldu, bedldv, bottom, rho, u, v, u_stokes, v_stokes, zeta, iconv, ubot, vbot, ur, vr, bustrc, bvstrc, bustrw, bvstrw, bustrcwmax, bvstrcwmax, bustr, bvstr)
Definition ssw_bbl.h:141
bbl_mod::sg_bstress
subroutine sg_bstress(sg_row, sg_zrozn, sg_phicw, sg_ubokur, sg_ubouc, sg_mu, sg_epsilon, sg_ro, sg_fofx)
Definition sg_bbl.h:704
bbl_mod::sg_kelvin8m
subroutine sg_kelvin8m(x, ber, bei, ker, kei, berp, beip, kerp, keip)
Definition sg_bbl.h:981
bbl_mod::madsen94
subroutine madsen94(ubr, wr, ucr, zr, phiwc, kn, dstp, ustrc, ustrwm, ustrr, fwc, zoa, dwc_va)
Definition ssw_bbl.h:1534
bc_2d_mod
Definition bc_2d.F:2
bc_2d_mod::bc_r2d_tile
subroutine bc_r2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
Definition bc_2d.F:44
bc_2d_mod::bc_v2d_tile
subroutine bc_v2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
Definition bc_2d.F:345
bc_2d_mod::bc_u2d_tile
subroutine bc_u2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
Definition bc_2d.F:167
mod_bbl
Definition mod_bbl.F:2
mod_bbl::bbl
type(t_bbl), dimension(:), allocatable bbl
Definition mod_bbl.F:62
mod_forces
Definition mod_forces.F:2
mod_forces::forces
type(t_forces), dimension(:), allocatable forces
Definition mod_forces.F:554
mod_grid
Definition mod_grid.F:2
mod_grid::grid
type(t_grid), dimension(:), allocatable grid
Definition mod_grid.F:365
mod_iounits
Definition mod_iounits.F:2
mod_iounits::stdout
integer stdout
Definition mod_iounits.F:195
mod_ocean
Definition mod_ocean.F:2
mod_ocean::ocean
type(t_ocean), dimension(:), allocatable ocean
Definition mod_ocean.F:351
mod_parallel
Definition mod_parallel.F:2
mod_parallel::master
logical master
Definition mod_parallel.F:40
mod_param
Definition mod_param.F:2
mod_param::inlm
integer, parameter inlm
Definition mod_param.F:662
mod_param::nghostpoints
integer nghostpoints
Definition mod_param.F:710
mod_scalars
Definition mod_scalars.F:2
mod_scalars::sg_kappa
real(dp) sg_kappa
Definition mod_scalars.F:1541
mod_scalars::sg_z100
real(dp) sg_z100
Definition mod_scalars.F:1545
mod_scalars::sg_mp
complex(c8) sg_mp
Definition mod_scalars.F:1551
mod_scalars::cdb_min
real(dp) cdb_min
Definition mod_scalars.F:772
mod_scalars::vonkar
real(dp) vonkar
Definition mod_scalars.F:469
mod_scalars::dt
real(dp), dimension(:), allocatable dt
Definition mod_scalars.F:269
mod_scalars::sg_n
integer, parameter sg_n
Definition mod_scalars.F:1534
mod_scalars::ewperiodic
logical, dimension(:), allocatable ewperiodic
Definition mod_scalars.F:510
mod_scalars::nsperiodic
logical, dimension(:), allocatable nsperiodic
Definition mod_scalars.F:511
mod_scalars::sg_tol
real(dp) sg_tol
Definition mod_scalars.F:1543
mod_scalars::sg_alpha
real(dp) sg_alpha
Definition mod_scalars.F:1539
mod_scalars::g
real(dp) g
Definition mod_scalars.F:464
mod_scalars::cdb_max
real(dp) cdb_max
Definition mod_scalars.F:773
mod_scalars::sg_z1p
real(dp) sg_z1p
Definition mod_scalars.F:1547
mod_scalars::sg_z1min
real(dp) sg_z1min
Definition mod_scalars.F:1546
mod_scalars::pi
real(dp), parameter pi
Definition mod_scalars.F:813
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::iwsed
integer, parameter iwsed
Definition sediment_mod.h:126
mod_sediment::irlen
integer, parameter irlen
Definition sediment_mod.h:128
mod_sediment::irhgt
integer, parameter irhgt
Definition sediment_mod.h:129
mod_sediment::iznik
integer, parameter iznik
Definition sediment_mod.h:133
mod_sediment::ibwav
integer, parameter ibwav
Definition sediment_mod.h:130
mod_sediment::idens
integer, parameter idens
Definition sediment_mod.h:125
mod_sediment::isd50
integer, parameter isd50
Definition sediment_mod.h:124
mod_sediment::izbld
integer, parameter izbld
Definition sediment_mod.h:136
mod_sediment::izapp
integer, parameter izapp
Definition sediment_mod.h:132
mod_sediment::izbio
integer, parameter izbio
Definition sediment_mod.h:134
mod_sediment::izdef
integer, parameter izdef
Definition sediment_mod.h:131
mod_sediment::itauc
integer, parameter itauc
Definition sediment_mod.h:127
mod_sediment::izwbl
integer, parameter izwbl
Definition sediment_mod.h:137
mod_sediment::izbfm
integer, parameter izbfm
Definition sediment_mod.h:135
mod_stepping
Definition mod_stepping.F:2
mod_stepping::nrhs
integer, dimension(:), allocatable nrhs
Definition mod_stepping.F:70
mp_exchange_mod
Definition mp_exchange.F:2
mp_exchange_mod::mp_exchange2d
subroutine mp_exchange2d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, nghost, ew_periodic, ns_periodic, a, b, c, d)
Definition mp_exchange.F:294
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