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