| blob | c81e67c33eaac0873f3d3d9b43e213fe70e84eb3 |
1 !WRF:model_layer:physics
2 !
3 module module_bl_gwdo
4 contains
5 !-------------------------------------------------------------------------------
6 subroutine gwdo(u3d,v3d,t3d,qv3d,p3d,p3di,pi3d,z, &
7 rublten,rvblten, &
8 dtaux3d,dtauy3d,dusfcg,dvsfcg, &
9 var2d,oc12d,oa2d1,oa2d2,oa2d3,oa2d4,ol2d1,ol2d2,ol2d3,ol2d4, &
10 sina,cosa,znu,znw,p_top, &
11 cp,g,rd,rv,ep1,pi, &
12 dt,dx,kpbl2d,itimestep, &
13 ids,ide, jds,jde, kds,kde, &
14 ims,ime, jms,jme, kms,kme, &
15 its,ite, jts,jte, kts,kte)
16 !-------------------------------------------------------------------------------
17 implicit none
18 !-------------------------------------------------------------------------------
19 !
20 !-- u3d 3d u-velocity interpolated to theta points (m/s)
21 !-- v3d 3d v-velocity interpolated to theta points (m/s)
22 !-- t3d temperature (k)
23 !-- qv3d 3d water vapor mixing ratio (kg/kg)
24 !-- p3d 3d pressure (pa)
25 !-- p3di 3d pressure (pa) at interface level
26 !-- pi3d 3d exner function (dimensionless)
27 !-- rublten u tendency due to pbl parameterization (m/s/s)
28 !-- rvblten v tendency due to pbl parameterization (m/s/s)
29 !-- sina sine rotation angle
30 !-- cosa cosine rotation angle
31 !-- znu eta values (sigma values)
32 !-- cp heat capacity at constant pressure for dry air (j/kg/k)
33 !-- g acceleration due to gravity (m/s^2)
34 !-- rd gas constant for dry air (j/kg/k)
35 !-- z height above sea level (m)
36 !-- rv gas constant for water vapor (j/kg/k)
37 !-- dt time step (s)
38 !-- dx model grid interval (m)
39 !-- ep1 constant for virtual temperature (r_v/r_d - 1) (dimensionless)
40 !-- ids start index for i in domain
41 !-- ide end index for i in domain
42 !-- jds start index for j in domain
43 !-- jde end index for j in domain
44 !-- kds start index for k in domain
45 !-- kde end index for k in domain
46 !-- ims start index for i in memory
47 !-- ime end index for i in memory
48 !-- jms start index for j in memory
49 !-- jme end index for j in memory
50 !-- kms start index for k in memory
51 !-- kme end index for k in memory
52 !-- its start index for i in tile
53 !-- ite end index for i in tile
54 !-- jts start index for j in tile
55 !-- jte end index for j in tile
56 !-- kts start index for k in tile
57 !-- kte end index for k in tile
58 !
59 !-------------------------------------------------------------------------------
60 integer, intent(in ) :: ids,ide, jds,jde, kds,kde, &
61 ims,ime, jms,jme, kms,kme, &
62 its,ite, jts,jte, kts,kte
63 integer, intent(in ) :: itimestep
64 !
65 real, intent(in ) :: dt,dx,cp,g,rd,rv,ep1,pi
66 !
67 real, dimension( ims:ime, kms:kme, jms:jme ) , &
68 intent(in ) :: qv3d, &
69 p3d, &
70 pi3d, &
71 t3d, &
72 z
73 real, dimension( ims:ime, kms:kme, jms:jme ) , &
74 intent(in ) :: p3di
75 !
76 real, dimension( ims:ime, kms:kme, jms:jme ) , &
77 intent(inout) :: rublten, &
78 rvblten
79 real, dimension( ims:ime, kms:kme, jms:jme ) , &
80 intent(inout) :: dtaux3d, &
81 dtauy3d
82 !
83 real, dimension( ims:ime, kms:kme, jms:jme ) , &
84 intent(in ) :: u3d, &
85 v3d
86 !
87 integer, dimension( ims:ime, jms:jme ) , &
88 intent(in ) :: kpbl2d
89 real, dimension( ims:ime, jms:jme ) , &
90 intent(inout ) :: dusfcg, &
91 dvsfcg
92 !
93 real, dimension( ims:ime, jms:jme ) , &
94 intent(in ) :: var2d, &
95 oc12d, &
96 oa2d1,oa2d2,oa2d3,oa2d4, &
97 ol2d1,ol2d2,ol2d3,ol2d4, &
98 sina,cosa
99 !
100 real, dimension( kms:kme ) , &
101 optional , &
102 intent(in ) :: znu, &
103 znw
104 !
105 real, optional, intent(in ) :: p_top
106 !
107 !local
108 !
109 real, dimension( its:ite, kts:kte ) :: delprsi, &
110 pdh
111 real, dimension( its:ite, kts:kte ) :: ugeo, vgeo, dudt, dvdt, dtaux, dtauy
112 real, dimension( its:ite ) :: dusfc, dvsfc
113 real, dimension( its:ite, kts:kte+1 ) :: pdhi
114 real, dimension( its:ite, 4 ) :: oa4, &
115 ol4
116 integer :: i,j,k,kpblmax
117 !
118 do k = kts,kte
119 if (znu(k).gt.0.6) kpblmax = k + 1
120 enddo
121 !
122 do j = jts,jte
123 do k = kts,kte+1
124 do i = its,ite
125 if (k.le.kte)pdh(i,k) = p3d(i,k,j)
126 pdhi(i,k) = p3di(i,k,j)
127 enddo
128 enddo
129 !
130 do k = kts,kte
131 do i = its,ite
132 delprsi(i,k) = pdhi(i,k)-pdhi(i,k+1)
133 ! rotate winds to zonal/meridional
134 ugeo(i,k) = u3d(i,k,j)*cosa(i,j) - v3d(i,k,j)*sina(i,j)
135 vgeo(i,k) = u3d(i,k,j)*sina(i,j) + v3d(i,k,j)*cosa(i,j)
136 dudt(i,k) = 0.0
137 dvdt(i,k) = 0.0
138 enddo
139 enddo
140 do i = its,ite
141 oa4(i,1) = oa2d1(i,j)
142 oa4(i,2) = oa2d2(i,j)
143 oa4(i,3) = oa2d3(i,j)
144 oa4(i,4) = oa2d4(i,j)
145 ol4(i,1) = ol2d1(i,j)
146 ol4(i,2) = ol2d2(i,j)
147 ol4(i,3) = ol2d3(i,j)
148 ol4(i,4) = ol2d4(i,j)
149 enddo
150 call gwdo2d(dudt=dudt(its,kts),dvdt=dvdt(its,kts) &
151 ,dtaux2d=dtaux(its,kts),dtauy2d=dtauy(its,kts) &
152 ,u1=ugeo(its,kts),v1=vgeo(its,kts) &
153 ,t1=t3d(ims,kms,j),q1=qv3d(ims,kms,j) &
154 ,del=delprsi(its,kts) &
155 ,prsi=pdhi(its,kts) &
156 ,prsl=pdh(its,kts),prslk=pi3d(ims,kms,j) &
157 ,zl=z(ims,kms,j) &
158 ,kpblmax=kpblmax &
159 ,var=var2d(ims,j),oc1=oc12d(ims,j) &
160 ,oa4=oa4,ol4=ol4 &
161 ,dusfc=dusfc(its),dvsfc=dvsfc(its) &
162 ,g_=g,cp_=cp,rd_=rd,rv_=rv,fv_=ep1,pi_=pi &
163 ,dxmeter=dx,deltim=dt &
164 ,kpbl=kpbl2d(ims,j),lat=j &
165 ,ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde &
166 ,ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme &
167 ,its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte )
168 do k = kts,kte
169 do i = its,ite
170 ! rotate tendencies from zonal/meridional to model grid
171 rublten(i,k,j) = rublten(i,k,j)+dudt(i,k)*cosa(i,j) + dvdt(i,k)*sina(i,j)
172 rvblten(i,k,j) = rvblten(i,k,j)-dudt(i,k)*sina(i,j) + dvdt(i,k)*cosa(i,j)
173 dtaux3d(i,k,j) = dtaux(i,k)*cosa(i,j) + dtauy(i,k)*sina(i,j)
174 dtauy3d(i,k,j) =-dtaux(i,k)*sina(i,j) + dtauy(i,k)*cosa(i,j)
175 if(k.eq.kts)then
176 dusfcg(i,j) = dusfc(i)*cosa(i,j) + dvsfc(i)*sina(i,j)
177 dvsfcg(i,j) =-dusfc(i)*sina(i,j) + dvsfc(i)*cosa(i,j)
178 endif
179 enddo
180 enddo
181 enddo
182 !
183 end subroutine gwdo
184 !-------------------------------------------------------------------------------
185 !-------------------------------------------------------------------------------
186 subroutine gwdo2d(dudt, dvdt, dtaux2d, dtauy2d, &
187 u1, v1, t1, q1, &
188 del, &
189 prsi, prsl, prslk, zl, &
190 kpblmax, &
191 var, oc1, oa4, ol4, dusfc, dvsfc, &
192 g_, cp_, rd_, rv_, fv_, pi_, &
193 dxmeter, deltim, kpbl, lat, &
194 ids, ide, jds, jde, kds, kde, &
195 ims, ime, jms, jme, kms, kme, &
196 its, ite, jts, jte, kts, kte)
197 !-------------------------------------------------------------------------------
198 !
199 ! abstract :
200 ! this code handles the time tendencies of u v due to the effect of
201 ! mountain induced gravity wave drag from sub-grid scale orography.
202 ! this routine not only treats the traditional upper-level wave breaking due
203 ! to mountain variance (alpert 1988), but also the enhanced
204 ! lower-tropospheric wave breaking due to mountain convexity and asymmetry
205 ! (kim and arakawa 1995). thus, in addition to the terrain height data
206 ! in a model grid gox, additional 10-2d topographic statistics files are
207 ! needed, including orographic standard deviation (var), convexity (oc1),
208 ! asymmetry (oa4) and ol (ol4). these data sets are prepared based on the
209 ! 30 sec usgs orography (hong 1999). the current scheme was implmented as in
210 ! choi and hong (2015), which names kim gwdo since it was developed by
211 ! kiaps staffs for kiaps integrated model system (kim). the scheme
212 ! additionally includes the effects of orographic anisotropy and
213 ! flow-blocking drag.
214 ! coded by song-you hong and young-joon kim and implemented by song-you hong
215 !
216 ! history log :
217 ! 2015-07-01 hyun-joo choi add flow-blocking drag and orographic anisotropy
218 !
219 ! references :
220 ! choi and hong (2015), j. geophys. res.
221 ! hong et al. (2008), wea. forecasting
222 ! kim and doyle (2005), q. j. r. meteor. soc.
223 ! kim and arakawa (1995), j. atmos. sci.
224 ! alpet et al. (1988), NWP conference
225 ! hong (1999), NCEP office note 424
226 !
227 ! input :
228 ! dudt, dvdt - non-lin tendency for u and v wind component
229 ! u1, v1 - zonal and meridional wind m/sec at t0-dt
230 ! t1 - temperature deg k at t0-dt
231 ! q1 - mixing ratio at t0-dt
232 ! deltim - time step (s)
233 ! del - positive increment of pressure across layer (pa)
234 ! kpblmax, kpbl - vertical index of pbl height
235 ! prslk, zl, prsl, prsi - pressure and height variables
236 ! oa4, ol4, omax, var, oc1 - orographic statistics
237 !
238 ! output :
239 ! dudt, dvdt - wind tendency due to gwdo
240 ! dtaux2d, dtauy2d - diagnoised orographic gwd
241 ! dusfc, dvsfc - gw stress
242 !
243 !-------------------------------------------------------------------------------
244 implicit none
245 !
246 integer , intent(in ) :: lat, kpblmax, &
247 ids, ide, jds, jde, &
248 kds, kde, ims, ime, &
249 jms, jme, kms, kme, &
250 its, ite, jts, jte, &
251 kts, kte
252 integer, dimension(ims:ime) , intent(in ) :: kpbl
253 real , intent(in ) :: g_, pi_, rd_, rv_, fv_,&
254 cp_, deltim
255 real , intent(in ) :: dxmeter
256 real, dimension(its:ite,kts:kte) , intent(inout) :: dudt, dvdt
257 real, dimension(its:ite,kts:kte) , intent( out) :: dtaux2d, dtauy2d
258 real, dimension(its:ite,kts:kte) , intent(in ) :: u1, v1
259 real, dimension(ims:ime,kms:kme) , intent(in ) :: t1, q1, prslk, zl
260 !
261 real, dimension(its:ite,kts:kte) , intent(in ) :: prsl, del
262 real, dimension(its:ite,kts:kte+1), intent(in ) :: prsi
263 real, dimension(its:ite,4) , intent(in ) :: oa4, ol4
264 !
265 real, dimension(ims:ime) , intent(in ) :: var, oc1
266 real, dimension(its:ite) , intent( out) :: dusfc, dvsfc
267 !
268 real, parameter :: ric = 0.25 ! critical richardson number
269 real, parameter :: dw2min = 1.
270 real, parameter :: rimin = -100.
271 real, parameter :: bnv2min = 1.0e-5
272 real, parameter :: efmin = 0.0
273 real, parameter :: efmax = 10.0
274 real, parameter :: xl = 4.0e4
275 real, parameter :: critac = 1.0e-5
276 real, parameter :: gmax = 1.
277 real, parameter :: veleps = 1.0
278 real, parameter :: frc = 1.0
279 real, parameter :: ce = 0.8
280 real, parameter :: cg = 0.5
281 integer,parameter :: kpblmin = 2
282 !
283 ! local variables
284 !
285 integer :: latd,lond
286 integer :: i,k,lcap,lcapp1,nwd,idir, &
287 klcap,kp1,ikount,kk
288 !
289 real :: fdir,cs,rcsks, &
290 wdir,ti,rdz,temp,tem2,dw2,shr2,bvf2,rdelks, &
291 wtkbj,tem,gfobnv,hd,fro,rim,temc,tem1,efact, &
292 temv,dtaux,dtauy
293 !
294 logical, dimension(its:ite) :: ldrag, icrilv, flag,kloop1
295 real, dimension(its:ite) :: coefm
296 !
297 real, dimension(its:ite) :: taub, xn, yn, ubar, vbar, fr, &
298 ulow, rulow, bnv, oa, ol, rhobar, &
299 dtfac, brvf, xlinv, delks,delks1, &
300 zlowtop,cleff
301 real, dimension(its:ite,kts:kte+1) :: taup
302 real, dimension(its:ite,kts:kte-1) :: velco
303 real, dimension(its:ite,kts:kte) :: bnv2, usqj, taud, rho, vtk, vtj
304 !
305 integer, dimension(its:ite) :: kbl, klowtop
306 integer, parameter :: mdir=8
307 integer, dimension(mdir) :: nwdir
308 data nwdir/6,7,5,8,2,3,1,4/
309 !
310 ! variables for flow-blocking drag
311 !
312 real, parameter :: frmax = 10.
313 real, parameter :: olmin = 1.0e-5
314 real, parameter :: odmin = 0.1
315 real, parameter :: odmax = 10.
316 !
317 real :: fbdcd
318 real :: zblk, tautem
319 real :: fbdpe, fbdke
320 real, dimension(its:ite) :: delx, dely
321 real, dimension(its:ite,4) :: dxy4, dxy4p
322 real, dimension(4) :: ol4p
323 real, dimension(its:ite) :: dxy, dxyp, olp, od
324 real, dimension(its:ite,kts:kte+1) :: taufb
325 !
326 integer, dimension(its:ite) :: komax
327 integer :: kblk
328 !-------------------------------------------------------------------------------
329 !
330 ! constants
331 !
332 lcap = kte
333 lcapp1 = lcap + 1
334 fdir = mdir / (2.0*pi_)
335 !
336 ! calculate length of grid for flow-blocking drag
337 !
338 delx(its:ite) = dxmeter
339 dely(its:ite) = dxmeter
340 dxy4(its:ite,1) = delx(its:ite)
341 dxy4(its:ite,2) = dely(its:ite)
342 dxy4(its:ite,3) = sqrt(delx(its:ite)**2. + dely(its:ite)**2.)
343 dxy4(its:ite,4) = dxy4(its:ite,3)
344 dxy4p(its:ite,1) = dxy4(its:ite,2)
345 dxy4p(its:ite,2) = dxy4(its:ite,1)
346 dxy4p(its:ite,3) = dxy4(its:ite,4)
347 dxy4p(its:ite,4) = dxy4(its:ite,3)
348 !
349 cleff(its:ite) = dxmeter
350 !
351 ! initialize arrays
352 !
353 ldrag = .false. ; icrilv = .false. ; flag = .true.
354 !
355 klowtop = 0 ; kbl = 0
356 !
357 dtaux = 0. ; dtauy = 0. ; xn = 0. ; yn = 0.
358 ubar = 0. ; vbar = 0. ; rhobar = 0. ; ulow = 0.
359 oa = 0. ; ol = 0. ; taub = 0.
360 !
361 usqj = 0. ; bnv2 = 0. ; vtj = 0. ; vtk = 0.
362 taup = 0. ; taud = 0. ; dtaux2d = 0. ; dtauy2d = 0.
363 !
364 dtfac = 1.0 ; xlinv = 1.0/xl
365 !
366 ! initialize arrays for flow-blocking drag
367 !
368 komax = 0
369 taufb = 0.0
370 !
371 do k = kts,kte
372 do i = its,ite
373 vtj(i,k) = t1(i,k) * (1.+fv_*q1(i,k))
374 vtk(i,k) = vtj(i,k) / prslk(i,k)
375 rho(i,k) = 1./rd_ * prsl(i,k) / vtj(i,k) ! density kg/m**3
376 enddo
377 enddo
378 !
379 do i = its,ite
380 zlowtop(i) = 2. * var(i)
381 enddo
382 !
383 do i = its,ite
384 kloop1(i) = .true.
385 enddo
386 !
387 do k = kts+1,kte
388 do i = its,ite
389 if (kloop1(i).and.zl(i,k)-zl(i,1).ge.zlowtop(i)) then
390 klowtop(i) = k+1
391 kloop1(i) = .false.
392 endif
393 enddo
394 enddo
395 !
396 do i = its,ite
397 !
398 ! determine reference level: 2*var
399 !
400 kbl(i) = klowtop(i)
401 kbl(i) = max(min(kbl(i),kpblmax),kpblmin)
402 enddo
403 !
404 ! determine the level of maximum orographic height
405 !
406 komax(:) = kbl(:)
407 !
408 do i = its,ite
409 delks(i) = 1.0 / (prsi(i,1) - prsi(i,kbl(i)))
410 delks1(i) = 1.0 / (prsl(i,1) - prsl(i,kbl(i)))
411 enddo
412 !
413 ! compute low level averages within pbl
414 !
415 do k = kts,kpblmax
416 do i = its,ite
417 if (k.lt.kbl(i)) then
418 rcsks = del(i,k) * delks(i)
419 rdelks = del(i,k) * delks(i)
420 ubar(i) = ubar(i) + rcsks * u1(i,k) ! pbl u mean
421 vbar(i) = vbar(i) + rcsks * v1(i,k) ! pbl v mean
422 rhobar(i) = rhobar(i) + rdelks * rho(i,k) ! pbl rho mean
423 endif
424 enddo
425 enddo
426 !
427 ! figure out low-level horizontal wind direction
428 !
429 ! nwd 1 2 3 4 5 6 7 8
430 ! wd w s sw nw e n ne se
431 !
432 do i = its,ite
433 wdir = atan2(ubar(i),vbar(i)) + pi_
434 idir = mod(nint(fdir*wdir),mdir) + 1
435 nwd = nwdir(idir)
436 oa(i) = (1-2*int( (nwd-1)/4 )) * oa4(i,mod(nwd-1,4)+1)
437 ol(i) = ol4(i,mod(nwd-1,4)+1)
438 !
439 ! compute orographic width along (ol) and perpendicular (olp) the wind direction
440 !
441 ol4p(1) = ol4(i,2)
442 ol4p(2) = ol4(i,1)
443 ol4p(3) = ol4(i,4)
444 ol4p(4) = ol4(i,3)
445 olp(i) = ol4p(mod(nwd-1,4)+1)
446 !
447 ! compute orographic direction (horizontal orographic aspect ratio)
448 !
449 od(i) = olp(i)/max(ol(i),olmin)
450 od(i) = min(od(i),odmax)
451 od(i) = max(od(i),odmin)
452 !
453 ! compute length of grid in the along(dxy) and cross(dxyp) wind directions
454 !
455 dxy(i) = dxy4(i,MOD(nwd-1,4)+1)
456 dxyp(i) = dxy4p(i,MOD(nwd-1,4)+1)
457 enddo
458 !
459 ! saving richardson number in usqj for migwdi
460 !
461 do k = kts,kte-1
462 do i = its,ite
463 ti = 2.0 / (t1(i,k)+t1(i,k+1))
464 rdz = 1./(zl(i,k+1) - zl(i,k))
465 tem1 = u1(i,k) - u1(i,k+1)
466 tem2 = v1(i,k) - v1(i,k+1)
467 dw2 = tem1*tem1 + tem2*tem2
468 shr2 = max(dw2,dw2min) * rdz * rdz
469 bvf2 = g_*(g_/cp_+rdz*(vtj(i,k+1)-vtj(i,k))) * ti
470 usqj(i,k) = max(bvf2/shr2,rimin)
471 bnv2(i,k) = 2.0*g_*rdz*(vtk(i,k+1)-vtk(i,k))/(vtk(i,k+1)+vtk(i,k))
472 enddo
473 enddo
474 !
475 ! compute the "low level" or 1/3 wind magnitude (m/s)
476 !
477 do i = its,ite
478 ulow(i) = max(sqrt(ubar(i)*ubar(i) + vbar(i)*vbar(i)), 1.0)
479 rulow(i) = 1./ulow(i)
480 enddo
481 !
482 do k = kts,kte-1
483 do i = its,ite
484 velco(i,k) = 0.5 * ((u1(i,k)+u1(i,k+1)) * ubar(i) &
485 + (v1(i,k)+v1(i,k+1)) * vbar(i))
486 velco(i,k) = velco(i,k) * rulow(i)
487 if ((velco(i,k).lt.veleps) .and. (velco(i,k).gt.0.)) then
488 velco(i,k) = veleps
489 endif
490 enddo
491 enddo
492 !
493 ! no drag when critical level in the base layer
494 !
495 do i = its,ite
496 ldrag(i) = velco(i,1).le.0.
497 enddo
498 !
499 ! no drag when velco.lt.0
500 !
501 do k = kpblmin,kpblmax
502 do i = its,ite
503 if (k .lt. kbl(i)) ldrag(i) = ldrag(i).or. velco(i,k).le.0.
504 enddo
505 enddo
506 !
507 ! the low level weighted average ri is stored in usqj(1,1; im)
508 ! the low level weighted average n**2 is stored in bnv2(1,1; im)
509 ! this is called bnvl2 in phy_gwd_alpert_sub not bnv2
510 ! rdelks (del(k)/delks) vert ave factor so we can * instead of /
511 !
512 do i = its,ite
513 wtkbj = (prsl(i,1)-prsl(i,2)) * delks1(i)
514 bnv2(i,1) = wtkbj * bnv2(i,1)
515 usqj(i,1) = wtkbj * usqj(i,1)
516 enddo
517 !
518 do k = kpblmin,kpblmax
519 do i = its,ite
520 if (k .lt. kbl(i)) then
521 rdelks = (prsl(i,k)-prsl(i,k+1)) * delks1(i)
522 bnv2(i,1) = bnv2(i,1) + bnv2(i,k) * rdelks
523 usqj(i,1) = usqj(i,1) + usqj(i,k) * rdelks
524 endif
525 enddo
526 enddo
527 !
528 do i = its,ite
529 ldrag(i) = ldrag(i) .or. bnv2(i,1).le.0.0
530 ldrag(i) = ldrag(i) .or. ulow(i).eq.1.0
531 ldrag(i) = ldrag(i) .or. var(i) .le. 0.0
532 enddo
533 !
534 ! set all ri low level values to the low level value
535 !
536 do k = kpblmin,kpblmax
537 do i = its,ite
538 if (k .lt. kbl(i)) usqj(i,k) = usqj(i,1)
539 enddo
540 enddo
541 !
542 do i = its,ite
543 if (.not.ldrag(i)) then
544 bnv(i) = sqrt( bnv2(i,1) )
545 fr(i) = bnv(i) * rulow(i) * var(i) * od(i)
546 fr(i) = min(fr(i),frmax)
547 xn(i) = ubar(i) * rulow(i)
548 yn(i) = vbar(i) * rulow(i)
549 endif
550 enddo
551 !
552 ! compute the base level stress and store it in taub
553 ! calculate enhancement factor, number of mountains & aspect
554 ! ratio const. use simplified relationship between standard
555 ! deviation & critical hgt
556 !
557 do i = its,ite
558 if (.not. ldrag(i)) then
559 efact = (oa(i) + 2.) ** (ce*fr(i)/frc)
560 efact = min( max(efact,efmin), efmax )
561 coefm(i) = (1. + ol(i)) ** (oa(i)+1.)
562 xlinv(i) = coefm(i) / cleff(i)
563 tem = fr(i) * fr(i) * oc1(i)
564 gfobnv = gmax * tem / ((tem + cg)*bnv(i))
565 taub(i) = xlinv(i) * rhobar(i) * ulow(i) * ulow(i) &
566 * ulow(i) * gfobnv * efact
567 else
568 taub(i) = 0.0
569 xn(i) = 0.0
570 yn(i) = 0.0
571 endif
572 enddo
573 !
574 ! now compute vertical structure of the stress.
575 !
576 do k = kts,kpblmax
577 do i = its,ite
578 if (k .le. kbl(i)) taup(i,k) = taub(i)
579 enddo
580 enddo
581 !
582 do k = kpblmin, kte-1 ! vertical level k loop!
583 kp1 = k + 1
584 do i = its,ite
585 !
586 ! unstablelayer if ri < ric
587 ! unstable layer if upper air vel comp along surf vel <=0 (crit lay)
588 ! at (u-c)=0. crit layer exists and bit vector should be set (.le.)
589 !
590 if (k .ge. kbl(i)) then
591 icrilv(i) = icrilv(i) .or. ( usqj(i,k) .lt. ric) &
592 .or. (velco(i,k) .le. 0.0)
593 brvf(i) = max(bnv2(i,k),bnv2min) ! brunt-vaisala frequency squared
594 brvf(i) = sqrt(brvf(i)) ! brunt-vaisala frequency
595 endif
596 enddo
597 !
598 do i = its,ite
599 if (k .ge. kbl(i) .and. (.not. ldrag(i))) then
600 if (.not.icrilv(i) .and. taup(i,k) .gt. 0.0 ) then
601 temv = 1.0 / velco(i,k)
602 tem1 = coefm(i)/dxy(i)*(rho(i,kp1)+rho(i,k))*brvf(i)*velco(i,k)*0.5
603 hd = sqrt(taup(i,k) / tem1)
604 fro = brvf(i) * hd * temv
605 !
606 ! rim is the minimum-richardson number by shutts (1985)
607 !
608 tem2 = sqrt(usqj(i,k))
609 tem = 1. + tem2 * fro
610 rim = usqj(i,k) * (1.-fro) / (tem * tem)
611 !
612 ! check stability to employ the 'saturation hypothesis'
613 ! of lindzen (1981) except at tropospheric downstream regions
614 !
615 if (rim .le. ric) then ! saturation hypothesis!
616 if ((oa(i) .le. 0.).or.(kp1 .ge. kpblmin )) then
617 temc = 2.0 + 1.0 / tem2
618 hd = velco(i,k) * (2.*sqrt(temc)-temc) / brvf(i)
619 taup(i,kp1) = tem1 * hd * hd
620 endif
621 else ! no wavebreaking!
622 taup(i,kp1) = taup(i,k)
623 endif
624 endif
625 endif
626 enddo
627 enddo
628 !
629 if (lcap.lt.kte) then
630 do klcap = lcapp1,kte
631 do i = its,ite
632 taup(i,klcap) = prsi(i,klcap) / prsi(i,lcap) * taup(i,lcap)
633 enddo
634 enddo
635 endif
636 do i = its,ite
637 if (.not.ldrag(i)) then
638 !
639 ! determine the height of flow-blocking layer
640 !
641 kblk = 0
642 fbdpe = 0.0
643 fbdke = 0.0
644 do k = kte, kpblmin, -1
645 if (kblk.eq.0 .and. k.le.kbl(i)) then
646 fbdpe = fbdpe + bnv2(i,k)*(zl(i,kbl(i))-zl(i,k)) &
647 *del(i,k)/g_/rho(i,k)
648 fbdke = 0.5*(u1(i,k)**2.+v1(i,k)**2.)
649 !
650 ! apply flow-blocking drag when fbdpe >= fbdke
651 !
652 if (fbdpe.ge.fbdke) then
653 kblk = k
654 kblk = min(kblk,kbl(i))
655 zblk = zl(i,kblk)-zl(i,kts)
656 endif
657 endif
658 enddo
659 if (kblk.ne.0) then
660 !
661 ! compute flow-blocking stress
662 !
663 fbdcd = max(2.0-1.0/od(i),0.0)
664 taufb(i,kts) = 0.5*rhobar(i)*coefm(i)/dxmeter**2*fbdcd*dxyp(i) &
665 *olp(i)*zblk*ulow(i)**2
666 tautem = taufb(i,kts)/real(kblk-kts)
667 do k = kts+1, kblk
668 taufb(i,k) = taufb(i,k-1) - tautem
669 enddo
670 !
671 ! sum orographic GW stress and flow-blocking stress
672 !
673 taup(i,:) = taup(i,:) + taufb(i,:)
674 endif
675 endif
676 enddo
677 !
678 ! calculate - (g)*d(tau)/d(pressure) and deceleration terms dtaux, dtauy
679 !
680 do k = kts,kte
681 do i = its,ite
682 taud(i,k) = 1. * (taup(i,k+1) - taup(i,k)) * g_ / del(i,k)
683 enddo
684 enddo
685 !
686 ! if the gravity wave drag would force a critical line
687 ! in the lower ksmm1 layers during the next deltim timestep,
688 ! then only apply drag until that critical line is reached.
689 !
690 do k = kts,kpblmax-1
691 do i = its,ite
692 if (k .le. kbl(i)) then
693 if (taud(i,k).ne.0.) &
694 dtfac(i) = min(dtfac(i),abs(velco(i,k)/(deltim*taud(i,k))))
695 endif
696 enddo
697 enddo
698 !
699 do i = its,ite
700 dusfc(i) = 0.
701 dvsfc(i) = 0.
702 enddo
703 !
704 do k = kts,kte
705 do i = its,ite
706 taud(i,k) = taud(i,k) * dtfac(i)
707 dtaux = taud(i,k) * xn(i)
708 dtauy = taud(i,k) * yn(i)
709 dtaux2d(i,k) = dtaux
710 dtauy2d(i,k) = dtauy
711 dudt(i,k) = dtaux + dudt(i,k)
712 dvdt(i,k) = dtauy + dvdt(i,k)
713 dusfc(i) = dusfc(i) + dtaux * del(i,k)
714 dvsfc(i) = dvsfc(i) + dtauy * del(i,k)
715 enddo
716 enddo
717 !
718 do i = its,ite
719 dusfc(i) = (-1./g_) * dusfc(i)
720 dvsfc(i) = (-1./g_) * dvsfc(i)
721 enddo
722 !
723 return
724 end subroutine gwdo2d
725 !-------------------------------------------------------------------------------
726 !-------------------------------------------------------------------------------
727 end module module_bl_gwdo