ROMS
Loading...
Searching...
No Matches
sg_bbl.h
Go to the documentation of this file.
1 MODULE bbl_mod
2!
3!git $Id$
4!================================================== Hernan G. Arango ===
5! Copyright (c) 2002-2025 The ROMS Group Richard Styles !
6! Licensed under a MIT/X style license !
7! See License_ROMS.md !
8!=======================================================================
9! !
10! This routine compute bottom stresses for the case when the wave !
11! solution in the wave boundary layer is based on a 2-layer eddy !
12! viscosity that is linear increasing above Zo and constant above !
13! Z1. !
14! !
15! Reference: !
16! !
17! Styles, R. and S.M. glenn, 2000: Modeling stratified wave and !
18! current bottom boundary layers in the continental shelf, JGR, !
19! 105, 24119-24139. !
20! !
21!=======================================================================
22!
23 implicit none
24!
25 PRIVATE
26 PUBLIC :: bblm
27!
28 CONTAINS
29!
30!***********************************************************************
31 SUBROUTINE bblm (ng, tile)
32!***********************************************************************
33!
34 USE mod_param
35 USE mod_bbl
36 USE mod_forces
37 USE mod_grid
38 USE mod_ocean
39 USE mod_sedbed
40 USE mod_stepping
41!
42! Imported variable declarations.
43!
44 integer, intent(in) :: ng, tile
45!
46! Local variable declarations.
47!
48 character (len=*), parameter :: MyFile = &
49 & __FILE__
50!
51#include "tile.h"
52!
53#ifdef PROFILE
54 CALL wclock_on (ng, inlm, 37, __line__, myfile)
55#endif
56 CALL sg_bbl_tile (ng, tile, &
57 & lbi, ubi, lbj, ubj, &
58 & imins, imaxs, jmins, jmaxs, &
59 & nrhs(ng), &
60 & grid(ng) % h, &
61 & grid(ng) % z_r, &
62 & grid(ng) % z_w, &
63 & grid(ng) % angler, &
64 & grid(ng) % ZoBot, &
65#if defined SG_CALC_UB
66 & forces(ng) % Hwave, &
67#else
68 & forces(ng) % Uwave_rms, &
69#endif
70 & forces(ng) % Dwave, &
71 & forces(ng) % Pwave_bot, &
72 & ocean(ng) % rho, &
73 & ocean(ng) % u, &
74 & ocean(ng) % v, &
75 & sedbed(ng) % bottom, &
76 & bbl(ng) % Iconv, &
77 & bbl(ng) % Ubot, &
78 & bbl(ng) % Vbot, &
79 & bbl(ng) % Ur, &
80 & bbl(ng) % Vr, &
81 & bbl(ng) % bustrc, &
82 & bbl(ng) % bvstrc, &
83 & bbl(ng) % bustrw, &
84 & bbl(ng) % bvstrw, &
85 & bbl(ng) % bustrcwmax, &
86 & bbl(ng) % bvstrcwmax, &
87 & forces(ng) % bustr, &
88 & forces(ng) % bvstr)
89#ifdef PROFILE
90 CALL wclock_off (ng, inlm, 37, __line__, myfile)
91#endif
92!
93 RETURN
94 END SUBROUTINE bblm
95!
96!***********************************************************************
97 SUBROUTINE sg_bbl_tile (ng, tile, &
98 & LBi, UBi, LBj, UBj, &
99 & IminS, ImaxS, JminS, JmaxS, &
100 & nrhs, &
101 & h, z_r, z_w, angler, ZoBot, &
102#if defined SG_CALC_UB
103 & Hwave, &
104#else
105 & Uwave_rms, &
106#endif
107 & Dwave, Pwave_bot, &
108 & rho, u, v, &
109 & bottom, &
110 & Iconv, &
111 & Ubot, Vbot, Ur, Vr, &
112 & bustrc, bvstrc, &
113 & bustrw, bvstrw, &
114 & bustrcwmax, bvstrcwmax, &
115 & bustr, bvstr)
116!***********************************************************************
117!
118 USE mod_param
119 USE mod_scalars
120 USE mod_sediment
121!
122 USE bc_2d_mod
123#ifdef DISTRIBUTE
124 USE mp_exchange_mod, ONLY : mp_exchange2d
125#endif
126!
127! Imported variable declarations.
128!
129 integer, intent(in) :: ng, tile
130 integer, intent(in) :: LBi, UBi, LBj, UBj
131 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
132 integer, intent(in) :: nrhs
133!
134#ifdef ASSUMED_SHAPE
135 integer, intent(inout) :: Iconv(LBi:,LBj:)
136
137 real(r8), intent(in) :: h(LBi:,LBj:)
138 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
139 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
140 real(r8), intent(in) :: angler(LBi:,LBj:)
141 real(r8), intent(in) :: ZoBot(LBi:,LBj:)
142# if defined SG_CALC_UB
143 real(r8), intent(in) :: Hwave(LBi:,LBj:)
144# else
145 real(r8), intent(in) :: Uwave_rms(LBi:,LBj:)
146# endif
147 real(r8), intent(in) :: Dwave(LBi:,LBj:)
148 real(r8), intent(in) :: Pwave_bot(LBi:,LBj:)
149 real(r8), intent(in) :: rho(LBi:,LBj:,:)
150 real(r8), intent(in) :: u(LBi:,LBj:,:,:)
151 real(r8), intent(in) :: v(LBi:,LBj:,:,:)
152
153 real(r8), intent(inout) :: bottom(LBi:,LBj:,:)
154
155 real(r8), intent(out) :: Ubot(LBi:,LBj:)
156 real(r8), intent(out) :: Vbot(LBi:,LBj:)
157 real(r8), intent(out) :: Ur(LBi:,LBj:)
158 real(r8), intent(out) :: Vr(LBi:,LBj:)
159 real(r8), intent(out) :: bustrc(LBi:,LBj:)
160 real(r8), intent(out) :: bvstrc(LBi:,LBj:)
161 real(r8), intent(out) :: bustrw(LBi:,LBj:)
162 real(r8), intent(out) :: bvstrw(LBi:,LBj:)
163 real(r8), intent(out) :: bustrcwmax(LBi:,LBj:)
164 real(r8), intent(out) :: bvstrcwmax(LBi:,LBj:)
165 real(r8), intent(out) :: bustr(LBi:,LBj:)
166 real(r8), intent(out) :: bvstr(LBi:,LBj:)
167#else
168 integer, intent(inout) :: Iconv(LBi:UBi,LBj:UBj)
169
170 real(r8), intent(in) :: h(LBi:UBi,LBj:UBj)
171 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
172 real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
173 real(r8), intent(in) :: angler(LBi:UBi,LBj:UBj)
174 real(r8), intent(in) :: ZoBot(LBi:UBi,LBj:UBj)
175# if defined SG_CALC_UB
176 real(r8), intent(in) :: Hwave(LBi:UBi,LBj:UBj)
177# else
178 real(r8), intent(in) :: Uwave_rms(LBi:UBi,LBj:UBj)
179# endif
180 real(r8), intent(in) :: Dwave(LBi:UBi,LBj:UBj)
181 real(r8), intent(in) :: Pwave_bot(LBi:UBi,LBj:UBj)
182 real(r8), intent(in) :: rho(LBi:UBi,LBj:UBj,N(ng))
183 real(r8), intent(in) :: u(LBi:UBi,LBj:UBj,N(ng),2)
184 real(r8), intent(in) :: v(LBi:UBi,LBj:UBj,N(ng),2)
185
186 real(r8), intent(inout) :: bottom(LBi:UBi,LBj:UBj,MBOTP)
187
188 real(r8), intent(out) :: Ubot(LBi:UBi,LBj:UBj)
189 real(r8), intent(out) :: Vbot(LBi:UBi,LBj:UBj)
190 real(r8), intent(out) :: Ur(LBi:UBi,LBj:UBj)
191 real(r8), intent(out) :: Vr(LBi:UBi,LBj:UBj)
192 real(r8), intent(out) :: bustrc(LBi:UBi,LBj:UBj)
193 real(r8), intent(out) :: bvstrc(LBi:UBi,LBj:UBj)
194 real(r8), intent(out) :: bustrw(LBi:UBi,LBj:UBj)
195 real(r8), intent(out) :: bvstrw(LBi:UBi,LBj:UBj)
196 real(r8), intent(out) :: bustrcwmax(LBi:UBi,LBj:UBj)
197 real(r8), intent(out) :: bvstrcwmax(LBi:UBi,LBj:UBj)
198 real(r8), intent(out) :: bustr(LBi:UBi,LBj:UBj)
199 real(r8), intent(out) :: bvstr(LBi:UBi,LBj:UBj)
200#endif
201!
202! Local variable declarations.
203!
204 logical :: ITERATE
205
206 integer :: Iter, i, j, k
207
208 real(r8), parameter :: eps = 1.0e-10_r8
209
210 real(r8) :: Fwave_bot, Kb, Kbh, KboKb0, Kb0, Kdelta, Ustr
211 real(r8) :: anglec, anglew
212 real(r8) :: cff, cff1, cff2, cff3, og, fac, fac1, fac2
213 real(r8) :: sg_ab, sg_abokb, sg_a1, sg_b1, sg_chi, sg_c1, sg_dd
214 real(r8) :: sg_epsilon, sg_eta, sg_fofa, sg_fofb, sg_fofc, sg_fwm
215 real(r8) :: sg_kbs, sg_lambda, sg_mu, sg_phicw, sg_ro, sg_row
216 real(r8) :: sg_shdnrm, sg_shld, sg_shldcr, sg_scf, sg_ss, sg_star
217 real(r8) :: sg_ub, sg_ubokur, sg_ubouc, sg_ubouwm, sg_ur
218 real(r8) :: sg_ustarc, sg_ustarcw, sg_ustarwm, sg_znot, sg_znotp
219 real(r8) :: sg_zr, sg_zrozn, sg_z1, sg_z1ozn, sg_z2, twopi, z1, z2
220
221 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Ab
222 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Tauc
223 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Tauw
224 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Taucwmax
225 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Ur_sg
226 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Vr_sg
227 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Ub
228 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Ucur
229 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Umag
230 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Vcur
231 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Zr
232 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: phic
233 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: phicw
234 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: rheight
235 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: rlength
236 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: u100
237 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: znot
238 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: znotc
239
240#include "set_bounds.h"
241!
242!-----------------------------------------------------------------------
243! Initalize to default values.
244!-----------------------------------------------------------------------
245!
246 DO j=jstrv-1,jend
247 DO i=istru-1,iend
248 tauc(i,j)=0.0_r8
249 tauw(i,j)=0.0_r8
250 taucwmax(i,j)=0.0_r8
251 u100(i,j)=0.0_r8
252 rheight(i,j)=0.0_r8
253 rlength(i,j)=0.0_r8
254 znot(i,j)=zobot(i,j)
255 znotc(i,j)=0.0_r8
256 END DO
257 END DO
258!
259!-----------------------------------------------------------------------
260! Set currents above bed.
261!-----------------------------------------------------------------------
262!
263 DO j=jstrv-1,jend+1
264 DO i=istru-1,iend+1
265 zr(i,j)=z_r(i,j,1)-z_w(i,j,0)
266 ur_sg(i,j)=u(i,j,1,nrhs)
267 vr_sg(i,j)=v(i,j,1,nrhs)
268#ifdef SG_LOGINT
269!
270! If current height is less than z1ur, interpolate logarithmically
271! to z1ur.
272!
273 IF (zr(i,j).lt.sg_z1min) THEN
274 DO k=2,n(ng)
275 z1=z_r(i,j,k-1)-z_w(i,j,0)
276 z2=z_r(i,j,k )-z_w(i,j,0)
277 IF ((z1.lt.sg_z1min).and.(sg_z1min.lt.z2)) THEN
278 fac=1.0_r8/log(z2/z1)
279 fac1=fac*log(z2/sg_z1min)
280 fac2=fac*log(sg_z1min/z1)
281 ur_sg(i,j)=fac1*u(i,j,k-1,nrhs)+fac2*u(i,j,k,nrhs)
282 vr_sg(i,j)=fac1*v(i,j,k-1,nrhs)+fac2*v(i,j,k,nrhs)
283 zr(i,j)=sg_z1min
284 END IF
285 END DO
286 END IF
287#endif
288 END DO
289 END DO
290!
291!-----------------------------------------------------------------------
292! Compute bed wave orbital velocity (m/s) and excursion amplitude
293! (m) from wind-induced waves. Use linear wave theory dispersion
294! relation for wave number.
295!-----------------------------------------------------------------------
296!
297 twopi=2.0_r8*pi
298 og=1.0_r8/g
299 DO j=jstrv-1,jend
300 DO i=istru-1,iend
301!
302! Compute first guess for wavenumber, Kb0. Use deep water (Kb0*h>1)
303! and shallow water (Kb0*H<1) approximations.
304!
305 fwave_bot=twopi/max(pwave_bot(i,j),0.05_r8)
306!
307! Compute bed wave orbital velocity and excursion amplitude.
308!
309#ifdef SG_CALC_UB
310 kb0=fwave_bot*fwave_bot*og
311 IF (kb0*h(i,j).ge.1.0_r8) THEN
312 kb=kb0
313 ELSE
314 kb=fwave_bot/sqrt(g*h(i,j))
315 END IF
316!
317! Compute bottom wave number via Newton-Raphson method.
318!
319 iterate=.true.
320 DO iter=1,sg_n
321 IF (iterate) THEN
322 kbh=kb*h(i,j)
323 kbokb0=kb/kb0
324 kdelta=(1.0_r8-kbokb0*tanh(kbh))/ &
325 & (1.0_r8+kbh*(kbokb0-1.0_r8/kbokb0))
326 iterate=abs(kb*kdelta) .ge. sg_tol
327 kb=kb*(1.0_r8+kdelta)
328 END IF
329 END DO
330 ab(i,j)=0.5_r8*hwave(i,j)/sinh(kb*h(i,j))+eps
331 ub(i,j)=fwave_bot*ab(i,j)+eps
332#else
333 ub(i,j)=abs(uwave_rms(i,j))+eps
334 ab(i,j)=ub(i,j)/fwave_bot+eps
335#endif
336!
337! Compute bottom current magnitude at RHO-points.
338!
339 ucur(i,j)=0.5_r8*(ur_sg(i,j)+ur_sg(i+1,j))
340 vcur(i,j)=0.5_r8*(vr_sg(i,j)+vr_sg(i,j+1))
341 umag(i,j)=sqrt(ucur(i,j)*ucur(i,j)+vcur(i,j)*vcur(i,j))+eps
342!
343! Compute angle between currents and waves (radians)
344!
345 IF (ucur(i,j).eq.0.0_r8) THEN
346 phic(i,j)=0.5_r8*pi*sign(1.0_r8,vcur(i,j))
347 ELSE
348 phic(i,j)=atan2(vcur(i,j),ucur(i,j))
349 ENDIF
350 phicw(i,j)=1.5_r8*pi-dwave(i,j)-phic(i,j)-angler(i,j)
351 END DO
352 END DO
353!
354!-----------------------------------------------------------------------
355! Set default logarithmic profile.
356!-----------------------------------------------------------------------
357!
358 DO j=jstrv-1,jend
359 DO i=istru-1,iend
360 IF (umag(i,j).gt.0.0_r8) THEN
361!! Ustr=MIN(sg_ustarcdef,Umag(i,j)*vonKar/ &
362!! & LOG(Zr(i,j)/ZoBot(i,j)))
363!! Tauc(i,j)=Ustr*Ustr
364 cff1=vonkar/log(zr(i,j)/zobot(i,j))
365 cff2=min(cdb_max,max(cdb_min,cff1*cff1))
366 tauc(i,j)=cff2*umag(i,j)*umag(i,j)
367 END IF
368 END DO
369 END DO
370!
371!-----------------------------------------------------------------------
372! Wave-current interaction case.
373!-----------------------------------------------------------------------
374!
375 DO j=jstrv-1,jend
376 DO i=istru-1,iend
377 sg_dd=bottom(i,j,isd50)
378 sg_ss=bottom(i,j,idens)/(rho(i,j,1)+1000.0_r8)
379 sg_ab=ab(i,j)
380 sg_ub=ub(i,j)
381 sg_phicw=phicw(i,j)
382 sg_ur=umag(i,j)
383 sg_zr=zr(i,j)
384!
385! Compute hydraulic roughness "Znot" (m), ripple height "eta" (m),
386! and ripple length "lambda" (m).
387!
388#ifdef SG_CALC_ZNOT
389 sg_star=sg_dd/(4.0_r8*sg_nu)*sqrt((sg_ss-1.0_r8)*sg_g*sg_dd)
390!
391! Compute critical shield parameter based on grain diameter.
392! (sg_scf is a correction factor).
393!
394 sg_scf=1.0_r8
395 IF (sg_star.le.1.5_r8) THEN
396 sg_shldcr=sg_scf*0.0932_r8*sg_star**(-0.707_r8)
397 ELSE IF ((1.5_r8.lt.sg_star).and.(sg_star.lt.4.0_r8)) THEN
398 sg_shldcr=sg_scf*0.0848_r8*sg_star**(-0.473_r8)
399 ELSE IF ((4.0_r8.le.sg_star).and.(sg_star.lt.10.0_r8)) THEN
400 sg_shldcr=sg_scf*0.0680_r8*sg_star**(-0.314_r8)
401 ELSE IF ((10.0_r8.le.sg_star).and.(sg_star.lt.34.0_r8)) THEN
402 sg_shldcr=sg_scf*0.033_r8
403 ELSE IF ((34.0_r8.le.sg_star).and.(sg_star.lt.270.0_r8)) THEN
404 sg_shldcr=sg_scf*0.0134_r8*sg_star**(0.255_r8)
405 ELSE
406 sg_shldcr=sg_scf*0.056_r8
407 END IF
408!
409! Calculate skin friction shear stress based on Ole Madsen (1994)
410! empirical formula. Check initiation of sediment motion criteria,
411! to see if we compute sg_znot based on the wave-formed ripples.
412! If the skin friction calculation indicates that sediment is NOT
413! in motion, the ripple model is invalid and take the default value,
414! ZoBot.
415!
416 sg_abokb=sg_ab/sg_dd
417 IF (sg_abokb.le.100.0_r8) THEN
418 sg_fwm=exp(7.02_r8*sg_abokb**(-0.078_r8)-8.82_r8)
419 ELSE
420 sg_fwm=exp(5.61_r8*sg_abokb**(-0.109_r8)-7.30_r8)
421 END IF
422 sg_ustarwm=sqrt(0.5_r8*sg_fwm)*sg_ub
423 sg_shdnrm=(sg_ss-1.0_r8)*sg_dd*sg_g
424 sg_shld=sg_ustarwm*sg_ustarwm/sg_shdnrm
425 IF ((sg_shld/sg_shldcr).le.1.0_r8) THEN
426 sg_znot=zobot(i,j)
427 sg_eta=0.0_r8
428 sg_lambda=0.0_r8
429 ELSE
430!
431! Calculate ripple height and length and bottom roughness
432!
433 sg_chi=4.0_r8*sg_nu*sg_ub*sg_ub/ &
434 & (sg_dd*((sg_ss-1.0_r8)*sg_g*sg_dd)**1.5_r8)
435 IF (sg_chi.le.2.0_r8) THEN
436 sg_eta=sg_ab*0.30_r8*sg_chi**(-0.39_r8)
437 sg_lambda=sg_ab*1.96_r8*sg_chi**(-0.28_r8)
438 ELSE
439 sg_eta=sg_ab*0.45_r8*sg_chi**(-0.99_r8)
440 sg_lambda=sg_ab*2.71_r8*sg_chi**(-0.75_r8)
441 END IF
442 sg_kbs=sg_ab*0.0655_r8* &
443 & (sg_ub*sg_ub/((sg_ss-1.0_r8)*sg_g*sg_ab))**1.4_r8
444 sg_znot=(sg_dd+2.3_r8*sg_eta+sg_kbs)/30.0_r8
445 END IF
446#else
447 sg_znot=zobot(i,j)
448 sg_chi=4.0_r8*sg_nu*sg_ub*sg_ub/ &
449 & (sg_dd*((sg_ss-1.0_r8)*sg_g*sg_dd)**1.5_r8)
450 IF (sg_chi.le.2.0_r8) THEN
451 sg_eta=sg_ab*0.32_r8*sg_chi**(-0.34_r8)
452 sg_lambda=sg_ab*2.04_r8*sg_chi**(-0.23_r8)
453 ELSE
454 sg_eta=sg_ab*0.52_r8*sg_chi**(-1.01_r8)
455 sg_lambda=sg_ab*2.7_r8*sg_chi**(-0.78_r8)
456 END IF
457#endif
458 znot(i,j)=sg_znot
459 rheight(i,j)=sg_eta
460 rlength(i,j)=sg_lambda
461!
462! Compute only when nonzero currents and waves.
463!
464 sg_zrozn=sg_zr/sg_znot
465 IF ((sg_ur.gt.0.0_r8).and.(sg_ub.gt.0.0_r8).and. &
466 & (sg_zrozn.gt.1.0_r8)) THEN
467!
468! Compute bottom stress based on ripple roughness.
469!
470 sg_ubokur=sg_ub/(sg_kappa*sg_ur)
471 sg_row=sg_ab/sg_znot
472 sg_a1=1.0e-6_r8
473 CALL sg_bstress (sg_row, sg_zrozn, sg_phicw, sg_ubokur, &
474 & sg_a1, sg_mu, sg_epsilon, sg_ro, sg_fofa)
475 sg_abokb=sg_ab/(30.0_r8*sg_znot)
476 IF (sg_abokb.le.100.0_r8) THEN
477 sg_fwm=exp(-8.82_r8+7.02_r8*sg_abokb**(-0.078_r8))
478 ELSE
479 sg_fwm=exp(-7.30_r8+5.61_r8*sg_abokb**(-0.109_r8))
480 END IF
481 sg_ubouwm=sqrt(2.0_r8/sg_fwm)
482!
483! Determine the maximum ratio of wave over combined shear stresses,
484! sg_ubouwm (ub/ustarwm).
485!
486 CALL sg_purewave (sg_row, sg_ubouwm, sg_znotp, sg_ro)
487!
488! Set initial guess of the ratio of wave over shear stress, sg_c1
489! (ub/ustarc).
490!
491 sg_b1=sg_ubouwm
492 sg_fofb=-sg_fofa
493 sg_c1=0.5_r8*(sg_a1+sg_b1)
494 CALL sg_bstress (sg_row, sg_zrozn, sg_phicw, sg_ubokur, &
495 & sg_c1, sg_mu, sg_epsilon, sg_ro, sg_fofc)
496!
497! Solve PDE via bi-section method.
498!
499 iterate=.true.
500 DO iter=1,sg_n
501 IF (iterate) THEN
502 IF ((sg_fofb*sg_fofc).lt.0.0_r8) THEN
503 sg_a1=sg_c1
504 ELSE
505 sg_b1=sg_c1
506 END IF
507 sg_c1=0.5_r8*(sg_a1+sg_b1)
508 CALL sg_bstress (sg_row, sg_zrozn, sg_phicw, sg_ubokur, &
509 & sg_c1, sg_mu, sg_epsilon, sg_ro, &
510 & sg_fofc)
511 iterate=(sg_b1-sg_c1) .ge. sg_tol
512 IF (iterate) iconv(i,j)=iter
513 END IF
514 END DO
515 sg_ubouc=sg_c1
516!
517! Compute bottom shear stress magnitude (m/s).
518!
519 sg_ustarcw=sg_ub/sg_ubouc
520 sg_ustarwm=sg_mu*sg_ustarcw
521!! sg_ustarc=MIN(sg_ustarcdef,sg_epsilon*sg_ustarcw)
522!! sg_ustarc=MIN(SQRT(Tauc(i,j)),sg_epsilon*sg_ustarcw)
523 sg_ustarc=max(sqrt(tauc(i,j)),sg_epsilon*sg_ustarcw)
524 tauc(i,j)=sg_ustarc*sg_ustarc
525 tauw(i,j)=sg_ustarwm*sg_ustarwm
526 taucwmax(i,j)=sqrt((tauc(i,j)+ &
527 & tauw(i,j)*cos(phicw(i,j)))**2+ &
528 & (tauw(i,j)*sin(phicw(i,j)))**2)
529!
530! Compute apparent hydraulic roughness (m).
531!
532 IF (sg_epsilon.gt.0.0_r8) THEN
533 sg_z1=sg_alpha*sg_kappa*sg_ab/sg_ubouc
534 sg_z2=sg_z1/sg_epsilon
535 sg_z1ozn=sg_z1/sg_znot
536 znotc(i,j)=sg_z2* &
537 & exp(-(1.0_r8-sg_epsilon+ &
538 & sg_epsilon*log(sg_z1ozn)))
539!
540! Compute mean (m/s) current at 100 cm above the bottom.
541!
542 IF (sg_z100.gt.sg_z2) THEN
543 u100(i,j)=sg_ustarc* &
544 & (log(sg_z100/sg_z2)+1.0_r8-sg_epsilon+ &
545 & sg_epsilon*log(sg_z1ozn))/sg_kappa
546 ELSE IF ((sg_z100.le.sg_z2).and.(sg_zr.gt.sg_z1)) THEN
547 u100(i,j)=sg_ustarc*sg_epsilon* &
548 & (sg_z100/sg_z1-1.0_r8+log(sg_z1ozn))/sg_kappa
549 ELSE
550 u100(i,j)=sg_ustarc*sg_epsilon* &
551 & log(sg_z100/sg_znot)/sg_kappa
552 END IF
553 END IF
554 END IF
555 END DO
556 END DO
557!
558!-----------------------------------------------------------------------
559! Compute kinematic bottom stress components due current and wind-
560! induced waves.
561!-----------------------------------------------------------------------
562!
563 DO j=jstr,jend
564 DO i=istru,iend
565 anglec=ur_sg(i,j)/(0.5*(umag(i-1,j)+umag(i,j)))
566 bustr(i,j)=0.5_r8*(tauc(i-1,j)+tauc(i,j))*anglec
567# ifdef WET_DRY
568 cff2=0.75_r8*0.5_r8*(z_w(i-1,j,1)+z_w(i,j,1)- &
569 & z_w(i-1,j,0)-z_w(i,j,0))
570 bustr(i,j)=sign(1.0_r8,bustr(i,j))*min(abs(bustr(i,j)), &
571 & abs(u(i,j,1,nrhs))*cff2/dt(ng))
572# endif
573 END DO
574 END DO
575 DO j=jstrv,jend
576 DO i=istr,iend
577 anglec=vr_sg(i,j)/(0.5_r8*(umag(i,j-1)+umag(i,j)))
578 bvstr(i,j)=0.5_r8*(tauc(i,j-1)+tauc(i,j))*anglec
579# ifdef WET_DRY
580 cff2=0.75_r8*0.5_r8*(z_w(i,j-1,1)+z_w(i,j,1)- &
581 & z_w(i,j-1,0)-z_w(i,j,0))
582 bvstr(i,j)=sign(1.0_r8,bvstr(i,j))*min(abs(bvstr(i,j)), &
583 & abs(v(i,j,1,nrhs))*cff2/dt(ng))
584# endif
585 END DO
586 END DO
587 DO j=jstr,jend
588 DO i=istr,iend
589 anglec=ucur(i,j)/umag(i,j)
590 anglew=cos(1.5_r8*pi-dwave(i,j)-angler(i,j))
591 bustrc(i,j)=tauc(i,j)*anglec
592 bustrw(i,j)=tauw(i,j)*anglew
593 bustrcwmax(i,j)=taucwmax(i,j)*anglew
594 ubot(i,j)=ub(i,j)*anglew
595 ur(i,j)=ucur(i,j)
596!
597 anglec=vcur(i,j)/umag(i,j)
598 anglew=sin(1.5_r8*pi-dwave(i,j)-angler(i,j))
599 bvstrc(i,j)=tauc(i,j)*anglec
600 bvstrw(i,j)=tauw(i,j)*anglew
601 bvstrcwmax(i,j)=taucwmax(i,j)*anglew
602 vbot(i,j)=ub(i,j)*anglew
603 vr(i,j)=vcur(i,j)
604!
605 bottom(i,j,ibwav)=ab(i,j)
606 bottom(i,j,irhgt)=rheight(i,j)
607 bottom(i,j,irlen)=rlength(i,j)
608 bottom(i,j,izdef)=znot(i,j)
609 bottom(i,j,izapp)=znotc(i,j)
610 END DO
611 END DO
612!
613! Apply periodic or gradient boundary conditions for output
614! purposes only.
615!
616 CALL bc_u2d_tile (ng, tile, &
617 & lbi, ubi, lbj, ubj, &
618 & bustr)
619 CALL bc_v2d_tile (ng, tile, &
620 & lbi, ubi, lbj, ubj, &
621 & bvstr)
622 CALL bc_r2d_tile (ng, tile, &
623 & lbi, ubi, lbj, ubj, &
624 & bustrc)
625 CALL bc_r2d_tile (ng, tile, &
626 & lbi, ubi, lbj, ubj, &
627 & bvstrc)
628 CALL bc_r2d_tile (ng, tile, &
629 & lbi, ubi, lbj, ubj, &
630 & bustrw)
631 CALL bc_r2d_tile (ng, tile, &
632 & lbi, ubi, lbj, ubj, &
633 & bvstrw)
634 CALL bc_u2d_tile (ng, tile, &
635 & lbi, ubi, lbj, ubj, &
636 & bustrcwmax)
637 CALL bc_r2d_tile (ng, tile, &
638 & lbi, ubi, lbj, ubj, &
639 & bvstrcwmax)
640 CALL bc_r2d_tile (ng, tile, &
641 & lbi, ubi, lbj, ubj, &
642 & ubot)
643 CALL bc_r2d_tile (ng, tile, &
644 & lbi, ubi, lbj, ubj, &
645 & vbot)
646 CALL bc_r2d_tile (ng, tile, &
647 & lbi, ubi, lbj, ubj, &
648 & ur)
649 CALL bc_r2d_tile (ng, tile, &
650 & lbi, ubi, lbj, ubj, &
651 & vr)
652 CALL bc_r2d_tile (ng, tile, &
653 & lbi, ubi, lbj, ubj, &
654 & bottom(:,:,ibwav))
655 CALL bc_r2d_tile (ng, tile, &
656 & lbi, ubi, lbj, ubj, &
657 & bottom(:,:,irhgt))
658 CALL bc_r2d_tile (ng, tile, &
659 & lbi, ubi, lbj, ubj, &
660 & bottom(:,:,irlen))
661 CALL bc_r2d_tile (ng, tile, &
662 & lbi, ubi, lbj, ubj, &
663 & bottom(:,:,izdef))
664 CALL bc_r2d_tile (ng, tile, &
665 & lbi, ubi, lbj, ubj, &
666 & bottom(:,:,izapp))
667#ifdef DISTRIBUTE
668 CALL mp_exchange2d (ng, tile, inlm, 4, &
669 & lbi, ubi, lbj, ubj, &
670 & nghostpoints, &
671 & ewperiodic(ng), nsperiodic(ng), &
672 & bustr, bvstr, bustrc, bvstrc)
673 CALL mp_exchange2d (ng, tile, inlm, 4, &
674 & lbi, ubi, lbj, ubj, &
675 & nghostpoints, &
676 & ewperiodic(ng), nsperiodic(ng), &
677 & bustrw, bvstrw, bustrcwmax, bvstrcwmax)
678 CALL mp_exchange2d (ng, tile, inlm, 4, &
679 & lbi, ubi, lbj, ubj, &
680 & nghostpoints, &
681 & ewperiodic(ng), nsperiodic(ng), &
682 & ubot, vbot, ur, vr)
683 CALL mp_exchange2d (ng, tile, inlm, 3, &
684 & lbi, ubi, lbj, ubj, &
685 & nghostpoints, &
686 & ewperiodic(ng), nsperiodic(ng), &
687 & bottom(:,:,ibwav), &
688 & bottom(:,:,irhgt), &
689 & bottom(:,:,irlen))
690 CALL mp_exchange2d (ng, tile, inlm, 4, &
691 & lbi, ubi, lbj, ubj, &
692 & nghostpoints, &
693 & ewperiodic(ng), nsperiodic(ng), &
694 & bottom(:,:,izdef), &
695 & bottom(:,:,izapp))
696#endif
697!
698 RETURN
699 END SUBROUTINE sg_bbl_tile
700!
701 SUBROUTINE sg_bstress (sg_row, sg_zrozn, sg_phicw, sg_ubokur, &
702 & sg_ubouc, sg_mu, sg_epsilon, sg_ro, &
703 & sg_fofx)
704!
705!=======================================================================
706! !
707! This routine computes bottom stresses via bottom boundary layer !
708! formulation of Styles and Glenn (1999). !
709! !
710! On Input: !
711! !
712! sg_row Ratio of wave excursion amplitude over roughness. !
713! sg_zrozn Ratio of height of current over roughness. !
714! sg_phiwc Angle between wave and currents (radians). !
715! sg_ubokur Ratio of wave over current velocity: !
716! ub/(vonKar*ur) !
717! sg_ubouc Ratio of bed wave orbital over bottom shear stress !
718! (ub/ustarc), first guess. !
719! !
720! On Output: !
721! !
722! sg_ubouc Ratio of bed wave orbital over bottom shear stress !
723! (ub/ustarc), iterated value. !
724! sg_mu Ratio between wave and current bottom shear !
725! stresses (ustarwm/ustarc). !
726! sg_epsilon Ratio between combined (wave and current) and !
727! current bottom shear stresses (ustarc/ustarcw). !
728! sg_ro Internal friction Rossby number: !
729! ustarc/(omega*znot) !
730! sg_fofx Root of PDE used for convergence. !
731! !
732!=======================================================================
733!
734 USE mod_param
735 USE mod_scalars
736!
737! Imported variable declarations.
738!
739 real(r8), intent(in) :: sg_row, sg_zrozn, sg_phicw, sg_ubokur
740
741 real(r8), intent(inout) :: sg_ubouc
742
743 real(r8), intent(out) :: sg_mu, sg_epsilon, sg_ro, sg_fofx
744!
745! Local variable declarations.
746!
747 logical :: ITERATE
748
749 integer :: Iter
750
751 real(r8) :: cff, sg_bei, sg_beip, sg_ber, sg_berp, sg_cosphi
752 real(r8) :: sg_eps2, sg_kei, sg_keip, sg_ker, sg_kerp, sg_mu2
753 real(r8) :: sg_phi, sg_ror, sg_x, sg_z2p, sg_znotp, sg_zroz1
754 real(r8) :: sg_zroz2, sg_z1ozn, sg_z2ozn
755
756 complex(c8) :: sg_argi, sg_bnot, sg_bnotp, sg_b1, sg_b1p
757 complex(c8) :: sg_gammai, sg_knot, sg_knotp, sg_k1, sg_k1p
758 complex(c8) :: sg_ll, sg_nn
759!
760!-----------------------------------------------------------------------
761! Compute bottom stresses.
762!-----------------------------------------------------------------------
763!
764! Compute nondimensional bottom wave shear, phi. Iterate to make
765! sure that there is an upper limit in "ubouc". It usually requires
766! only one pass.
767!
768 iterate=.true.
769 DO iter=1,sg_n
770 IF (iterate) THEN
771 sg_ro=sg_row/sg_ubouc
772 sg_znotp=1.0_r8/(sg_kappa*sg_ro)
773 IF ((sg_z1p/sg_znotp).gt.1.0_r8) THEN
774 sg_x=2.0_r8*sqrt(sg_znotp)
775 IF (sg_x.le.8.0_r8) THEN
776 CALL sg_kelvin8m (sg_x, sg_ber, sg_bei, sg_ker, sg_kei, &
777 & sg_berp, sg_beip, sg_kerp, sg_keip)
778 ELSE
779 CALL sg_kelvin8p (sg_x, sg_ker, sg_kei, sg_ber, sg_bei, &
780 & sg_kerp, sg_keip, sg_berp, sg_beip)
781 END IF
782 cff=1.0_r8/sqrt(sg_znotp)
783 sg_bnot =cmplx(sg_ber,sg_bei)
784 sg_knot =cmplx(sg_ker,sg_kei)
785 sg_bnotp=cmplx(sg_berp,sg_beip)*cff
786 sg_knotp=cmplx(sg_kerp,sg_keip)*cff
787!
788 sg_x=2.0_r8*sqrt(sg_z1p)
789 IF (sg_x.le.8.0_r8) THEN
790 CALL sg_kelvin8m (sg_x, sg_ber, sg_bei, sg_ker, sg_kei, &
791 & sg_berp, sg_beip, sg_kerp, sg_keip)
792 ELSE
793 CALL sg_kelvin8p (sg_x, sg_ker, sg_kei, sg_ber, sg_bei, &
794 & sg_kerp, sg_keip, sg_berp, sg_beip)
795 END IF
796 cff=1.0_r8/sqrt(sg_z1p)
797 sg_b1 =cmplx(sg_ber,sg_bei)
798 sg_k1 =cmplx(sg_ker,sg_kei)
799 sg_b1p=cmplx(sg_berp,sg_beip)*cff
800 sg_k1p=cmplx(sg_kerp,sg_keip)*cff
801!
802 sg_ll=sg_mp*sg_b1+sg_b1p
803 sg_nn=sg_mp*sg_k1+sg_k1p
804 sg_argi=sg_bnotp*sg_nn/(sg_bnot*sg_nn-sg_knot*sg_ll)+ &
805 & sg_knotp*sg_ll/(sg_knot*sg_ll-sg_bnot*sg_nn)
806 sg_gammai=-sg_kappa*sg_znotp*sg_argi
807 sg_phi=cabs(sg_gammai)
808 ELSE
809 sg_gammai=-sg_kappa*sg_z1p*sg_mp
810 sg_phi=cabs(sg_gammai)
811 END IF
812!
813 IF (sg_ubouc.gt.(1.0_r8/sg_phi)) THEN
814 sg_ubouc=1.0_r8/sg_phi
815 ELSE
816 iterate=.false.
817 END IF
818 END IF
819 END DO
820!
821! Compute ratio of wave over current bottom shear stresses.
822!
823 sg_mu=sqrt(sg_ubouc*sg_phi)
824!
825! Compute ratio of current over combined bottom shear stresses.
826!
827 IF (sg_mu.eq.1.0_r8) THEN
828 sg_epsilon=0.0_r8
829 ELSE
830 sg_mu2=sg_mu*sg_mu
831 sg_cosphi=abs(cos(sg_phicw))
832 sg_eps2=-sg_mu2*sg_cosphi+ &
833 & sqrt(1.0_r8+sg_mu2*sg_mu2*(sg_cosphi*sg_cosphi-1.0_r8))
834 sg_epsilon=sqrt(sg_eps2)
835 END IF
836!
837! Determine root of PDE used for convergence.
838!
839 IF (sg_epsilon.ne.0.0_r8) THEN
840 sg_z2p=sg_z1p/sg_epsilon
841 sg_ror=sg_ro/sg_zrozn
842 sg_zroz1=1.0_r8/(sg_alpha*sg_kappa*sg_ror)
843 sg_zroz2=sg_epsilon*sg_zroz1
844 sg_z1ozn=sg_alpha*sg_kappa*sg_ro
845 sg_z2ozn=sg_z1ozn/sg_epsilon
846!
847 IF ((sg_zroz2.gt.1.0_r8).and.(sg_z1ozn.gt.1.0_r8)) THEN
848 sg_fofx=-sg_ubouc+sg_ubokur*sg_epsilon* &
849 & (log(sg_zroz2)+1.0_r8-sg_epsilon+ &
850 & sg_epsilon*log(sg_z1ozn))
851 ELSE IF ((sg_zroz2.le.1.0_r8).and.(sg_zroz1.gt.1.0_r8).and. &
852 & (sg_z1ozn.gt.1.0_r8)) THEN
853 sg_fofx=-sg_ubouc+sg_ubokur*sg_epsilon*sg_epsilon* &
854 & (sg_zroz1-1.0_r8+log(sg_z1ozn))
855 ELSE IF ((sg_zroz1.le.1.0_r8).and.(sg_z1ozn.gt.1.0_r8)) THEN
856 sg_fofx=-sg_ubouc+sg_ubokur*sg_epsilon*sg_epsilon* &
857 & log(sg_zrozn)
858 ELSE IF ((sg_zroz2.gt.1.0_r8).and.(sg_z1ozn.le.1.0_r8).and. &
859 & (sg_z2ozn.gt.1.0_r8)) THEN
860 sg_fofx=-sg_ubouc+sg_ubokur*sg_epsilon* &
861 & (log(sg_zroz2)+1.0_r8-1.0_r8/sg_z2ozn)
862 ELSE IF ((sg_zroz2.le.1.0_r8).and.(sg_zroz1.gt.1.0_r8).and. &
863 & (sg_z1ozn.le.1.0_r8).and.(sg_z2ozn.gt.1.0_r8)) THEN
864 sg_fofx=-sg_ubouc+sg_ubokur*sg_epsilon*sg_epsilon* &
865 & (sg_zroz1-1.0_r8/sg_z1ozn)
866 ELSE IF ((sg_zroz2.gt.1.0_r8).and.(sg_z2ozn.le.1.0_r8)) THEN
867 sg_fofx=-sg_ubouc+sg_ubokur*sg_epsilon*log(sg_zrozn)
868 END IF
869 END IF
870!
871 RETURN
872 END SUBROUTINE sg_bstress
873!
874 SUBROUTINE sg_purewave (sg_row, sg_ubouwm, sg_znotp, sg_ro)
875!
876!=======================================================================
877! !
878! This routine determines the maximum ratio of waves over combined !
879! bottom shear stress. !
880! !
881! On Input: !
882! !
883! sg_row Ratio of wave excursion amplitude over roughness. !
884! sg_ubouwm Maximum ratio of waves over combined bottom shear !
885! stress. !
886! !
887! On Output: !
888! !
889! sg_ubouwm Maximum ratio of waves over combined bottom shear !
890! stress. !
891! sg_znotp Ratio of hydraulic roughness over scaled height !
892! of bottom boundary layer. !
893! sg_ro Internal friction Rossby number: !
894! ustarc/(omega*znot) !
895! !
896!=======================================================================
897!
898 USE mod_param
899 USE mod_scalars
900!
901! Imported variable declarations.
902!
903 real(r8), intent(in) :: sg_row
904
905 real(r8), intent(inout) :: sg_ubouwm
906
907 real(r8), intent(out) :: sg_znotp, sg_ro
908!
909! Local variable declarations.
910!
911 integer :: Iter
912
913 real(r8) :: cff, sg_bei, sg_beip, sg_ber, sg_berp, sg_kei
914 real(r8) :: sg_keip, sg_ker, sg_kerp, sg_phi, sg_ubouwmn, sg_x
915
916 complex(c8) :: sg_argi, sg_bnot, sg_bnotp, sg_b1, sg_b1p
917 complex(c8) :: sg_gammai, sg_knot, sg_knotp, sg_k1, sg_k1p
918 complex(c8) :: sg_ll, sg_nn
919!
920!-----------------------------------------------------------------------
921! Compute wind-induced wave stress.
922!-----------------------------------------------------------------------
923!
924 DO iter=1,sg_n
925 sg_ro=sg_row/sg_ubouwm
926 sg_znotp=1.0_r8/(sg_kappa*sg_ro)
927 IF (sg_z1p/sg_znotp.gt.1.0_r8) THEN
928 sg_x=2.0_r8*sqrt(sg_znotp)
929 IF (sg_x.le.8.0_r8) THEN
930 CALL sg_kelvin8m (sg_x, sg_ber, sg_bei, sg_ker, sg_kei, &
931 & sg_berp, sg_beip, sg_kerp, sg_keip)
932 ELSE
933 CALL sg_kelvin8p (sg_x, sg_ker, sg_kei, sg_ber, sg_bei, &
934 & sg_kerp, sg_keip, sg_berp, sg_beip)
935 END IF
936 cff=1.0_r8/sqrt(sg_znotp)
937 sg_bnot =cmplx(sg_ber,sg_bei)
938 sg_knot =cmplx(sg_ker,sg_kei)
939 sg_bnotp=cmplx(sg_berp,sg_beip)*cff
940 sg_knotp=cmplx(sg_kerp,sg_keip)*cff
941!
942 sg_x=2.0*sqrt(sg_z1p)
943 IF (sg_x.le.8.0_r8) THEN
944 CALL sg_kelvin8m (sg_x, sg_ber, sg_bei, sg_ker, sg_kei, &
945 & sg_berp, sg_beip, sg_kerp, sg_keip)
946 ELSE
947 CALL sg_kelvin8p (sg_x, sg_ker, sg_kei, sg_ber, sg_bei, &
948 & sg_kerp, sg_keip, sg_berp, sg_beip)
949 END IF
950 cff=1.0_r8/sqrt(sg_z1p)
951 sg_b1 =cmplx(sg_ber,sg_bei)
952 sg_k1 =cmplx(sg_ker,sg_kei)
953 sg_b1p=cmplx(sg_berp,sg_beip)*cff
954 sg_k1p=cmplx(sg_kerp,sg_keip)*cff
955!
956 sg_ll=sg_mp*sg_b1+sg_b1p
957 sg_nn=sg_mp*sg_k1+sg_k1p
958 sg_argi=sg_bnotp*sg_nn/(sg_bnot*sg_nn-sg_knot*sg_ll)+ &
959 & sg_knotp*sg_ll/(sg_knot*sg_ll-sg_bnot*sg_nn)
960 sg_gammai=-sg_kappa*sg_znotp*sg_argi
961 sg_phi=cabs(sg_gammai)
962 ELSE
963 sg_gammai=-sg_kappa*sg_z1p*sg_mp
964 sg_phi=cabs(sg_gammai)
965 END IF
966!
967 sg_ubouwmn=1.0_r8/sg_phi
968 IF (abs((sg_ubouwmn-sg_ubouwm)/sg_ubouwmn).le.sg_tol) THEN
969 sg_ubouwm=sg_ubouwmn
970 RETURN
971 ELSE
972 sg_ubouwm=sg_ubouwmn
973 END IF
974 END DO
975!
976 RETURN
977 END SUBROUTINE sg_purewave
978!
979 SUBROUTINE sg_kelvin8m (x, ber, bei, ker, kei, berp, beip, &
980 & kerp, keip)
981!
982!=======================================================================
983! !
984! This rotuine computes the Kelvin functions for arguments less !
985! than eight. !
986! !
987!=======================================================================
988!
989 USE mod_scalars
990!
991! Imported variable declarations.
992!
993 real(r8), intent(in) :: x
994 real(r8), intent(out) :: ber, bei, ker, kei
995 real(r8), intent(out) :: berp, beip, kerp, keip
996!
997! Local variable declarations.
998!
999 integer :: i
1000
1001 real(r8) :: cff, xhalf
1002
1003 real(r8), dimension(28) :: xp
1004!
1005!-----------------------------------------------------------------------
1006! Compute Kelvin functions.
1007!-----------------------------------------------------------------------
1008!
1009 cff=0.125_r8*x
1010 xp(1)=cff
1011 DO i=2,28
1012 xp(i)=xp(i-1)*cff
1013 END DO
1014 xhalf=0.5_r8*x
1015!
1016 ber=1.0_r8- &
1017 & 64.0_r8*xp(4)+113.77777774_r8*xp(8)- &
1018 & 32.36345652_r8*xp(12)+2.64191397_r8*xp(16)- &
1019 & 0.08349609_r8*xp(20)+0.00122552_r8*xp(24)- &
1020 & 0.00000901_r8*xp(28)
1021 bei=16.0_r8*xp(2)-113.77777774_r8*xp(6)+ &
1022 & 72.81777742*xp(10)-10.56765779_r8*xp(14)+ &
1023 & 0.52185615_r8*xp(18)-0.01103667_r8*xp(22)+ &
1024 & 0.00011346*xp(26)
1025!
1026 ker=-ber*log(xhalf)+0.25_r8*pi*bei- &
1027 & 0.57721566_r8-59.05819744*xp(4)+ &
1028 & 171.36272133_r8*xp(8)-60.60977451_r8*xp(12)+ &
1029 & 5.65539121_r8*xp(16)-0.19636347_r8*xp(20)+ &
1030 & 0.00309699_r8*xp(24)-0.00002458_r8*xp(28)
1031 kei=-bei*log(xhalf)-0.25_r8*pi*ber+ &
1032 & 6.76454936_r8*xp(2)-142.91827687_r8*xp(6)+ &
1033 & 124.23569650_r8*xp(10)-21.30060904_r8*xp(14)+ &
1034 & 1.17509064_r8*xp(18)-0.02695875_r8*xp(22)+ &
1035 & 0.00029532_r8*xp(26)
1036!
1037 berp=x*(-4.0_r8*xp(2)+14.22222222_r8*xp(6)- &
1038 & 6.06814810_r8*xp(10)+0.66047849_r8*xp(14)- &
1039 & 0.02609253_r8*xp(18)+0.00045957_r8*xp(22)- &
1040 & 0.00000394_r8*xp(26))
1041 beip=x*(0.5_r8-10.66666666_r8*xp(4)+11.37777772_r8*xp(8)- &
1042 & 2.31167514_r8*xp(12)+0.14677204_r8*xp(16)- &
1043 & 0.00379386_r8*xp(20)+0.00004609_r8*xp(24))
1044!
1045 kerp=-berp*log(xhalf)-ber/x+0.25*pi*beip+ &
1046 & x*(-3.69113734_r8*xp(2)+21.42034017_r8*xp(6)- &
1047 & 11.36433272_r8*xp(10)+1.41384780_r8*xp(14)- &
1048 & 0.06136358_r8*xp(18)+0.00116137_r8*xp(22)- &
1049 & 0.00001075*xp(26))
1050 keip=-beip*log(xhalf)-bei/x-0.25_r8*pi*berp+ &
1051 & x*(0.21139217_r8-13.39858846_r8*xp(4)+ &
1052 & 19.41182758_r8*xp(8)-4.65950823_r8*xp(12)+ &
1053 & 0.33049424_r8*xp(16)-0.00926707_r8*xp(20)+ &
1054 & 0.00011997_r8*xp(24))
1055!
1056 RETURN
1057 END SUBROUTINE sg_kelvin8m
1058!
1059 SUBROUTINE sg_kelvin8p (x, ker, kei, ber, bei, kerp, keip, &
1060 & berp, beip)
1061!
1062!=======================================================================
1063! !
1064! This rotuine computes the Kelvin functions for arguments greater !
1065! than eight. !
1066! !
1067!=======================================================================
1068!
1069 USE mod_scalars
1070!
1071! Imported variable declarations.
1072!
1073 real(r8), intent(in) :: x
1074 real(r8), intent(out) :: ker, kei, ber, bei
1075 real(r8), intent(out) :: kerp, keip, berp, beip
1076!
1077! Local variable declarations.
1078!
1079 integer :: i
1080
1081 real(r8) :: cff, xhalf
1082
1083 real(r8), dimension(6) :: xm, xp
1084
1085 complex(c8) :: argm, argp, fofx, gofx, phim, phip, thetam, thetap
1086!
1087!-----------------------------------------------------------------------
1088! Compute Kelvin functions.
1089!-----------------------------------------------------------------------
1090!
1091 cff=8.0_r8/x
1092 xp(1)=cff
1093 xm(1)=-cff
1094 DO i=2,6
1095 xp(i)=xp(i-1)*cff
1096 xm(i)=-xm(i-1)*cff
1097 END DO
1098!
1099 thetap=cmplx(0.0_r8,-0.3926991_r8)+ &
1100 & cmplx(0.0110486_r8,-0.0110485_r8)*xp(1)+ &
1101 & cmplx(0.0_r8,-0.0009765_r8)*xp(2)+ &
1102 & cmplx(-0.0000906_r8,-0.0000901_r8)*xp(3)+ &
1103 & cmplx(-0.0000252_r8,0.0_r8)*xp(4)+ &
1104 & cmplx(-0.0000034_r8,0.0000051_r8)*xp(5)+ &
1105 & cmplx(0.0000006,0.0000019)*xp(6)
1106 thetam=cmplx(0.0_r8,-0.3926991_r8)+ &
1107 & cmplx(0.0110486_r8,-0.0110485_r8)*xm(1)+ &
1108 & cmplx(0.0_r8,-0.0009765_r8)*xm(2)+ &
1109 & cmplx(-0.0000906_r8,-0.0000901_r8)*xm(3)+ &
1110 & cmplx(-0.0000252_r8,0.0_r8)*xm(4)+ &
1111 & cmplx(-0.0000034_r8,0.0000051_r8)*xm(5)+ &
1112 & cmplx(0.0000006_r8,0.0000019_r8)*xm(6)
1113!
1114 phip=cmplx(0.7071068_r8,0.7071068_r8)+ &
1115 & cmplx(-0.0625001_r8,-0.0000001_r8)*xp(1)+ &
1116 & cmplx(-0.0013813_r8,0.0013811_r8)*xp(2)+ &
1117 & cmplx(0.0000005_r8,0.0002452_r8)*xp(3)+ &
1118 & cmplx(0.0000346_r8,0.0000338_r8)*xp(4)+ &
1119 & cmplx(0.0000117_r8,-0.0000024_r8)*xp(5)+ &
1120 & cmplx(0.0000016_r8,-0.0000032_r8)*xp(6)
1121 phim=cmplx(0.7071068_r8,0.7071068_r8)+ &
1122 & cmplx(-0.0625001_r8,-0.0000001_r8)*xm(1)+ &
1123 & cmplx(-0.0013813_r8,0.0013811_r8)*xm(2)+ &
1124 & cmplx(0.0000005_r8,0.0002452_r8)*xm(3)+ &
1125 & cmplx(0.0000346_r8,0.0000338_r8)*xm(4)+ &
1126 & cmplx(0.0000117_r8,-0.0000024_r8)*xm(5)+ &
1127 & cmplx(0.0000016_r8,-0.0000032_r8)*xm(6)
1128!
1129 cff=x/sqrt(2.0_r8)
1130 argm=-cff*cmplx(1.0_r8,1.0_r8)+thetam
1131 fofx=sqrt(pi/(2.0_r8*x))*cexp(argm)
1132 ker=real(fofx)
1133 kei=aimag(fofx)
1134!
1135 argp=cff*cmplx(1.0_r8,1.0_r8)+thetap
1136 gofx=1.0_r8/sqrt(2.0_r8*pi*x)*cexp(argp)
1137 ber=real(gofx)-kei/pi
1138 bei=aimag(gofx)+ker/pi
1139!
1140 kerp=real(-fofx*phim)
1141 keip=aimag(-fofx*phim)
1142!
1143 berp=real(gofx*phip)-keip/pi
1144 beip=aimag(gofx*phip)+kerp/pi
1145!
1146 RETURN
1147 END SUBROUTINE sg_kelvin8p
1148
1149 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::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::sg_bbl_tile
subroutine sg_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, rho, u, v, bottom, iconv, ubot, vbot, ur, vr, bustrc, bvstrc, bustrw, bvstrw, bustrcwmax, bvstrcwmax, bustr, bvstr)
Definition sg_bbl.h:116
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
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_ocean
Definition mod_ocean.F:2
mod_ocean::ocean
type(t_ocean), dimension(:), allocatable ocean
Definition mod_ocean.F:351
mod_param
Definition mod_param.F:2
mod_param::inlm
integer, parameter inlm
Definition mod_param.F:662
mod_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_g
real(dp) sg_g
Definition mod_scalars.F:1540
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::sg_nu
real(dp) sg_nu
Definition mod_scalars.F:1542
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::irlen
integer, parameter irlen
Definition sediment_mod.h:128
mod_sediment::irhgt
integer, parameter irhgt
Definition sediment_mod.h:129
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::izapp
integer, parameter izapp
Definition sediment_mod.h:132
mod_sediment::izdef
integer, parameter izdef
Definition sediment_mod.h:131
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