| blob | 99212b533b344f5e26f499d6c25e2426a0dd230c |
1 !WRF:model_layer:physics
2 !
3 !
4 !
5 !
6 module module_bl_gwdo_gsl
7 contains
8 !-------------------------------------------------------------------------------
9 subroutine gwdo_gsl(u3d,v3d,t3d,qv3d,p3d,p3di,pi3d,z, &
10 rublten,rvblten,rthblten, &
11 dtaux3d_ls,dtauy3d_ls,dtaux3d_bl,dtauy3d_bl, &
12 dtaux3d_ss,dtauy3d_ss,dtaux3d_fd,dtauy3d_fd, &
13 dusfcg_ls,dvsfcg_ls,dusfcg_bl,dvsfcg_bl,dusfcg_ss,dvsfcg_ss, &
14 dusfcg_fd,dvsfcg_fd,xland,br, &
15 var2d,oc12d,oa2d1,oa2d2,oa2d3,oa2d4,ol2d1,ol2d2,ol2d3,ol2d4, &
16 var2dss,oc12dss, &
17 oa2d1ss,oa2d2ss,oa2d3ss,oa2d4ss, &
18 ol2d1ss,ol2d2ss,ol2d3ss,ol2d4ss, &
19 sina,cosa,znu,znw,p_top,dz,pblh, &
20 cp,g,rd,rv,ep1,pi, &
21 dt,dx,kpbl2d,itimestep,gwd_opt,gwd_diags, &
22 spp_pbl,pattern_spp_pbl, &
23 ids,ide, jds,jde, kds,kde, &
24 ims,ime, jms,jme, kms,kme, &
25 its,ite, jts,jte, kts,kte)
26 !-------------------------------------------------------------------------------
27 implicit none
28 !-------------------------------------------------------------------------------
29 !
30 !-- u3d 3d u-velocity interpolated to theta points (m/s)
31 !-- v3d 3d v-velocity interpolated to theta points (m/s)
32 !-- t3d temperature (k)
33 !-- qv3d 3d water vapor mixing ratio (kg/kg)
34 !-- p3d 3d pressure (pa)
35 !-- p3di 3d pressure (pa) at interface level
36 !-- pi3d 3d exner function (dimensionless)
37 !-- rublten u tendency due to pbl parameterization (m/s/s)
38 !-- rvblten v tendency due to pbl parameterization (m/s/s)
39 !-- rthblten theta tendency due to pbl parameterization (K/s)
40 !-- sina sine rotation angle
41 !-- cosa cosine rotation angle
42 !-- znu eta values (sigma values)
43 !-- cp heat capacity at constant pressure for dry air (j/kg/k)
44 !-- g acceleration due to gravity (m/s^2)
45 !-- rd gas constant for dry air (j/kg/k)
46 !-- z height above sea level (m)
47 !-- rv gas constant for water vapor (j/kg/k)
48 !-- dt time step (s)
49 !-- dx model grid interval (m)
50 !-- dz height of model layers (m)
51 !-- xland land mask (1 for land, 2 for water)
52 !-- br bulk richardson number in surface layer
53 !-- pblh planetary boundary layer height (m)
54 !-- ep1 constant for virtual temperature (r_v/r_d - 1) (dimensionless)
55 !-- ids start index for i in domain
56 !-- ide end index for i in domain
57 !-- jds start index for j in domain
58 !-- jde end index for j in domain
59 !-- kds start index for k in domain
60 !-- kde end index for k in domain
61 !-- ims start index for i in memory
62 !-- ime end index for i in memory
63 !-- jms start index for j in memory
64 !-- jme end index for j in memory
65 !-- kms start index for k in memory
66 !-- kme end index for k in memory
67 !-- its start index for i in tile
68 !-- ite end index for i in tile
69 !-- jts start index for j in tile
70 !-- jte end index for j in tile
71 !-- kts start index for k in tile
72 !-- kte end index for k in tile
73 !
74 !-------------------------------------------------------------------------------
75 integer, intent(in ) :: ids,ide, jds,jde, kds,kde, &
76 ims,ime, jms,jme, kms,kme, &
77 its,ite, jts,jte, kts,kte
78 integer, intent(in ) :: itimestep,gwd_opt,gwd_diags
79 !
80 real, intent(in ) :: dt,dx,cp,g,rd,rv,ep1,pi
81 !
82 real, dimension( ims:ime, kms:kme, jms:jme ) , &
83 intent(in ) :: qv3d, &
84 p3d, &
85 pi3d, &
86 t3d, &
87 z, &
88 dz
89 real, dimension( ims:ime, kms:kme, jms:jme ) , &
90 intent(in ) :: p3di
91 !
92 real, dimension( ims:ime, kms:kme, jms:jme ) , &
93 intent(inout) :: rublten, &
94 rvblten, &
95 rthblten
96 real, dimension( ims:ime, kms:kme, jms:jme ), optional , &
97 intent(inout) :: dtaux3d_ls,dtauy3d_ls,dtaux3d_bl,dtauy3d_bl, &
98 dtaux3d_ss,dtauy3d_ss,dtaux3d_fd,dtauy3d_fd
99 !
100 real, dimension( ims:ime, kms:kme) :: &
101 dtaux2d_ls,dtauy2d_ls,dtaux2d_bl,dtauy2d_bl, &
102 dtaux2d_ss,dtauy2d_ss,dtaux2d_fd,dtauy2d_fd
104 real, dimension( ims:ime, kms:kme, jms:jme ) , &
105 intent(in ) :: u3d, &
106 v3d
107 !
108 integer, dimension( ims:ime, jms:jme ) , &
109 intent(in ) :: kpbl2d
111 real, dimension( ims:ime, jms:jme ) , &
112 intent(in ) :: pblh, &
113 br, &
114 xland
116 real, dimension( ims:ime, jms:jme ), optional , &
117 intent(inout ) :: dusfcg_ls,dvsfcg_ls,dusfcg_bl,dvsfcg_bl, &
118 dusfcg_ss,dvsfcg_ss,dusfcg_fd,dvsfcg_fd
120 real, dimension( ims:ime ) :: dusfc_ls,dvsfc_ls,dusfc_bl,dvsfc_bl, &
121 dusfc_ss,dvsfc_ss,dusfc_fd,dvsfc_fd
122 !
123 real, dimension( ims:ime, jms:jme ) , &
124 intent(in ) :: var2d, &
125 oc12d, &
126 oa2d1,oa2d2,oa2d3,oa2d4, &
127 ol2d1,ol2d2,ol2d3,ol2d4, &
128 ! additional topographic statistics based on 1km orographic dataset
129 ! for small-scale orographic drag schemes
130 var2dss,oc12dss, &
131 oa2d1ss,oa2d2ss,oa2d3ss,oa2d4ss, &
132 ol2d1ss,ol2d2ss,ol2d3ss,ol2d4ss, &
133 ! info for converting velocities to/from zonal and meridional components
134 sina,cosa
135 !
136 real, dimension( kms:kme ) , &
137 optional , &
138 intent(in ) :: znu, &
139 znw
140 !
141 real, optional, intent(in ) :: p_top
143 ! Stochastic fields
144 INTEGER, INTENT(IN) ::spp_pbl
145 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN),OPTIONAL ::pattern_spp_pbl
146 REAL, DIMENSION(ITS:ITE) :: rstoch1D
148 !
149 !local
150 !
151 real, dimension( its:ite, kts:kte ) :: delprsi, &
152 pdh
153 real, dimension( its:ite, kts:kte+1 ) :: pdhi
154 real, dimension( its:ite, 4 ) :: oa4, &
155 ol4, &
156 oa4ss, &
157 ol4ss
158 real, dimension( its:ite, kts:kte ) :: ugeo, vgeo, dudt, dvdt
160 integer :: i,j,k,kdt,kpblmax
161 !
162 do k = kts,kte
163 if(znu(k).gt.0.6) kpblmax = k + 1
164 enddo
165 !
166 do j = jts,jte
168 do k = kts,kte+1
169 do i = its,ite
170 if(k.le.kte)pdh(i,k) = p3d(i,k,j)
171 pdhi(i,k) = p3di(i,k,j)
172 enddo
173 enddo
174 !
175 do k = kts,kte
176 do i = its,ite
177 delprsi(i,k) = pdhi(i,k)-pdhi(i,k+1)
178 ! rotate winds to zonal/meridional
179 ugeo(i,k) = u3d(i,k,j)*cosa(i,j) - v3d(i,k,j)*sina(i,j)
180 vgeo(i,k) = u3d(i,k,j)*sina(i,j) + v3d(i,k,j)*cosa(i,j)
181 dudt(i,k) = 0.0
182 dvdt(i,k) = 0.0
183 enddo
184 enddo
186 !Add SPP
187 if (spp_pbl==1) then
188 do i = its,ite
189 rstoch1D(i)=pattern_spp_pbl(i,kts,j)
190 enddo
191 else
192 do i = its,ite
193 rstoch1D(i)=0.0
194 enddo
195 endif
196 !end SPP
197 do i = its,ite
198 oa4(i,1) = oa2d1(i,j)
199 oa4(i,2) = oa2d2(i,j)
200 oa4(i,3) = oa2d3(i,j)
201 oa4(i,4) = oa2d4(i,j)
202 ol4(i,1) = ol2d1(i,j)
203 ol4(i,2) = ol2d2(i,j)
204 ol4(i,3) = ol2d3(i,j)
205 ol4(i,4) = ol2d4(i,j)
206 oa4ss(i,1) = oa2d1ss(i,j)
207 oa4ss(i,2) = oa2d2ss(i,j)
208 oa4ss(i,3) = oa2d3ss(i,j)
209 oa4ss(i,4) = oa2d4ss(i,j)
210 ol4ss(i,1) = ol2d1ss(i,j)
211 ol4ss(i,2) = ol2d2ss(i,j)
212 ol4ss(i,3) = ol2d3ss(i,j)
213 ol4ss(i,4) = ol2d4ss(i,j)
214 enddo
215 call gwdo2d(dudt=dudt(its,kts),dvdt=dvdt(its,kts) &
216 ,dthdt=rthblten(ims,kms,j) &
217 ,dtaux2d_ls=dtaux2d_ls,dtauy2d_ls=dtauy2d_ls &
218 ,dtaux2d_bl=dtaux2d_bl,dtauy2d_bl=dtauy2d_bl &
219 ,dtaux2d_ss=dtaux2d_ss,dtauy2d_ss=dtauy2d_ss &
220 ,dtaux2d_fd=dtaux2d_fd,dtauy2d_fd=dtauy2d_fd &
221 ,u1=ugeo(its,kts),v1=vgeo(its,kts) &
222 ,t1=t3d(ims,kms,j),q1=qv3d(ims,kms,j) &
223 ,del=delprsi(its,kts) &
224 ,prsi=pdhi(its,kts) &
225 ,prsl=pdh(its,kts),prslk=pi3d(ims,kms,j) &
226 ,zl=z(ims,kms,j),rcl=1.0 &
227 ,xland1=xland(ims,j),br1=br(ims,j),hpbl=pblh(ims,j) &
228 ,dz2=dz(ims,kms,j) &
229 ,kpblmax=kpblmax &
230 ,dusfc_ls=dusfc_ls,dvsfc_ls=dvsfc_ls &
231 ,dusfc_bl=dusfc_bl,dvsfc_bl=dvsfc_bl &
232 ,dusfc_ss=dusfc_ss,dvsfc_ss=dvsfc_ss &
233 ,dusfc_fd=dusfc_fd,dvsfc_fd=dvsfc_fd &
234 ,var=var2d(ims,j),oc1=oc12d(ims,j) &
235 ,oa4=oa4,ol4=ol4 &
236 ,varss=var2dss(ims,j),oc1ss=oc12dss(ims,j) &
237 ,oa4ss=oa4ss,ol4ss=ol4ss &
238 ,g=g,cp=cp,rd=rd,rv=rv,fv=ep1,pi=pi &
239 ,dxmeter=dx,deltim=dt &
240 ,kpbl=kpbl2d(ims,j),kdt=itimestep,lat=j &
241 ,rstoch=rstoch1d &
242 ,ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde &
243 ,ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme &
244 ,its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte )
246 do k = kts,kte
247 do i = its,ite
248 ! rotate tendencies from zonal/meridional to model grid
249 rublten(i,k,j) = rublten(i,k,j)+dudt(i,k)*cosa(i,j) + dvdt(i,k)*sina(i,j)
250 rvblten(i,k,j) = rvblten(i,k,j)-dudt(i,k)*sina(i,j) + dvdt(i,k)*cosa(i,j)
251 end do
252 end do
255 IF (gwd_diags == 1) then !research mode
256 do k = kts,kte
257 do i = its,ite
258 ! rotate tendencies from zonal/meridional to model grid
259 dtaux3d_ls(i,k,j)=dtaux2d_ls(i,k)*cosa(i,j)+dtauy2d_ls(i,k)*sina(i,j)
260 dtaux3d_bl(i,k,j)=dtaux2d_bl(i,k)*cosa(i,j)+dtauy2d_bl(i,k)*sina(i,j)
261 dtaux3d_ss(i,k,j)=dtaux2d_ss(i,k)*cosa(i,j)+dtauy2d_ss(i,k)*sina(i,j)
262 dtaux3d_fd(i,k,j)=dtaux2d_fd(i,k)*cosa(i,j)+dtauy2d_fd(i,k)*sina(i,j)
263 dtauy3d_ls(i,k,j)=-dtaux2d_ls(i,k)*sina(i,j)+dtauy2d_ls(i,k)*cosa(i,j)
264 dtauy3d_bl(i,k,j)=-dtaux2d_bl(i,k)*sina(i,j)+dtauy2d_bl(i,k)*cosa(i,j)
265 dtauy3d_ss(i,k,j)=-dtaux2d_ss(i,k)*sina(i,j)+dtauy2d_ss(i,k)*cosa(i,j)
266 dtauy3d_fd(i,k,j)=-dtaux2d_fd(i,k)*sina(i,j)+dtauy2d_fd(i,k)*cosa(i,j)
267 enddo
268 enddo
269 do i = its,ite
270 ! rotate tendencies from zonal/meridional to model grid
271 dusfcg_ls(i,j)=dusfc_ls(i)*cosa(i,j)+dvsfc_ls(i)*sina(i,j)
272 dusfcg_bl(i,j)=dusfc_bl(i)*cosa(i,j)+dvsfc_bl(i)*sina(i,j)
273 dusfcg_ss(i,j)=dusfc_ss(i)*cosa(i,j)+dvsfc_ss(i)*sina(i,j)
274 dusfcg_fd(i,j)=dusfc_fd(i)*cosa(i,j)+dvsfc_fd(i)*sina(i,j)
275 dvsfcg_ls(i,j)=-dusfc_ls(i)*sina(i,j)+dvsfc_ls(i)*cosa(i,j)
276 dvsfcg_bl(i,j)=-dusfc_bl(i)*sina(i,j)+dvsfc_bl(i)*cosa(i,j)
277 dvsfcg_ss(i,j)=-dusfc_ss(i)*sina(i,j)+dvsfc_ss(i)*cosa(i,j)
278 dvsfcg_fd(i,j)=-dusfc_fd(i)*sina(i,j)+dvsfc_fd(i)*cosa(i,j)
279 enddo
280 ENDIF
282 enddo !end-j
283 !
284 end subroutine gwdo_gsl
285 !-------------------------------------------------------------------------------
286 !
287 !-------------------------------------------------------------------------------
288 subroutine gwdo2d(dudt,dvdt,dthdt,dtaux2d_ls,dtauy2d_ls, &
289 dtaux2d_bl,dtauy2d_bl,dtaux2d_ss,dtauy2d_ss, &
290 dtaux2d_fd,dtauy2d_fd,u1,v1,t1,q1, &
291 del, &
292 prsi,prsl,prslk,zl,rcl, &
293 xland1,br1,hpbl,dz2, &
294 kpblmax,dusfc_ls,dvsfc_ls,dusfc_bl,dvsfc_bl, &
295 dusfc_ss,dvsfc_ss,dusfc_fd,dvsfc_fd,var,oc1,oa4,ol4, &
296 varss,oc1ss,oa4ss,ol4ss, &
297 g,cp,rd,rv,fv,pi,dxmeter,deltim,kpbl,kdt,lat,rstoch, &
298 ids,ide, jds,jde, kds,kde, &
299 ims,ime, jms,jme, kms,kme, &
300 its,ite, jts,jte, kts,kte)
301 !-------------------------------------------------------------------------------
302 !
303 ! this code handles the time tendencies of u v due to the effect of mountain
304 ! induced gravity wave drag from sub-grid scale orography. this routine
305 ! not only treats the traditional upper-level wave breaking due to mountain
306 ! variance (alpert 1988), but also the enhanced lower-tropospheric wave
307 ! breaking due to mountain convexity and asymmetry (kim and arakawa 1995).
308 ! thus, in addition to the terrain height data in a model grid gox,
309 ! additional 10-2d topographic statistics files are needed, including
310 ! orographic standard deviation (var), convexity (oc1), asymmetry (oa4)
311 ! and ol (ol4). these data sets are prepared based on the 30 sec usgs orography
312 ! hong (1999). the current scheme was implmented as in hong et al.(2008)
313 !
314 ! coded by song-you hong and young-joon kim and implemented by song-you hong
315 !
316 ! program history log:
317 ! 2014-10-01 Hyun-Joo Choi (from KIAPS) flow-blocking drag of kim and doyle
318 ! with blocked height by dividing streamline theory
319 ! 2017-04-06 Joseph Olson (from Gert-Jan Steeneveld) added small-scale
320 ! orographic grabity wave drag:
321 ! 2017-09-15 Joseph Olson, with some bug fixes from Michael Toy: added the
322 ! topographic form drag of Beljaars et al. (2004, QJRMS)
323 ! Activation of each component is done by specifying the integer-parameters
324 ! (defined below) to 0: inactive or 1: active
325 ! gsl_gwd_ls = 0 or 1: large-scale
326 ! gsl_gwd_bl = 0 or 1: blocking drag
327 ! gsl_gwd_ss = 0 or 1: small-scale gravity wave drag
328 ! gsl_gwd_fd = 0 or 1: topographic form drag
329 ! 2017-09-25 Michael Toy (from NCEP GFS model) added dissipation heating
330 ! gsl_diss_ht_opt = 0: dissipation heating off
331 ! gsl_diss_ht_opt = 1: dissipation heating on
332 !
333 ! references:
334 ! hong et al. (2008), wea. and forecasting
335 ! kim and doyle (2005), Q. J. R. Meteor. Soc.
336 ! kim and arakawa (1995), j. atmos. sci.
337 ! alpet et al. (1988), NWP conference.
338 ! hong (1999), NCEP office note 424.
339 ! steeneveld et al (2008), JAMC
340 ! Tsiringakis et al. (2017), Q. J. R. Meteor. Soc.
341 !
342 ! notice : comparible or lower resolution orography files than model resolution
343 ! are desirable in preprocess (wps) to prevent weakening of the drag
344 !-------------------------------------------------------------------------------
345 !
346 ! input
347 ! dudt (its:ite,kts:kte) non-lin tendency for u wind component
348 ! dvdt (its:ite,kts:kte) non-lin tendency for v wind component
349 ! u1(its:ite,kts:kte) zonal wind / sqrt(rcl) m/sec at t0-dt
350 ! v1(its:ite,kts:kte) meridional wind / sqrt(rcl) m/sec at t0-dt
351 ! t1(ims:ime,kms:kme) temperature deg k at t0-dt
352 ! q1(ims:ime,kms:kme) specific humidity at t0-dt
353 !
354 ! rcl a scaling factor = reciprocal of square of cos(lat)
355 ! for gmp. rcl=1 if u1 and v1 are wind components.
356 ! deltim time step secs
357 ! del(kts:kte) positive increment of pressure across layer (pa)
358 !
359 ! output
360 ! dudt, dvdt wind tendency due to gwdo
361 !
362 !-------------------------------------------------------------------------------
363 implicit none
364 !-------------------------------------------------------------------------------
365 integer :: kdt,lat,latd,lond,kpblmax, &
366 ids,ide, jds,jde, kds,kde, &
367 ims,ime, jms,jme, kms,kme, &
368 its,ite, jts,jte, kts,kte
369 !
370 real :: g,rd,rv,fv,cp,pi,dxmeter,deltim,rcl
371 real :: dudt(its:ite,kts:kte),dvdt(its:ite,kts:kte), &
372 dthdt(ims:ime,kms:kme), &
373 dtaux2d_ls(ims:ime,kms:kme),dtauy2d_ls(ims:ime,kms:kme), &
374 dtaux2d_bl(ims:ime,kms:kme),dtauy2d_bl(ims:ime,kms:kme), &
375 dtaux2d_ss(ims:ime,kms:kme),dtauy2d_ss(ims:ime,kms:kme), &
376 dtaux2d_fd(ims:ime,kms:kme),dtauy2d_fd(ims:ime,kms:kme), &
377 t1(ims:ime,kms:kme),q1(ims:ime,kms:kme), &
378 zl(ims:ime,kms:kme),prsl(its:ite,kts:kte), &
379 prslk(ims:ime,kms:kme)
380 real :: u1(its:ite,kts:kte),v1(its:ite,kts:kte)
381 real :: prsi(its:ite,kts:kte+1),del(its:ite,kts:kte)
382 real :: oa4(its:ite,4),ol4(its:ite,4), &
383 oa4ss(its:ite,4),ol4ss(its:ite,4)
384 !
385 ! GSL surface drag options to regulate specific components
386 ! Each component is tapered off automatically as a function of dx, so best to
387 ! keep them activated (=1).
388 integer, parameter :: &
389 gsl_gwd_ls = 1, & ! large-scale gravity wave drag
390 gsl_gwd_bl = 1, & ! blocking drag
391 gsl_gwd_ss = 1, & ! small-scale gravity wave drag (Steeneveld et al. 2008)
392 gsl_gwd_fd = 1, & ! form drag (Beljaars et al. 2004, QJRMS)
393 gsl_diss_ht_opt = 0
394 !
395 ! added for small-scale orographic wave drag
396 real, dimension(its:ite,kts:kte) :: utendwave,vtendwave,thx,thvx,za
397 real, dimension(ims:ime), intent(in) :: br1,hpbl,xland1
398 real, dimension(its:ite) :: govrth
399 real, dimension(ims:ime,kms:kme), intent(in) :: dz2
400 real, dimension(its:ite,kts:kte+1) :: zq
401 real :: tauwavex0,tauwavey0,XNBV,density, &
402 tvcon,hpbl2
403 integer :: kpbl2,kvar
404 ! added for SPP
405 real, dimension(its:ite) :: rstoch
406 real :: var_stoch(ims:ime), &
407 varss_stoch(ims:ime)
408 !
409 integer :: kpbl(ims:ime)
410 real :: var(ims:ime),oc1(ims:ime), &
411 varss(ims:ime),oc1ss(ims:ime), &
412 dusfc_ls(ims:ime),dvsfc_ls(ims:ime), &
413 dusfc_bl(ims:ime),dvsfc_bl(ims:ime), &
414 dusfc_ss(ims:ime),dvsfc_ss(ims:ime), &
415 dusfc_fd(ims:ime),dvsfc_fd(ims:ime)
416 ! Variables for scale-awareness:
417 ! Small-scale GWD + turbulent form drag
418 real, parameter :: dxmin_ss = 1000., dxmax_ss = 12000. ! min,max range of tapering (m)
419 ! Large-scale GWD
420 real, parameter :: dxmin_ls = 3000., dxmax_ls = 13000. ! min,max range of tapering (m)
421 real :: ss_taper, ls_taper ! small- and large-scale tapering factors (-)
424 !
425 ! added Beljaars orographic form drag
426 real, dimension(its:ite,kts:kte) :: utendform,vtendform
427 real :: a1,a2,wsp
428 real :: H_efold
429 ! critical richardson number for wave breaking : ! larger drag with larger value
430 !
431 real,parameter :: ric = 0.25
432 !
433 real,parameter :: dw2min = 1.
434 real,parameter :: rimin = -100.
435 real,parameter :: bnv2min = 1.0e-5
436 real,parameter :: efmin = 0.0
437 real,parameter :: efmax = 10.0
438 real,parameter :: xl = 4.0e4
439 real,parameter :: critac = 1.0e-5
440 real,parameter :: gmax = 1.
441 real,parameter :: veleps = 1.0
442 real,parameter :: factop = 0.5
443 real,parameter :: frc = 1.0
444 real,parameter :: ce = 0.8
445 real,parameter :: cg = 0.5
446 integer,parameter :: kpblmin = 2
448 !
449 ! Variables for limiting topographic standard deviation (var)
450 real, parameter :: varmax_ss = 35., &
451 varmax_fd = 160., &
452 beta_ss = 0.1, &
453 beta_fd = 0.2
454 real :: var_temp
455 real, dimension(its:ite) :: &
456 varmax_ss_stoch, &
457 varmax_fd_stoch
458 !
459 ! local variables
460 !
461 integer :: i,k,lcap,lcapp1,nwd,idir, &
462 klcap,kp1,ikount,kk
463 !
464 real :: rcs,rclcs,csg,fdir,cleff,cleff_ss,cs,rcsks, &
465 wdir,ti,rdz,temp,tem2,dw2,shr2,bvf2,rdelks, &
466 wtkbj,tem,gfobnv,hd,fro,rim,temc,tem1,efact, &
467 temv,dtaux,dtauy,eng0,eng1
468 !
469 logical :: ldrag(its:ite),icrilv(its:ite), &
470 flag(its:ite),kloop1(its:ite)
471 !
472 real :: taub(its:ite),taup(its:ite,kts:kte+1), &
473 xn(its:ite),yn(its:ite), &
474 ubar(its:ite),vbar(its:ite), &
475 fr(its:ite),ulow(its:ite), &
476 rulow(its:ite),bnv(its:ite), &
477 oa(its:ite),ol(its:ite), &
478 oass(its:ite),olss(its:ite), &
479 roll(its:ite),dtfac(its:ite), &
480 brvf(its:ite),xlinv(its:ite), &
481 delks(its:ite),delks1(its:ite), &
482 bnv2(its:ite,kts:kte),usqj(its:ite,kts:kte), &
483 taud_ls(its:ite,kts:kte),taud_bl(its:ite,kts:kte), &
484 ro(its:ite,kts:kte), &
485 vtk(its:ite,kts:kte),vtj(its:ite,kts:kte), &
486 zlowtop(its:ite),velco(its:ite,kts:kte-1), &
487 coefm(its:ite),coefm_ss(its:ite)
488 !
489 integer :: kbl(its:ite),klowtop(its:ite)
490 !
491 logical :: iope
492 integer,parameter :: mdir=8
493 integer :: nwdir(mdir)
494 data nwdir/6,7,5,8,2,3,1,4/
495 !
496 ! variables for flow-blocking drag
497 !
498 real,parameter :: frmax = 10.
499 real,parameter :: olmin = 1.0e-5
500 real,parameter :: odmin = 0.1
501 real,parameter :: odmax = 10.
502 integer :: komax(its:ite)
503 integer :: kblk
504 real :: cd
505 real :: zblk,tautem
506 real :: pe,ke
507 real :: delx,dely,dxy4(4),dxy4p(4)
508 real :: dxy(its:ite),dxyp(its:ite)
509 real :: ol4p(4),olp(its:ite),od(its:ite)
510 real :: taufb(its:ite,kts:kte+1)
512 !
513 ! misc. coefficients
514 !
515 real, parameter :: &
516 clf_coeff = 3.4E+07, & ! "effective length" coeff or large-scale GWD
517 ! Note: Can be used as "tuning knob"
518 clf_coeff_ss = 0.1, & ! Adjustment for small-scale GWD
519 a1_coeff = 0.00026615161, & ! Coefficient for TOFD from Beljaars et al. (2004)
520 a2_coeff = 0.005363, & ! ""
521 TOFD_coeff = 0.0759, & ! ""
522 Hefold_nom = 1500. ! Nominal TOFD e-folding height (m)
525 !
526 !---- constants
527 !
528 rcs = sqrt(rcl)
529 cs = 1. / sqrt(rcl)
530 csg = cs * g
531 lcap = kte
532 lcapp1 = lcap + 1
533 fdir = mdir / (2.0*pi)
534 !
535 !--- calculate scale-aware tapering factors
536 !
537 if ( dxmeter .ge. dxmax_ls ) then
538 ls_taper = 1.
539 else
540 if ( dxmeter .le. dxmin_ls) then
541 ls_taper = 0.
542 else
543 ls_taper = 0.5 * ( SIN(pi*(dxmeter-0.5*(dxmax_ls+dxmin_ls))/ &
544 (dxmax_ls-dxmin_ls)) + 1. )
545 end if
546 end if
547 if ( dxmeter .ge. dxmax_ss ) then
548 ss_taper = 1.
549 else
550 if ( dxmeter .le. dxmin_ss) then
551 ss_taper = 0.
552 else
553 ss_taper = dxmax_ss * (1. - dxmin_ss/dxmeter)/(dxmax_ss-dxmin_ss)
554 end if
555 end if
556 !
557 !--- calculate length of grid for flow-blocking drag
558 !
559 delx = dxmeter
560 dely = dxmeter
561 dxy4(1) = delx
562 dxy4(2) = dely
563 dxy4(3) = sqrt(delx*delx + dely*dely)
564 dxy4(4) = dxy4(3)
565 dxy4p(1) = dxy4(2)
566 dxy4p(2) = dxy4(1)
567 dxy4p(3) = dxy4(4)
568 dxy4p(4) = dxy4(3)
569 !
570 !
571 !-----initialize arrays
572 !
573 dtaux = 0.0
574 dtauy = 0.0
575 do i = its,ite
576 klowtop(i) = 0
577 kbl(i) = 0
578 enddo
579 !
580 do i = its,ite
581 xn(i) = 0.0
582 yn(i) = 0.0
583 ubar (i) = 0.0
584 vbar (i) = 0.0
585 roll (i) = 0.0
586 taub (i) = 0.0
587 oa(i) = 0.0
588 ol(i) = 0.0
589 oass(i) = 0.0
590 olss(i) = 0.0
591 ulow (i) = 0.0
592 dtfac(i) = 1.0
593 ldrag(i) = .false.
594 icrilv(i) = .false.
595 flag(i) = .true.
596 enddo
597 !
598 do k = kts,kte
599 do i = its,ite
600 usqj(i,k) = 0.0
601 bnv2(i,k) = 0.0
602 vtj(i,k) = 0.0
603 vtk(i,k) = 0.0
604 taup(i,k) = 0.0
605 taud_ls(i,k) = 0.0
606 taud_bl(i,k) = 0.0
607 dtaux2d_ls(i,k)= 0.0
608 dtauy2d_ls(i,k)= 0.0
609 dtaux2d_bl(i,k)= 0.0
610 dtauy2d_bl(i,k)= 0.0
611 dtaux2d_ss(i,k)= 0.0
612 dtauy2d_ss(i,k)= 0.0
613 dtaux2d_fd(i,k)= 0.0
614 dtauy2d_fd(i,k)= 0.0
615 enddo
616 enddo
617 !
618 do i = its,ite
619 dusfc_ls(i) = 0.0
620 dvsfc_ls(i) = 0.0
621 dusfc_bl(i) = 0.0
622 dvsfc_bl(i) = 0.0
623 dusfc_ss(i) = 0.0
624 dvsfc_ss(i) = 0.0
625 dusfc_fd(i) = 0.0
626 dvsfc_fd(i) = 0.0
627 enddo
628 !
629 do i = its,ite
630 taup(i,kte+1) = 0.0
631 xlinv(i) = 1.0/xl
632 enddo
633 !
634 ! SPP (rstoch will be 0 is spp_pbl = 0)
635 do i = its,ite
636 var_stoch(i) = var(i) + var(i)*0.666*rstoch(i)
637 varss_stoch(i) = varss(i) + varss(i)*0.666*rstoch(i)
638 varmax_ss_stoch(i) = varmax_ss + varmax_ss*0.5*rstoch(i)
639 varmax_fd_stoch(i) = varmax_fd + varmax_fd*0.5*rstoch(i)
640 enddo
641 !
642 ! initialize array for flow-blocking drag
643 !
644 taufb(its:ite,kts:kte+1) = 0.0
645 komax(its:ite) = 0
646 !
647 do k = kts,kte
648 do i = its,ite
649 vtj(i,k) = t1(i,k) * (1.+fv*q1(i,k))
650 vtk(i,k) = vtj(i,k) / prslk(i,k)
651 ro(i,k) = 1./rd * prsl(i,k) / vtj(i,k) ! density kg/m**3
652 enddo
653 enddo
654 !
655 ! determine reference level: maximum of 2*var and pbl heights
656 !
657 do i = its,ite
658 zlowtop(i) = 2. * var_stoch(i)
659 enddo
660 !
661 do i = its,ite
662 kloop1(i) = .true.
663 enddo
664 !
665 do k = kts+1,kte
666 do i = its,ite
667 if(kloop1(i).and.zl(i,k)-zl(i,1).ge.zlowtop(i)) then
668 klowtop(i) = k+1
669 kloop1(i) = .false.
670 endif
671 enddo
672 enddo
673 !
674 do i = its,ite
675 kbl(i) = max(kpbl(i), klowtop(i))
676 kbl(i) = max(min(kbl(i),kpblmax),kpblmin)
677 enddo
678 !
679 ! determine the level of maximum orographic height
680 !
681 ! komax(:) = kbl(:)
682 komax(:) = klowtop(:) - 1 ! modification by NOAA/GSL March 2018
683 !
684 do i = its,ite
685 delks(i) = 1.0 / (prsi(i,1) - prsi(i,kbl(i)))
686 delks1(i) = 1.0 / (prsl(i,1) - prsl(i,kbl(i)))
687 enddo
688 !
689 ! compute low level averages within pbl
690 !
691 do k = kts,kpblmax
692 do i = its,ite
693 if (k.lt.kbl(i)) then
694 rcsks = rcs * del(i,k) * delks(i)
695 rdelks = del(i,k) * delks(i)
696 ubar(i) = ubar(i) + rcsks * u1(i,k) ! pbl u mean
697 vbar(i) = vbar(i) + rcsks * v1(i,k) ! pbl v mean
698 roll(i) = roll(i) + rdelks * ro(i,k) ! ro mean
699 endif
700 enddo
701 enddo
702 !
703 ! figure out low-level horizontal wind direction
704 !
705 ! nwd 1 2 3 4 5 6 7 8
706 ! wd w s sw nw e n ne se
707 !
708 do i = its,ite
709 wdir = atan2(ubar(i),vbar(i)) + pi
710 idir = mod(nint(fdir*wdir),mdir) + 1
711 nwd = nwdir(idir)
712 oa(i) = (1-2*int( (nwd-1)/4 )) * oa4(i,mod(nwd-1,4)+1)
713 ol(i) = ol4(i,mod(nwd-1,4)+1)
714 ! Repeat for small-scale gwd
715 oass(i) = (1-2*int( (nwd-1)/4 )) * oa4ss(i,mod(nwd-1,4)+1)
716 olss(i) = ol4ss(i,mod(nwd-1,4)+1)
718 !
719 !----- compute orographic width along (ol) and perpendicular (olp)
720 !----- the direction of wind
721 !
722 ol4p(1) = ol4(i,2)
723 ol4p(2) = ol4(i,1)
724 ol4p(3) = ol4(i,4)
725 ol4p(4) = ol4(i,3)
726 olp(i) = ol4p(mod(nwd-1,4)+1)
727 !
728 !----- compute orographic direction (horizontal orographic aspect ratio)
729 !
730 od(i) = olp(i)/max(ol(i),olmin)
731 od(i) = min(od(i),odmax)
732 od(i) = max(od(i),odmin)
733 !
734 !----- compute length of grid in the along(dxy) and cross(dxyp) wind directions
735 !
736 dxy(i) = dxy4(MOD(nwd-1,4)+1)
737 dxyp(i) = dxy4p(MOD(nwd-1,4)+1)
738 enddo
739 !
740 ! END INITIALIZATION; BEGIN GWD CALCULATIONS:
741 !
742 IF ( ((gsl_gwd_ls .EQ. 1).or.(gsl_gwd_bl .EQ. 1)).and. &
743 (ls_taper .GT. 1.E-02) ) THEN !====
744 !
745 !--- saving richardson number in usqj for migwdi
746 !
747 do k = kts,kte-1
748 do i = its,ite
749 ti = 2.0 / (t1(i,k)+t1(i,k+1))
750 rdz = 1./(zl(i,k+1) - zl(i,k))
751 tem1 = u1(i,k) - u1(i,k+1)
752 tem2 = v1(i,k) - v1(i,k+1)
753 dw2 = rcl*(tem1*tem1 + tem2*tem2)
754 shr2 = max(dw2,dw2min) * rdz * rdz
755 bvf2 = g*(g/cp+rdz*(vtj(i,k+1)-vtj(i,k))) * ti
756 usqj(i,k) = max(bvf2/shr2,rimin)
757 bnv2(i,k) = 2.0*g*rdz*(vtk(i,k+1)-vtk(i,k))/(vtk(i,k+1)+vtk(i,k))
758 bnv2(i,k) = max( bnv2(i,k), bnv2min )
759 enddo
760 enddo
761 !
762 !----compute the "low level" or 1/3 wind magnitude (m/s)
763 !
764 do i = its,ite
765 ulow(i) = max(sqrt(ubar(i)*ubar(i) + vbar(i)*vbar(i)), 1.0)
766 rulow(i) = 1./ulow(i)
767 enddo
768 !
769 do k = kts,kte-1
770 do i = its,ite
771 velco(i,k) = (0.5*rcs) * ((u1(i,k)+u1(i,k+1)) * ubar(i) &
772 + (v1(i,k)+v1(i,k+1)) * vbar(i))
773 velco(i,k) = velco(i,k) * rulow(i)
774 if ((velco(i,k).lt.veleps) .and. (velco(i,k).gt.0.)) then
775 velco(i,k) = veleps
776 endif
777 enddo
778 enddo
779 !
780 ! no drag when critical level in the base layer
781 !
782 do i = its,ite
783 ldrag(i) = velco(i,1).le.0.
784 enddo
785 !
786 ! no drag when velco.lt.0
787 !
788 do k = kpblmin,kpblmax
789 do i = its,ite
790 if (k .lt. kbl(i)) ldrag(i) = ldrag(i).or. velco(i,k).le.0.
791 enddo
792 enddo
793 !
794 ! no drag when bnv2.lt.0
795 !
796 do k = kts,kpblmax
797 do i = its,ite
798 if (k .lt. kbl(i)) ldrag(i) = ldrag(i).or. bnv2(i,k).lt.0.
799 enddo
800 enddo
801 !
802 !-----the low level weighted average ri is stored in usqj(1,1; im)
803 !-----the low level weighted average n**2 is stored in bnv2(1,1; im)
804 !---- this is called bnvl2 in phys_gwd_alpert_sub not bnv2
805 !---- rdelks (del(k)/delks) vert ave factor so we can * instead of /
806 !
807 do i = its,ite
808 wtkbj = (prsl(i,1)-prsl(i,2)) * delks1(i)
809 bnv2(i,1) = wtkbj * bnv2(i,1)
810 usqj(i,1) = wtkbj * usqj(i,1)
811 enddo
812 !
813 do k = kpblmin,kpblmax
814 do i = its,ite
815 if (k .lt. kbl(i)) then
816 rdelks = (prsl(i,k)-prsl(i,k+1)) * delks1(i)
817 bnv2(i,1) = bnv2(i,1) + bnv2(i,k) * rdelks
818 usqj(i,1) = usqj(i,1) + usqj(i,k) * rdelks
819 endif
820 enddo
821 enddo
822 !
823 do i = its,ite
824 ldrag(i) = ldrag(i) .or. bnv2(i,1).le.0.0
825 ldrag(i) = ldrag(i) .or. ulow(i).eq.1.0
826 ldrag(i) = ldrag(i) .or. var_stoch(i) .le. 0.0
827 enddo
828 !
829 ! set all ri low level values to the low level value
830 !
831 do k = kpblmin,kpblmax
832 do i = its,ite
833 if (k .lt. kbl(i)) usqj(i,k) = usqj(i,1)
834 enddo
835 enddo
836 !
837 do i = its,ite
838 if (.not.ldrag(i)) then
839 bnv(i) = sqrt( bnv2(i,1) )
840 fr(i) = bnv(i) * rulow(i) * 2. * var_stoch(i) * od(i)
841 fr(i) = min(fr(i),frmax)
842 xn(i) = ubar(i) * rulow(i)
843 yn(i) = vbar(i) * rulow(i)
844 endif
845 enddo
846 !
847 ! compute the base level stress and store it in taub
848 ! calculate enhancement factor, number of mountains & aspect
849 ! ratio const. use simplified relationship between standard
850 ! deviation & critical hgt
851 !
852 do i = its,ite
853 if (.not. ldrag(i)) then
854 efact = (oa(i) + 2.) ** (ce*fr(i)/frc)
855 efact = min( max(efact,efmin), efmax )
856 !!!!!!! cleff ("effective grid length") is highly tunable parameter
857 ! Based on tuning by NOAA on the GFS, the value of cleff should
858 ! be inversely proportional to the square root of the grid length
859 ! (dxmeter). The larger (smaller) the grid size the smaller (larger)
860 ! the value of cleff.
861 cleff = clf_coeff / sqrt(dxmeter)
862 coefm(i) = (1. + ol(i)) ** (oa(i)+1.)
863 xlinv(i) = coefm(i) / cleff
864 tem = fr(i) * fr(i) * oc1(i)
865 gfobnv = gmax * tem / ((tem + cg)*bnv(i))
866 if ( gsl_gwd_ls .NE. 0 ) then
867 taub(i) = xlinv(i) * roll(i) * ulow(i) * ulow(i) &
868 * ulow(i) * gfobnv * efact
869 else ! We've gotten what we need for the blocking scheme
870 taub(i) = 0.0
871 end if
872 else
873 taub(i) = 0.0
874 xn(i) = 0.0
875 yn(i) = 0.0
876 endif
877 enddo
879 ENDIF ! (gsl_gwd_ls .EQ. 1).or.(gsl_gwd_bl .EQ. 1)
880 !=========================================================
881 ! add small-scale wavedrag for stable boundary layer
882 !=========================================================
883 XNBV=0.
884 tauwavex0=0.
885 tauwavey0=0.
886 density=1.2
887 utendwave=0.
888 vtendwave=0.
889 zq=0.
890 !
891 IF ( (gsl_gwd_ss .EQ. 1).and.(ss_taper.GT.1.E-02) ) THEN
892 !
893 ! declaring potential temperature
894 !
895 do k = kts,kte
896 do i = its,ite
897 thx(i,k) = t1(i,k)/prslk(i,k)
898 enddo
899 enddo
900 !
901 do k = kts,kte
902 do i = its,ite
903 tvcon = (1.+fv*q1(i,k))
904 thvx(i,k) = thx(i,k)*tvcon
905 enddo
906 enddo
907 !
908 ! Defining layer height
909 !
910 do k = kts,kte
911 do i = its,ite
912 zq(i,k+1) = dz2(i,k)+zq(i,k)
913 enddo
914 enddo
915 !
916 do k = kts,kte
917 do i = its,ite
918 za(i,k) = 0.5*(zq(i,k)+zq(i,k+1))
919 enddo
920 enddo
922 do i=its,ite
923 hpbl2 = hpbl(i)+10.
924 kpbl2 = kpbl(i)
925 kvar = 1
926 do k=kts+1,MAX(kpbl(i),kts+1)
927 IF (za(i,k)>300.) then
928 kpbl2 = k
929 IF (k == kpbl(i)) then
930 hpbl2 = hpbl(i)+10.
931 ELSE
932 hpbl2 = za(i,k)+10.
933 ENDIF
934 exit
935 ENDIF
936 enddo
937 if((xland1(i)-1.5).le.0. .and. 2.*varss_stoch(i).le.hpbl(i))then
938 if(br1(i).gt.0. .and. thvx(i,kpbl2)-thvx(i,kts) > 0.)then
939 cleff_ss = sqrt(dxy(i)**2 + dxyp(i)**2)
940 cleff_ss = clf_coeff_ss * max(dxmax_ss,cleff_ss)
941 coefm_ss(i) = (1. + olss(i)) ** (oass(i)+1.)
942 xlinv(i) = coefm_ss(i) / cleff_ss
943 govrth(i)=g/(0.5*(thvx(i,kpbl2)+thvx(i,kts)))
944 XNBV=sqrt(govrth(i)*(thvx(i,kpbl2)-thvx(i,kts))/hpbl2)
945 !
946 if(abs(XNBV/u1(i,kpbl2)).gt.xlinv(i))then
947 var_temp = MIN(varss_stoch(i),varmax_ss_stoch(i)) + &
948 MAX(0.,beta_ss*(varss_stoch(i)-varmax_ss_stoch(i)))
949 tauwavex0=0.5*XNBV*xlinv(i)*(2.*var_temp)**2*ro(i,kvar)*u1(i,kvar)
950 tauwavex0=tauwavex0*ss_taper ! "Scale-awareness"
951 else
952 tauwavex0=0.
953 endif
954 !
955 if(abs(XNBV/v1(i,kpbl2)).gt.xlinv(i))then
956 var_temp = MIN(varss_stoch(i),varmax_ss_stoch(i)) + &
957 MAX(0.,beta_ss*(varss_stoch(i)-varmax_ss_stoch(i)))
958 tauwavey0=0.5*XNBV*xlinv(i)*(2.*var_temp)**2*ro(i,kvar)*v1(i,kvar)
959 tauwavey0=tauwavey0*ss_taper ! "Scale-awareness"
960 else
961 tauwavey0=0.
962 endif
963 !
964 do k=kts,kpbl(i) !MIN(kpbl2+1,kte-1)
965 utendwave(i,k)=-1.*tauwavex0*2.*max((1.-za(i,k)/hpbl2),0.)/hpbl2
966 vtendwave(i,k)=-1.*tauwavey0*2.*max((1.-za(i,k)/hpbl2),0.)/hpbl2
967 enddo
968 endif
969 endif
970 enddo ! end i loop
972 do k = kts,kte
973 do i = its,ite
974 dudt(i,k) = dudt(i,k) + utendwave(i,k)
975 dvdt(i,k) = dvdt(i,k) + vtendwave(i,k)
976 dtaux2d_ss(i,k) = utendwave(i,k)
977 dtauy2d_ss(i,k) = vtendwave(i,k)
978 dusfc_ss(i) = dusfc_ss(i) + utendwave(i,k) * del(i,k)
979 dvsfc_ss(i) = dvsfc_ss(i) + vtendwave(i,k) * del(i,k)
980 enddo
981 enddo
983 ENDIF ! end if gsl_gwd_ss == 1
984 !================================================================
985 ! Topographic Form Drag from Beljaars et al. (2004, QJRMS, equ. 16):
986 !================================================================
987 IF ( (gsl_gwd_fd .EQ. 1).and.(ss_taper.GT.1.E-02) ) THEN
989 utendform=0.
990 vtendform=0.
991 zq=0.
993 IF ( (gsl_gwd_ss .NE. 1).and.(ss_taper.GT.1.E-02) ) THEN
994 ! Defining layer height. This is already done above is small-scale GWD is used
995 do k = kts,kte
996 do i = its,ite
997 zq(i,k+1) = dz2(i,k)+zq(i,k)
998 enddo
999 enddo
1001 do k = kts,kte
1002 do i = its,ite
1003 za(i,k) = 0.5*(zq(i,k)+zq(i,k+1))
1004 enddo
1005 enddo
1006 ENDIF
1008 DO i=its,ite
1009 IF ((xland1(i)-1.5) .le. 0.) then
1010 var_temp = MIN(varss_stoch(i),varmax_fd_stoch(i)) + &
1011 MAX(0.,beta_fd*(varss_stoch(i)-varmax_fd_stoch(i)))
1012 a1=a1_coeff*var_temp**2
1013 a2=a1*a2_coeff
1014 ! Revise e-folding height based on PBL height and topographic std. dev. -- M. Toy 3/12/2018
1015 H_efold = max(2*varss_stoch(i),hpbl(i))
1016 H_efold = min(H_efold,Hefold_nom)
1017 DO k=kts,kte
1018 wsp=SQRT(u1(i,k)**2 + v1(i,k)**2)
1019 ! Eqn. (16) of Beljaars et al. (2004)
1020 utendform(i,k)=-TOFD_coeff*wsp*u1(i,k)* &
1021 EXP(-(za(i,k)/H_efold)**1.5)*a2*za(i,k)**(-1.2)*ss_taper
1022 vtendform(i,k)=-TOFD_coeff*wsp*v1(i,k)* &
1023 EXP(-(za(i,k)/H_efold)**1.5)*a2*za(i,k)**(-1.2)*ss_taper
1024 ENDDO
1025 ENDIF
1026 ENDDO
1028 do k = kts,kte
1029 do i = its,ite
1030 dudt(i,k) = dudt(i,k) + utendform(i,k)
1031 dvdt(i,k) = dvdt(i,k) + vtendform(i,k)
1032 dtaux2d_fd(i,k) = utendform(i,k)
1033 dtauy2d_fd(i,k) = vtendform(i,k)
1034 dusfc_fd(i) = dusfc_fd(i) + utendform(i,k) * del(i,k)
1035 dvsfc_fd(i) = dvsfc_fd(i) + vtendform(i,k) * del(i,k)
1036 enddo
1037 enddo
1039 ENDIF ! end if gsl_gwd_fd == 1
1040 !=======================================================
1041 ! More for the large-scale gwd component
1042 IF ( (gsl_gwd_ls .EQ. 1).and.(ls_taper.GT.1.E-02) ) THEN
1043 !
1044 ! now compute vertical structure of the stress.
1045 !
1046 do k = kts,kpblmax
1047 do i = its,ite
1048 if (k .le. kbl(i)) taup(i,k) = taub(i)
1049 enddo
1050 enddo
1051 !
1052 do k = kpblmin, kte-1 ! vertical level k loop!
1053 kp1 = k + 1
1054 do i = its,ite
1055 !
1056 ! unstablelayer if ri < ric
1057 ! unstable layer if upper air vel comp along surf vel <=0 (crit lay)
1058 ! at (u-c)=0. crit layer exists and bit vector should be set (.le.)
1059 !
1060 if (k .ge. kbl(i)) then
1061 icrilv(i) = icrilv(i) .or. ( usqj(i,k) .lt. ric) &
1062 .or. (velco(i,k) .le. 0.0)
1063 brvf(i) = max(bnv2(i,k),bnv2min) ! brunt-vaisala frequency squared
1064 brvf(i) = sqrt(brvf(i)) ! brunt-vaisala frequency
1065 endif
1066 enddo
1067 !
1068 do i = its,ite
1069 if (k .ge. kbl(i) .and. (.not. ldrag(i))) then
1070 if (.not.icrilv(i) .and. taup(i,k) .gt. 0.0 ) then
1071 temv = 1.0 / velco(i,k)
1072 tem1 = coefm(i)/dxy(i)*(ro(i,kp1)+ro(i,k))*brvf(i)*velco(i,k)*0.5
1073 hd = sqrt(taup(i,k) / tem1)
1074 fro = brvf(i) * hd * temv
1075 !
1076 ! rim is the minimum-richardson number by shutts (1985)
1077 !
1078 tem2 = sqrt(usqj(i,k))
1079 tem = 1. + tem2 * fro
1080 rim = usqj(i,k) * (1.-fro) / (tem * tem)
1081 !
1082 ! check stability to employ the 'saturation hypothesis'
1083 ! of lindzen (1981) except at tropospheric downstream regions
1084 !
1085 if (rim .le. ric) then ! saturation hypothesis!
1086 if ((oa(i) .le. 0.).or.(kp1 .ge. kpblmin )) then
1087 temc = 2.0 + 1.0 / tem2
1088 hd = velco(i,k) * (2.*sqrt(temc)-temc) / brvf(i)
1089 taup(i,kp1) = tem1 * hd * hd
1090 endif
1091 else ! no wavebreaking!
1092 taup(i,kp1) = taup(i,k)
1093 endif
1094 endif
1095 endif
1096 enddo
1097 enddo
1098 !
1099 if(lcap.lt.kte) then
1100 do klcap = lcapp1,kte
1101 do i = its,ite
1102 taup(i,klcap) = prsi(i,klcap) / prsi(i,lcap) * taup(i,lcap)
1103 enddo
1104 enddo
1105 endif
1107 ENDIF !END LARGE-SCALE TAU CALCULATION
1109 !===============================================================
1110 !COMPUTE BLOCKING COMPONENT
1111 !===============================================================
1112 IF ( (gsl_gwd_bl .EQ. 1) .and. (ls_taper .GT. 1.E-02) ) THEN
1114 do i = its,ite
1115 if(.not.ldrag(i)) then
1116 !
1117 !------- determine the height of flow-blocking layer
1118 !
1119 kblk = 0
1120 pe = 0.0
1121 do k = kte, kpblmin, -1
1122 if(kblk.eq.0 .and. k.le.komax(i)) then
1123 pe = pe + bnv2(i,k)*(zl(i,komax(i))-zl(i,k))*del(i,k)/g/ro(i,k)
1124 ke = 0.5*((rcs*u1(i,k))**2.+(rcs*v1(i,k))**2.)
1125 !
1126 !---------- apply flow-blocking drag when pe >= ke
1127 !
1128 if(pe.ge.ke) then
1129 kblk = k
1130 kblk = min(kblk,kbl(i))
1131 zblk = zl(i,kblk)-zl(i,kts)
1132 endif
1133 endif
1134 enddo
1135 if(kblk.ne.0) then
1136 !
1137 !--------- compute flow-blocking stress
1138 !
1139 cd = max(2.0-1.0/od(i),0.0)
1140 taufb(i,kts) = 0.5 * roll(i) * coefm(i) / max(dxmax_ls,dxy(i))**2 * cd * dxyp(i) &
1141 * olp(i) * zblk * ulow(i)**2
1142 tautem = taufb(i,kts)/float(kblk-kts)
1143 do k = kts+1, kblk
1144 taufb(i,k) = taufb(i,k-1) - tautem
1145 enddo
1146 !
1147 !----------sum orographic GW stress and flow-blocking stress
1148 !
1149 ! taup(i,:) = taup(i,:) + taufb(i,:) ! Keep taup and taufb separate for now
1150 endif
1151 endif
1152 enddo
1154 ENDIF ! end blocking drag
1155 !===========================================================
1156 IF ( (gsl_gwd_ls .EQ. 1 .OR. gsl_gwd_bl .EQ. 1) .and. (ls_taper .GT. 1.E-02) ) THEN
1157 !
1158 ! calculate - (g)*d(tau)/d(pressure) and deceleration terms dtaux, dtauy
1159 !
1160 do k = kts,kte
1161 do i = its,ite
1162 taud_ls(i,k) = 1. * (taup(i,k+1) - taup(i,k)) * csg / del(i,k)
1163 taud_bl(i,k) = 1. * (taufb(i,k+1) - taufb(i,k)) * csg / del(i,k)
1164 enddo
1165 enddo
1166 !
1167 ! limit de-acceleration (momentum deposition ) at top to 1/2 value
1168 ! the idea is some stuff must go out the 'top'
1169 !
1170 do klcap = lcap,kte
1171 do i = its,ite
1172 taud_ls(i,klcap) = taud_ls(i,klcap) * factop
1173 taud_bl(i,klcap) = taud_bl(i,klcap) * factop
1174 enddo
1175 enddo
1176 !
1177 ! if the gravity wave drag would force a critical line
1178 ! in the lower ksmm1 layers during the next deltim timestep,
1179 ! then only apply drag until that critical line is reached.
1180 !
1181 do k = kts,kpblmax-1
1182 do i = its,ite
1183 if (k .le. kbl(i)) then
1184 if((taud_ls(i,k)+taud_bl(i,k)).ne.0.) &
1185 dtfac(i) = min(dtfac(i),abs(velco(i,k) &
1186 /(deltim*rcs*(taud_ls(i,k)+taud_bl(i,k)))))
1187 endif
1188 enddo
1189 enddo
1190 !
1191 do k = kts,kte
1192 do i = its,ite
1193 taud_ls(i,k) = taud_ls(i,k) * dtfac(i) * ls_taper
1194 taud_bl(i,k) = taud_bl(i,k) * dtfac(i) * ls_taper
1195 ! dtaux = taud(i,k) * xn(i)
1196 ! dtauy = taud(i,k) * yn(i)
1197 dtaux2d_ls(i,k) = taud_ls(i,k) * xn(i)
1198 dtauy2d_ls(i,k) = taud_ls(i,k) * yn(i)
1199 dtaux2d_bl(i,k) = taud_bl(i,k) * xn(i)
1200 dtauy2d_bl(i,k) = taud_bl(i,k) * yn(i)
1201 dudt(i,k) = dtaux2d_ls(i,k) + dtaux2d_bl(i,k) + dudt(i,k)
1202 dvdt(i,k) = dtauy2d_ls(i,k) + dtauy2d_bl(i,k) + dvdt(i,k)
1203 dusfc_ls(i) = dusfc_ls(i) + dtaux2d_ls(i,k) * del(i,k)
1204 dvsfc_ls(i) = dvsfc_ls(i) + dtauy2d_ls(i,k) * del(i,k)
1205 dusfc_bl(i) = dusfc_bl(i) + dtaux2d_bl(i,k) * del(i,k)
1206 dvsfc_bl(i) = dvsfc_bl(i) + dtauy2d_bl(i,k) * del(i,k)
1207 if ( gsl_diss_ht_opt .EQ. 1 ) then
1208 ! Calculate dissipation heating
1209 ! Initial kinetic energy (at t0-dt)
1210 eng0 = 0.5*( (rcs*u1(i,k))**2. + (rcs*v1(i,k))**2. )
1211 ! Kinetic energy after wave-breaking/flow-blocking
1212 eng1 = 0.5*( (rcs*(u1(i,k)+(dtaux2d_ls(i,k)+dtaux2d_bl(i,k))*deltim))**2. + &
1213 (rcs*(v1(i,k)+(dtauy2d_ls(i,k)+dtauy2d_bl(i,k))*deltim))**2. )
1214 ! Modify theta tendency
1215 dthdt(i,k) = dthdt(i,k) + max((eng0-eng1),0.0)/cp/deltim/prslk(i,k)
1216 end if
1217 enddo
1218 enddo
1220 ENDIF
1222 ! Finalize dusfc and dvsfc diagnoses
1223 do i = its,ite
1224 dusfc_ls(i) = (-1./g*rcs) * dusfc_ls(i)
1225 dvsfc_ls(i) = (-1./g*rcs) * dvsfc_ls(i)
1226 dusfc_bl(i) = (-1./g*rcs) * dusfc_bl(i)
1227 dvsfc_bl(i) = (-1./g*rcs) * dvsfc_bl(i)
1228 dusfc_ss(i) = (-1./g*rcs) * dusfc_ss(i)
1229 dvsfc_ss(i) = (-1./g*rcs) * dvsfc_ss(i)
1230 dusfc_fd(i) = (-1./g*rcs) * dusfc_fd(i)
1231 dvsfc_fd(i) = (-1./g*rcs) * dvsfc_fd(i)
1232 enddo
1234 !
1235 return
1236 end subroutine gwdo2d
1237 !-------------------------------------------------------------------
1238 end module module_bl_gwdo_gsl