| blob | 6e97351f8246b795fec171c9350a2ea75c4ceb97 |
1 !WRF:model_layer:physics\r
2 !\r
3 !####################e-epsilon scheme##################################\r
4 ! taken from the IPRC IRAM - Yuqing Wang, university of hawaii\r
5 ! added and improved by Chunxi Zhang and Yuqing Wang to wrf4.2, 2020\r
6 ! refenrences: Langland and Liou (1996, mwr, 124, 905-918)\r
7 ! Zhang et al. (2020, mwr, 148, 1121-1145)\r
8 !######################################################################\r
9 module module_bl_eeps\r
10 \r
11 use module_model_constants, only:rgas=>r_d, &\r
12 & t13=>Prandtl,ep1=>EP_1,ep2=>EP_2,grv=>g, &\r
13 & akv=>KARMAN,cp,rcp,xlv,xls,xlf,p1000mb\r
14 \r
15 real, private, parameter :: cs2=17.67,cs3=273.15,cs4=29.65\r
16 real, private, parameter :: ci2=22.514,ci3=273.16,ci4=0.0\r
17 real, private, parameter :: hgfr=233.15,t0=273.15\r
18 real, private, parameter :: epsl=1.e-20\r
19 real, private, parameter :: us2min=0.1\r
20 real, private, parameter :: akvmx = 1600.\r
21 real, private, parameter :: ricr=0.25,bt=5.0\r
22 real, private, parameter :: minpek = 1.0e-4\r
23 real, private, parameter :: minpep = 1.0e-6\r
24 real, private, parameter :: maxpek = 80.\r
25 real, private, parameter :: maxpep = 2.0\r
26 integer,private,parameter:: bl_eeps_tkeadv = 1\r
27 !---------------------------------------------------------------\r
28 ! ci--- five constants used in TKE tubrbulence closure scheme\r
29 !---------------------------------------------------------------\r
30 ! Detering and Etling (1986); Langland and Liou (1996)\r
31 ! real,parameter :: c1=1.35,c2=0.026,c3=1.13,c4=1.90,c5=0.77\r
32 ! Duynkerke and Driedonks (1987); Stubley and Rooney (1986)\r
33 real,parameter :: c1=1.35,c2=0.09,c3=1.44,c4=1.92,c5=0.77\r
34 ! Beljaars et al. (1987)\r
35 ! real,parameter :: c1=1.35,c2=0.032,c3=1.44,c4=1.92,c5=0.54\r
36 !\r
37 contains\r
38 !-------------------------------------------------------------------------------\r
39 !\r
40 subroutine eeps(u3d,v3d,t3d,qv3d,qc3d,qi3d,qr3d,qs3d,qg3d,p3d,pi3d, &\r
41 rublten,rvblten,rthblten, &\r
42 rqvblten,rqcblten,rqiblten, &\r
43 pek,pep, &\r
44 dz8w,psfc,qsfc,tsk,ust,rmol,wspd, &\r
45 xland,hfx,qfx, &\r
46 dt,dx,itimestep,exch_h,exch_m,pblh,kpbl, &\r
47 pek_adv,pep_adv, &\r
48 ids,ide, jds,jde, kds,kde, &\r
49 ims,ime, jms,jme, kms,kme, &\r
50 its,ite, jts,jte, kts,kte &\r
51 )\r
52 !-------------------------------------------------------------------------------\r
53 implicit none\r
54 !-------------------------------------------------------------------------------\r
55 !-- u3d 3d u-velocity interpolated to theta points (m/s)\r
56 !-- v3d 3d v-velocity interpolated to theta points (m/s)\r
57 !-- t3d temperature (k)\r
58 !-- qv3d 3d water vapor mixing ratio (kg/kg)\r
59 !-- qc3d 3d cloud mixing ratio (kg/kg)\r
60 !-- qi3d 3d ice mixing ratio (kg/kg)\r
61 !-- qr3d 3d rain mixing ratio (kg/kg)\r
62 !-- qs3d 3d snow mixing ratio (kg/kg)\r
63 !-- qg3d 3d graupel mixing ratio (kg/kg)\r
64 ! (note: if P_QI<PARAM_FIRST_SCALAR this should be zero filled)\r
65 !-- p3d 3d pressure (pa)\r
66 !-- pi3d 3d exner function (dimensionless)\r
67 !-- rublten u tendency due to\r
68 ! pbl parameterization (m/s/s)\r
69 !-- rvblten v tendency due to\r
70 ! pbl parameterization (m/s/s)\r
71 !-- rthblten theta tendency due to\r
72 ! pbl parameterization (K/s)\r
73 !-- rqvblten qv tendency due to\r
74 ! pbl parameterization (kg/kg/s)\r
75 !-- rqcblten qc tendency due to\r
76 ! pbl parameterization (kg/kg/s)\r
77 !-- rqiblten qi tendency due to\r
78 ! pbl parameterization (kg/kg/s)\r
79 !-- dz8w dz between full levels (m)\r
80 !-- pek 3d turbulent kinetic energy (TKE, m2/s2)\r
81 !-- pep 3d dissipation rate of TKE \r
82 !-- psfc pressure at the surface (pa)\r
83 !-- qsfc ground saturated mixing ratio\r
84 !-- ust u* in similarity theory (m/s)\r
85 !-- rmol inverse (1/L) of Monin-Obukhov length (m)\r
86 !-- xland land mask (1 for land, 2 for water)\r
87 !-- tsk surface skin temperature\r
88 !-- hfx upward heat flux at the surface (w/m^2)\r
89 !-- qfx upward moisture flux at the surface (kg/m^2/s)\r
90 !-- dt time step (s)\r
91 !-- ids start index for i in domain\r
92 !-- ide end index for i in domain\r
93 !-- jds start index for j in domain\r
94 !-- jde end index for j in domain\r
95 !-- kds start index for k in domain\r
96 !-- kde end index for k in domain\r
97 !-- ims start index for i in memory\r
98 !-- ime end index for i in memory\r
99 !-- jms start index for j in memory\r
100 !-- jme end index for j in memory\r
101 !-- kms start index for k in memory\r
102 !-- kme end index for k in memory\r
103 !-- its start index for i in tile\r
104 !-- ite end index for i in tile\r
105 !-- jts start index for j in tile\r
106 !-- jte end index for j in tile\r
107 !-- kts start index for k in tile\r
108 !-- kte end index for k in tile\r
109 !-------------------------------------------------------------------------------\r
110 !\r
111 integer, intent(in ) :: ids,ide, jds,jde, kds,kde, &\r
112 ims,ime, jms,jme, kms,kme, &\r
113 its,ite, jts,jte, kts,kte\r
114 \r
115 integer, intent(in ) :: itimestep\r
116 \r
117 real, intent(in ) :: dt, dx\r
118 !\r
119 !\r
120 real, dimension( ims:ime, kms:kme, jms:jme ) , &\r
121 intent(in ) :: qv3d, &\r
122 qc3d, &\r
123 qi3d, &\r
124 p3d, &\r
125 pi3d, &\r
126 t3d, &\r
127 dz8w, &\r
128 u3d, &\r
129 v3d\r
130 \r
131 real, dimension( ims:ime, kms:kme, jms:jme ) , &\r
132 optional , &\r
133 intent(in ) :: qr3d, &\r
134 qs3d, &\r
135 qg3d\r
136 \r
137 real, dimension( ims:ime, kms:kme, jms:jme ) , &\r
138 intent(inout) :: rublten, &\r
139 rvblten, &\r
140 rthblten, &\r
141 rqvblten, &\r
142 rqcblten, &\r
143 rqiblten\r
144 !\r
145 real, dimension( ims:ime, jms:jme ) , &\r
146 intent(in ) :: xland, &\r
147 hfx, &\r
148 qfx, &\r
149 psfc, &\r
150 qsfc, &\r
151 tsk, &\r
152 ust, &\r
153 rmol,&\r
154 wspd\r
155 !\r
156 real, dimension( ims:ime, kms:kme, jms:jme ) , &\r
157 intent(inout) :: pek, pep\r
158 \r
159 real, dimension( ims:ime, kms:kme, jms:jme ) , &\r
160 intent(inout) :: pek_adv,pep_adv\r
161 \r
162 real, dimension( ims:ime, kms:kme, jms:jme ) , &\r
163 intent(out) :: exch_h,exch_m\r
164 \r
165 real, dimension( ims:ime, jms:jme ) , &\r
166 intent(out ) :: pblh\r
167 \r
168 integer, dimension( ims:ime, jms:jme ) , &\r
169 intent(out ) :: kpbl\r
170 !\r
171 !local\r
172 real, dimension( its:ite, kts:kte+1):: zi,ghti \r
173 real, dimension( its:ite, kts:kte ) :: zl \r
174 \r
175 real, dimension( its:ite, kts:kte ) :: u1, &\r
176 v1, &\r
177 q1, &\r
178 q2, &\r
179 q3, &\r
180 q4, &\r
181 q5, &\r
182 q6, &\r
183 t1, &\r
184 ghtl, &\r
185 prsl\r
186 real, dimension( its:ite, kts:kte ) :: pek2d,pep2d\r
187 real, dimension( its:ite, kts:kte+1 ) :: exh2d,exm2d\r
188 real, dimension( its:ite) :: psfc1,ust1,rmol1, &\r
189 xland1,qsfc1,hfx1, &\r
190 qfx1,tsk1,pblh1,wspd1\r
191 integer, dimension( its:ite) :: kpbl1\r
192 \r
193 real :: rdelt\r
194 integer :: i,j,k,zz,im,kx,pp\r
195 \r
196 im=ite-its+1\r
197 kx=kte-kts+1\r
198 rdelt = 1./dt\r
199 \r
200 if(itimestep .eq. 1) then\r
201 pek_adv = pek\r
202 pep_adv = pep\r
203 end if\r
204 \r
205 do j=jts,jte\r
206 ! --------------- compute zi and zl -----------------------------------------\r
207 do i=its,ite\r
208 zi(i,kts)=0.0\r
209 enddo\r
210 !\r
211 do k=kts,kte\r
212 do i=its,ite\r
213 zi(i,k+1)=zi(i,k)+dz8w(i,k,j)\r
214 enddo\r
215 enddo\r
216 !\r
217 do k=kts,kte\r
218 do i=its,ite\r
219 zl(i,k)=0.5*(zi(i,k)+zi(i,k+1))\r
220 enddo\r
221 enddo\r
222 \r
223 ! reverse the vertical levels\r
224 pp = 0\r
225 do k=kts,kte\r
226 zz = kte-pp\r
227 do i=its,ite\r
228 u1(i,zz)=u3d(i,k,j)\r
229 v1(i,zz)=v3d(i,k,j)\r
230 t1(i,zz)=t3d(i,k,j)\r
231 q1(i,zz)=qv3d(i,k,j)/(1.+qv3d(i,k,j))\r
232 q2(i,zz)=qc3d(i,k,j)/(1.+qc3d(i,k,j))\r
233 q3(i,zz)=qi3d(i,k,j)/(1.+qi3d(i,k,j))\r
234 if(present(qr3d)) then\r
235 q4(i,zz)=qr3d(i,k,j)/(1.+qr3d(i,k,j))\r
236 else\r
237 q4(i,zz)=0.\r
238 end if\r
239 if(present(qs3d)) then\r
240 q5(i,zz)=qs3d(i,k,j)/(1.+qs3d(i,k,j))\r
241 else\r
242 q5(i,zz)=0.\r
243 end if\r
244 if(present(qg3d)) then\r
245 q6(i,zz)=qg3d(i,k,j)/(1.+qg3d(i,k,j))\r
246 else\r
247 q6(i,zz)=0.\r
248 end if\r
249 ghtl(i,zz)=zl(i,k)\r
250 prsl(i,zz) = p3d(i,k,j)\r
251 if(bl_eeps_tkeadv==1)then\r
252 if(k.ne.1) then \r
253 pek(i,k,j) = 0.5*(pek_adv(i,k,j)+pek_adv(i,k-1,j))\r
254 pep(i,k,j) = 0.5*(pep_adv(i,k,j)+pep_adv(i,k-1,j))\r
255 else\r
256 pek(i,k,j) = minpek\r
257 pep(i,k,j) = minpep\r
258 end if\r
259 end if\r
260 pek2d(i,zz)= min(maxpek,max(minpek,pek(i,k,j)))\r
261 pep2d(i,zz)= min(maxpep,max(minpep,pep(i,k,j)))\r
262 end do\r
263 pp = pp + 1\r
264 end do\r
265 \r
266 pp = 0\r
267 do k=kts,kte+1\r
268 zz = kte+1-pp\r
269 do i=its,ite\r
270 ghti(i,zz) = zi(i,k)\r
271 enddo\r
272 pp = pp + 1\r
273 enddo\r
274 \r
275 do i=its,ite\r
276 psfc1(i) = psfc(i,j)\r
277 ust1(i) = ust(i,j)\r
278 rmol1(i) = rmol(i,j)\r
279 wspd1(i) = wspd(i,j)\r
280 xland1(i)= xland(i,j)\r
281 qsfc1(i) = qsfc(i,j)\r
282 hfx1(i) = hfx(i,j)\r
283 qfx1(i) = qfx(i,j)\r
284 tsk1(i) = tsk(i,j)\r
285 end do\r
286 ! \r
287 call eeps2d(u1,v1,t1,q1,q2,q3,q4,q5,q6,prsl,ghtl,ghti,pek2d,pep2d &\r
288 ,psfc1,ust1,rmol1,wspd1,xland1,qsfc1,hfx1,qfx1,tsk1 &\r
289 ,dt,dx,itimestep,exh2d,exm2d,pblh1,kpbl1,im,kx)\r
290 !\r
291 pp = 0\r
292 do k=kts,kte\r
293 zz = kte-pp\r
294 do i=its,ite\r
295 rthblten(i,k,j)=(t1(i,zz)-t3d(i,k,j))/pi3d(i,k,j)*rdelt\r
296 rqvblten(i,k,j)=(q1(i,zz)/(1.-q1(i,zz))-qv3d(i,k,j))*rdelt\r
297 rqcblten(i,k,j)=(q2(i,zz)/(1.-q2(i,zz))-qc3d(i,k,j))*rdelt\r
298 rqiblten(i,k,j)=(q3(i,zz)/(1.-q3(i,zz))-qi3d(i,k,j))*rdelt\r
299 rublten(i,k,j) =(u1(i,zz)-u3d(i,k,j))*rdelt\r
300 rvblten(i,k,j) =(v1(i,zz)-v3d(i,k,j))*rdelt\r
301 \r
302 pek(i,k,j) = pek2d(i,zz)\r
303 pep(i,k,j) = pep2d(i,zz)\r
304 enddo\r
305 pp = pp + 1\r
306 enddo\r
307 \r
308 pp = 0\r
309 do k=kts,kte+1\r
310 zz = kte+1-pp\r
311 do i=its,ite\r
312 exch_h(i,k,j) = exh2d(i,zz)\r
313 exch_m(i,k,j) = exm2d(i,zz)\r
314 enddo\r
315 pp = pp + 1\r
316 enddo\r
317 \r
318 do i=its,ite\r
319 pblh(i,j) = pblh1(i)\r
320 kpbl(i,j) = kpbl1(i)\r
321 end do\r
322 \r
323 if(bl_eeps_tkeadv==1)then\r
324 do k=kts,kte\r
325 do i=its,ite\r
326 if(k.ne.kte) then\r
327 pek_adv(i,k,j) = 0.5*(pek(i,k,j)+pek(i,k+1,j))\r
328 pep_adv(i,k,j) = 0.5*(pep(i,k,j)+pep(i,k+1,j))\r
329 else\r
330 pek_adv(i,k,j) = 0.5*pek(i,k,j)\r
331 pep_adv(i,k,j) = 0.5*pep(i,k,j)\r
332 end if\r
333 end do\r
334 end do\r
335 end if\r
336 !\r
337 enddo\r
338 \r
339 !\r
340 end subroutine eeps\r
341 !\r
342 !-------------------------------------------------------------------------------\r
343 !\r
344 subroutine eeps2d(pu,pv,tz,pqv,pqc,pqi,pqr,pqs,pqg,prs,poz,zz,pek,pep, &\r
345 psfcpa,ust,rmol,wspd,xland,qsfc,hfx,qfx,tsk, &\r
346 dt,dx,itimestep,k_h,k_m,pblh,kpbl,lq,km)\r
347 !-------------------------------------------------------------------------------\r
348 implicit none\r
349 !-------------------------------------------------------------------------------\r
350 !\r
351 integer, intent(in ) :: itimestep,lq,km\r
352 real, intent(in ) :: dt, dx\r
353 !\r
354 real, dimension( lq ) , &\r
355 intent(in) :: ust, &\r
356 xland, &\r
357 hfx, &\r
358 qfx, &\r
359 rmol,&\r
360 wspd,&\r
361 psfcpa,qsfc,&\r
362 tsk\r
363 !\r
364 real, dimension( lq,km ),intent(inout) :: pu, &\r
365 pv, &\r
366 tz, &\r
367 pqv, &\r
368 pqc, &\r
369 pqi, &\r
370 pek, &\r
371 pep\r
372 \r
373 real, dimension( lq,km ),intent(in) :: prs, &\r
374 poz, &\r
375 pqr, &\r
376 pqs, &\r
377 pqg\r
378 real, dimension( lq,km+1 ),intent(in) :: zz\r
379 \r
380 real, dimension( lq,km+1 ),intent(out) :: k_h,k_m\r
381 real, dimension( lq ),intent(out) :: pblh \r
382 integer, dimension( lq ),intent(out) :: kpbl\r
383 \r
384 ! Local Vars\r
385 real, dimension( lq,km ) :: fac,psp,pt,ztv,szv,disht,up2\r
386 real, dimension( lq,km ) :: qvs,st,sh,rich,richb,dum2,dzz,doz\r
387 real, dimension( lq,km ) :: am,src,src1,ax,bx,cx,ac,ab,ba,alt\r
388 real, dimension( lq,km-1 ) :: ax1,bx1,cx1,yy1,amt\r
389 real, dimension( lq,km+1 ) :: akm,akh,pnt\r
390 \r
391 real, dimension( lq) :: tkm,zkm,qkm,us2,us,dens,wstar2,sflux\r
392 real, dimension( lq) :: rm2,sfa,hs,hq\r
393 real, dimension( lq) :: gamat,gamaq\r
394 real, dimension( lq,km) :: sgk1,sgk2 \r
395 integer :: j,k,mdt,lvl\r
396 real :: dtt,dum,dum0,alh,eqp,qvsw,qvsi,faf,rdz,du,dv,dss\r
397 real :: tk,qvk,alk,af,ag,aks,aco\r
398 real :: duj,alv0,richf,rch,alv,ak,cpm\r
399 real :: govrth,bfx0,wstar3,coef\r
400 \r
401 mdt = min(6,max(3,int(dt/10.0+0.1)))\r
402 dtt = dt/mdt\r
403 \r
404 ! dzz for dz between full levels, doz for dz between half levels\r
405 do k=1,km\r
406 do j=1,lq\r
407 dzz(j,k) = zz(j,k)-zz(j,k+1)\r
408 alt(j,k) = dt/dzz(j,k)\r
409 up2(j,k)=grv/max(us2min,pu(j,k)*pu(j,k)+pv(j,k)*pv(j,k))\r
410 enddo\r
411 enddo\r
412 \r
413 do k=1,km-1\r
414 do j=1,lq\r
415 doz(j,k) = poz(j,k)-poz(j,k+1)\r
416 end do\r
417 end do\r
418 !\r
419 do j=1,lq\r
420 doz(j,km)= 2.0*poz(j,km) ! not used\r
421 ! wstar2(j)= 0.0 ! it will be assigned values later\r
422 enddo\r
423 !\r
424 do k=1,km-1\r
425 do j=1,lq\r
426 amt(j,k) = dtt/doz(j,k)\r
427 enddo\r
428 enddo\r
429 \r
430 do k=1,km\r
431 do j=1,lq\r
432 sgk1(j,k) = 0.1\r
433 sgk2(j,k) = 0.01\r
434 end do\r
435 end do\r
436 !\r
437 do k=1,km\r
438 do j=1,lq\r
439 psp(j,k)=(prs(j,k)/p1000mb)**rcp\r
440 pt(j,k)=tz(j,k)/psp(j,k)\r
441 dum=1.0+ep1*pqv(j,k)\r
442 ztv(j,k)=tz(j,k)*dum\r
443 szv(j,k)=pt(j,k)*dum\r
444 fac(j,k) = pqc(j,k) + pqi(j,k) + pqr(j,k) + pqs(j,k) + pqg(j,k)\r
445 cx(j,k)=1.0+pqv(j,k)+fac(j,k)\r
446 \r
447 alh=3.1484e6-2.37e3*tz(j,k)\r
448 eqp=611.2*exp(cs2*(tz(j,k)-cs3)/(tz(j,k)-cs4))\r
449 eqp=min(0.5*prs(j,k),eqp)\r
450 qvsw=0.622*eqp/(prs(j,k)-eqp)\r
451 eqp=611.0*exp(ci2*(tz(j,k)-ci3)/(tz(j,k)-ci4))\r
452 eqp=min(0.5*prs(j,k),eqp)\r
453 qvsi=0.622*eqp/(prs(j,k)-eqp)\r
454 if(tz(j,k).ge.t0) then\r
455 faf=0.0\r
456 else if(tz(j,k).le.hgfr) then\r
457 faf=1.0\r
458 else\r
459 if(pqi(j,k).le.epsl) then\r
460 faf=0.0\r
461 else if(pqc(j,k).le.epsl) then\r
462 faf=1.0\r
463 else\r
464 faf=pqi(j,k)/(pqc(j,k)+pqi(j,k))\r
465 endif\r
466 endif\r
467 qvs(j,k)=(1.0-faf)*qvsw+faf*qvsi\r
468 ax(j,k)=(1.0-faf)*alh+faf*xls\r
469 end do\r
470 end do\r
471 \r
472 !\r
473 ! to calculate the Richardson number and vertical shear of\r
474 ! the horizontal wind above the surface layer\r
475 !\r
476 do k=1,km-1\r
477 do j=1,lq\r
478 rdz=1.0/(poz(j,k)-poz(j,k+1))\r
479 du=pu(j,k)-pu(j,k+1)\r
480 dv=pv(j,k)-pv(j,k+1)\r
481 sh(j,k)=(du*du+dv*dv)*(rdz*rdz)\r
482 if((pqv(j,k).lt.qvs(j,k)).or.(pqv(j,k+1).lt.qvs(j,k+1))) then\r
483 dss=log(szv(j,k)/szv(j,k+1))\r
484 st(j,k)=grv*dss*rdz\r
485 else\r
486 dss=log(pt(j,k)/pt(j,k+1))\r
487 tk=0.5*(tz(j,k+1)+tz(j,k))\r
488 qvk=0.5*(pqv(j,k+1)+pqv(j,k))\r
489 alk=0.5*(ax(j,k+1)+ax(j,k))\r
490 af=alk*qvk/(rgas*tk)\r
491 ag=alk/(cp*tk)\r
492 aks=ep2*af*ag\r
493 aco=(1.0+af)/(1.0+aks)\r
494 st(j,k)=(aco*(dss+ag*(pqv(j,k)-pqv(j,k+1))) &\r
495 -(cx(j,k)-cx(j,k+1)))*grv*rdz\r
496 endif\r
497 rich(j,k)=st(j,k)/max(1.0e-7,sh(j,k))\r
498 end do\r
499 end do\r
500 \r
501 ! calculate first guess pblh\r
502 do k=1,km\r
503 do j=1,lq\r
504 richb(j,k)=poz(j,k)*(szv(j,k)/szv(j,km)-1.0)*up2(j,k)\r
505 enddo\r
506 enddo\r
507 !\r
508 do j=1,lq\r
509 pblh(j)=poz(j,km)\r
510 kpbl(j)=1\r
511 lvl=km\r
512 do k=km-1,5,-1\r
513 if(richb(j,k).ge.ricr) then\r
514 lvl=k\r
515 exit\r
516 endif\r
517 enddo\r
518 \r
519 if(lvl.lt.km) then\r
520 dum=(poz(j,lvl)-poz(j,lvl+1))/(richb(j,lvl)-richb(j,lvl+1))\r
521 pblh(j)=max(poz(j,lvl+1),poz(j,lvl)-(richb(j,lvl)-ricr)*dum)\r
522 pblh(j)=min(pblh(j),poz(j,lvl))\r
523 endif\r
524 end do\r
525 \r
526 ! preparing for tke/ep calculations\r
527 !\r
528 do j=1,lq\r
529 tkm(j) = tz(j,km)\r
530 zkm(j) = poz(j,km)\r
531 qkm(j) = pqv(j,km)\r
532 sfa(j) = tsk(j)*(p1000mb/psfcpa(j))**rcp\r
533 dens(j) = psfcpa(j)/(rgas*tkm(j)*(1.+ep1*qkm(j)))\r
534 cpm = cp * (1.+0.81*qkm(j))\r
535 \r
536 rm2(j)=ust(j)*ust(j)/wspd(j)\r
537 hs(j)=hfx(j)/(dens(j)*cpm)\r
538 hq(j)=qfx(j)/dens(j)\r
539 !\r
540 sh(j,km)=0.0\r
541 st(j,km)=0.0\r
542 duj=ust(j)*ust(j)\r
543 pep(j,km)=duj*ust(j)/(akv*zkm(j))\r
544 pep(j,km) = min(maxpep,max(minpep,pep(j,km)))\r
545 govrth = grv/pt(j,km)\r
546 sflux(j) = hfx(j)/dens(j)/cpm + &\r
547 qfx(j)/dens(j)*ep1*pt(j,km)\r
548 bfx0 = max(sflux(j),0.)\r
549 wstar3 = govrth*bfx0*pblh(j)\r
550 wstar2(j) = (wstar3**2)**t13\r
551 \r
552 rich(j,km) = grv*zkm(j)* &\r
553 log(szv(j,km)/(sfa(j)*(1.+ep1*qsfc(j))))/(wspd(j)**2.)\r
554 \r
555 if(rmol(j).lt.0. .and. sflux(j).gt.0) then\r
556 pek(j,km)=3.75*duj+0.2*wstar2(j)+duj*((-zkm(j)*rmol(j))**2.)**t13\r
557 else\r
558 pek(j,km)=3.75*duj\r
559 endif\r
560 pek(j,km) = min(maxpek,max(minpek,pek(j,km)))\r
561 \r
562 end do\r
563 \r
564 ! to recalculate pblh here\r
565 do j=1,lq\r
566 if(rmol(j).lt.0. .and. sflux(j).gt.0) then\r
567 dss=szv(j,km)+min(5.0,bt*sflux(j)/(max(0.1,sqrt(wstar2(j)))))\r
568 do k=1,km\r
569 richb(j,k)=poz(j,k)*(szv(j,k)/dss-1.0)*up2(j,k)\r
570 end do\r
571 \r
572 pblh(j)=poz(j,km)\r
573 lvl=km\r
574 do k=km-1,5,-1\r
575 if(richb(j,k).ge.ricr) then\r
576 lvl=k\r
577 exit\r
578 end if\r
579 enddo\r
580 \r
581 if(lvl.lt.km) then\r
582 dum=(poz(j,lvl)-poz(j,lvl+1))/(richb(j,lvl)-richb(j,lvl+1))\r
583 pblh(j)=max(poz(j,lvl+1),poz(j,lvl)-(richb(j,lvl)-ricr)*dum)\r
584 pblh(j)=min(pblh(j),poz(j,lvl))\r
585 kpbl(j)=km-lvl+1\r
586 end if\r
587 end if\r
588 end do\r
589 \r
590 ! to add the countergradient term\r
591 do j=1,lq\r
592 if(rmol(j).lt.0.0 .and. sflux(j).gt.0) then\r
593 du=max(0.1,sqrt(wstar2(j)))\r
594 gamat(j)=max(bt*hs(j)/(pblh(j)*du),0.)\r
595 gamaq(j)=max(bt*hq(j)/(pblh(j)*du),0.)\r
596 else\r
597 gamat(j)=0.0\r
598 gamaq(j)=0.0\r
599 endif\r
600 end do\r
601 !\r
602 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
603 ! time integration of tke and ep equations\r
604 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
605 do k=1,km\r
606 do j=1,lq\r
607 richf=0.1912408\r
608 if(rich(j,k).lt.0.195) then\r
609 richf=0.6588*(rich(j,k)+0.1776-sqrt(rich(j,k)*rich(j,k)- &\r
610 0.3221*rich(j,k)+0.03156))\r
611 richf=max(-1000.0,min(0.1912408,richf))\r
612 endif\r
613 dum2(j,k)=1.122346\r
614 if(richf.lt.0.16) dum2(j,k)=1.318*(0.2231-richf) &\r
615 /(0.2341-richf)\r
616 end do\r
617 end do\r
618 !\r
619 lvl=1\r
620 do while(lvl .le. mdt)\r
621 \r
622 do k=1,km-1\r
623 do j=1,lq\r
624 akm(j,k)=c2*pek(j,k)*pek(j,k)/pep(j,k)\r
625 akh(j,k)=dum2(j,k)*akm(j,k)\r
626 akm(j,k)=sgk1(j,k)+akm(j,k)\r
627 akh(j,k)=sgk2(j,k)+akh(j,k)\r
628 akm(j,k)=min(akvmx,akm(j,k))\r
629 akh(j,k)=min(akvmx,akh(j,k))\r
630 if(zz(j,k).gt.pblh(j))then\r
631 src(j,k)=akm(j,k)*sh(j,k)-akh(j,k)*st(j,k)\r
632 src1(j,k)=max(akm(j,k)*sh(j,k),akm(j,k)*sh(j,k)-akh(j,k)*st(j,k))\r
633 else\r
634 src(j,k)=akm(j,k)*sh(j,k)-akh(j,k)*(st(j,k) &\r
635 -2.0*grv*gamat(j)/(szv(j,k+1)+szv(j,k)))\r
636 src1(j,k)=max(akm(j,k)*sh(j,k),akm(j,k)*sh(j,k)-akh(j,k)*(st(j,k) &\r
637 -2.0*grv*gamat(j)/(szv(j,k+1)+szv(j,k))))\r
638 end if\r
639 end do\r
640 end do\r
641 !\r
642 do j=1,lq\r
643 akm(j,km)=c2*pek(j,km)*pek(j,km)/pep(j,km)\r
644 akh(j,km)=dum2(j,km)*akm(j,km)\r
645 end do\r
646 !\r
647 ! to calculate the vertical diffusion coefficients at u,v,levels\r
648 !\r
649 do k=2,km\r
650 do j=1,lq\r
651 am(j,k)=0.5*(akm(j,k-1)+akm(j,k))/dzz(j,k)\r
652 end do\r
653 end do\r
654 \r
655 do j=1,lq\r
656 am(j,1)=0.5*akm(j,1)/dzz(j,1)\r
657 end do\r
658 !\r
659 ! to calculate the cofficients for the triangle matrix\r
660 ! for tke dissipation and solve the matrix\r
661 !\r
662 do k=1,km-1\r
663 do j=1,lq\r
664 dum=(c3*src1(j,k)-c4*pep(j,k))*dtt/pek(j,k)\r
665 ax1(j,k)=-c5*amt(j,k)*am(j,k+1)\r
666 cx1(j,k)=-c5*amt(j,k)*am(j,k)\r
667 if(dum.le.0.5) then\r
668 bx1(j,k)=1.0+c5*amt(j,k)*(am(j,k)+am(j,k+1))-dum\r
669 yy1(j,k)=pep(j,k)\r
670 else\r
671 bx1(j,k)=1.0+c5*amt(j,k)*(am(j,k)+am(j,k+1))\r
672 yy1(j,k)=pep(j,k)+dum*pep(j,k)\r
673 endif\r
674 end do\r
675 end do\r
676 !\r
677 do j=1,lq\r
678 cx1(j,1)=0.0\r
679 ax1(j,km-1)=0.0\r
680 yy1(j,km-1)=yy1(j,km-1)+c5*amt(j,km-1)*am(j,km)*pep(j,km)\r
681 end do\r
682 !\r
683 call tridiag(ax1,bx1,cx1,yy1,lq,km-1)\r
684 !\r
685 ! to calculate the cofficients for the triangle matrix\r
686 ! for tke and solve the matrix\r
687 !\r
688 do k=1,km-1\r
689 do j=1,lq\r
690 pep(j,k)=min(maxpep,max(minpep, yy1(j,k)))\r
691 dum=(src(j,k)-pep(j,k))*dtt\r
692 ax1(j,k)=-c1*amt(j,k)*am(j,k+1)\r
693 bx1(j,k)=1.0+c1*amt(j,k)*(am(j,k)+am(j,k+1))\r
694 cx1(j,k)=-c1*amt(j,k)*am(j,k)\r
695 yy1(j,k)=pek(j,k)+dum\r
696 end do\r
697 end do\r
698 \r
699 do j=1,lq\r
700 cx1(j,1)=0.0\r
701 ax1(j,km-1)=0.0\r
702 yy1(j,km-1)=yy1(j,km-1)+c1*amt(j,km-1)*am(j,km)*pek(j,km)\r
703 end do\r
704 !\r
705 call tridiag(ax1,bx1,cx1,yy1,lq,km-1)\r
706 !\r
707 do k=1,km-1\r
708 do j=1,lq\r
709 pek(j,k)=min(maxpek,max(minpek, yy1(j,k)))\r
710 end do\r
711 end do\r
712 \r
713 lvl = lvl + 1\r
714 end do ! end while\r
715 !\r
716 !-------------------------------------------------------------\r
717 ! to calculate the vertical diffusion coeficient for momentum \r
718 ! heat, moisture\r
719 !-------------------------------------------------------------\r
720 ! note: akm and akh at the first level (km) will not be used \r
721 do k=2,km\r
722 do j=1,lq\r
723 akm(j,k)=c2*pek(j,k-1)*pek(j,k-1)/pep(j,k-1)\r
724 akh(j,k)=dum2(j,k)*akm(j,k)\r
725 akm(j,k)=sgk1(j,k)+akm(j,k)\r
726 akh(j,k)=sgk2(j,k)+akh(j,k)\r
727 akm(j,k)=min(akvmx,akm(j,k))\r
728 akh(j,k)=min(akvmx,akh(j,k))\r
729 pnt(j,k)=akh(j,k)\r
730 if(zz(j,k).gt.pblh(j))pnt(j,k)=0.\r
731 akm(j,k)=akm(j,k)/doz(j,k-1)\r
732 akh(j,k)=akh(j,k)/doz(j,k-1)\r
733 end do\r
734 end do\r
735 !\r
736 do j=1,lq\r
737 akm(j,1)=0.0\r
738 akh(j,1)=0.0\r
739 pnt(j,1)=0.0\r
740 akm(j,km+1)=0.0\r
741 akh(j,km+1)=0.0\r
742 pnt(j,km+1)=0.0\r
743 end do\r
744 !\r
745 ! Momentum fluxes calculation\r
746 !\r
747 do k=1,km\r
748 do j=1,lq\r
749 ax(j,k)=-alt(j,k)*akm(j,k+1)\r
750 bx(j,k)=1.0+alt(j,k)*(akm(j,k)+akm(j,k+1))\r
751 cx(j,k)=-alt(j,k)*akm(j,k)\r
752 ba(j,k)=ax(j,k)\r
753 ab(j,k)=bx(j,k)\r
754 ac(j,k)=cx(j,k)\r
755 end do\r
756 end do\r
757 \r
758 do j=1,lq\r
759 bx(j,km) = bx(j,km)+alt(j,km)*rm2(j)\r
760 ab(j,km)=bx(j,km)\r
761 end do\r
762 \r
763 !\r
764 call tridiag(ba,ab,ac,pu,lq,km)\r
765 !\r
766 call tridiag(ax,bx,cx,pv,lq,km)\r
767 !\r
768 ! Heat and moisture fluxes calculation\r
769 !\r
770 do k=1,km\r
771 do j=1,lq\r
772 ax(j,k)=-alt(j,k)*akh(j,k+1)\r
773 bx(j,k)=1.0+alt(j,k)*(akh(j,k)+akh(j,k+1))\r
774 cx(j,k)=-alt(j,k)*akh(j,k)\r
775 ba(j,k)=ax(j,k)\r
776 ab(j,k)=bx(j,k)\r
777 ac(j,k)=cx(j,k)\r
778 dum0=0.0 \r
779 if(k.gt.1) dum0=pep(j,k-1)\r
780 dum=pep(j,k)\r
781 if(dum0.le.1.e-6) dum0=0.0\r
782 if(dum.le.1.e-6) dum=0.0\r
783 disht(j,k)=0.5*(dum0+dum)/(cp*(1.+.81*pqv(j,k)))\r
784 pt(j,k)=pt(j,k)+disht(j,k)*dt/psp(j,k)+alt(j,k)* &\r
785 gamat(j)*(pnt(j,k+1)-pnt(j,k))\r
786 end do\r
787 end do\r
788 !\r
789 do j=1,lq\r
790 pt(j,km)=pt(j,km)+alt(j,km)*hs(j)\r
791 end do\r
792 !\r
793 call tridiag(ba,ab,ac,pt,lq,km)\r
794 !\r
795 do k=1,km\r
796 do j=1,lq\r
797 ba(j,k)=ax(j,k)\r
798 ab(j,k)=bx(j,k)\r
799 ac(j,k)=cx(j,k)\r
800 end do\r
801 end do\r
802 \r
803 do k=1,km\r
804 do j=1,lq\r
805 pqv(j,k)=pqv(j,k)+alt(j,k)*gamaq(j)* &\r
806 (pnt(j,k+1)-pnt(j,k))\r
807 end do\r
808 end do\r
809 \r
810 do j=1,lq\r
811 pqv(j,km)=pqv(j,km)+alt(j,km)*hq(j)\r
812 end do\r
813 !\r
814 call tridiag(ba,ab,ac,pqv,lq,km)\r
815 !\r
816 do k=1,km\r
817 do j=1,lq\r
818 ba(j,k)=ax(j,k)\r
819 ab(j,k)=bx(j,k)\r
820 ac(j,k)=cx(j,k)\r
821 end do\r
822 end do\r
823 \r
824 call tridiag(ba,ab,ac,pqc,lq,km)\r
825 !\r
826 call tridiag(ax,bx,cx,pqi,lq,km)\r
827 !\r
828 do k=1,km\r
829 do j=1,lq\r
830 tz(j,k)=pt(j,k)*psp(j,k)\r
831 pqv(j,k)=max(epsl,pqv(j,k))\r
832 pqc(j,k)=max(epsl,pqc(j,k))\r
833 pqi(j,k)=max(epsl,pqi(j,k))\r
834 !-----------------------------------------------------\r
835 ! immediate melting of cloud ice below freezing level\r
836 !-----------------------------------------------------\r
837 if(tz(j,k).ge.t0.and.pqi(j,k).gt.epsl) then\r
838 if(pqc(j,k).le.epsl) then\r
839 pqc(j,k)=pqi(j,k)\r
840 else\r
841 pqc(j,k)=pqc(j,k)+pqi(j,k)\r
842 endif\r
843 cpm=cp*(1.0+0.81*pqv(j,k))\r
844 alh=3.1484e6-2.37e3*tz(j,k)\r
845 tz(j,k)=tz(j,k)-(xls-alh)*pqi(j,k)/cpm\r
846 pqi(j,k)=epsl\r
847 endif\r
848 end do\r
849 end do\r
850 \r
851 do k=2,km\r
852 do j=1,lq\r
853 k_h(j,k)=akh(j,k)*doz(j,k-1)\r
854 k_m(j,k)=akm(j,k)*doz(j,k-1)\r
855 end do\r
856 end do\r
857 \r
858 do j=1,lq\r
859 k_h(j,1)=0.\r
860 k_m(j,1)=0.\r
861 k_m(j,km+1)= 0. !c2*pek(j,km)*pek(j,km)/pep(j,km) ! we don't need it\r
862 k_h(j,km+1)= 0. !dum2(j,km)*k_m(j,km+1) ! we don't need it\r
863 end do\r
864 \r
865 end subroutine eeps2d\r
866 \r
867 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
868 !\r
869 ! ==================================================================\r
870 subroutine tridiag(a,b,c,d,lq,km)\r
871 !! to solve system of linear eqs on tridiagonal matrix n times n\r
872 !! after Peaceman and Rachford, 1955\r
873 !! a,b,c,d - are vectors of order n \r
874 !! a,b,c - are coefficients on the LHS\r
875 !! d - is initially RHS on the output becomes a solution vector\r
876 \r
877 !-------------------------------------------------------------------\r
878 implicit none\r
879 INTEGER, INTENT(in):: lq,km\r
880 REAL, DIMENSION(lq,km), INTENT(in) :: a,b\r
881 REAL, DIMENSION(lq,km), INTENT(inout) :: c,d\r
882 \r
883 INTEGER :: j,k\r
884 REAL :: p\r
885 REAL, DIMENSION(lq,km) :: q\r
886 \r
887 do j=1,lq\r
888 c(j,1)=0.\r
889 q(j,km)=-c(j,km)/b(j,km)\r
890 d(j,km)=d(j,km)/b(j,km)\r
891 end do\r
892 \r
893 DO k=km-1,1,-1\r
894 do j=1,lq\r
895 p=1./(b(j,k)+a(j,k)*q(j,k+1))\r
896 q(j,k)=-c(j,k)*p\r
897 d(j,k)=(d(j,k)-a(j,k)*d(j,k+1))*p\r
898 end do\r
899 ENDDO\r
900 \r
901 DO k=2,km\r
902 do j=1,lq\r
903 d(j,k)=d(j,k)+q(j,k)*d(j,k-1)\r
904 end do\r
905 ENDDO\r
906 \r
907 end subroutine tridiag\r
908 \r
909 !-------------------------------------------------------------------------------\r
910 subroutine eepsinit(rublten,rvblten,rthblten,rqvblten, &\r
911 rqcblten,rqiblten,p_qi,p_first_scalar,pek,pep, &\r
912 restart, allowed_to_read, &\r
913 ids, ide, jds, jde, kds, kde, &\r
914 ims, ime, jms, jme, kms, kme, &\r
915 its, ite, jts, jte, kts, kte )\r
916 !-------------------------------------------------------------------------------\r
917 implicit none\r
918 !-------------------------------------------------------------------------------\r
919 !\r
920 logical , intent(in) :: restart, allowed_to_read\r
921 integer , intent(in) :: ids, ide, jds, jde, kds, kde, &\r
922 ims, ime, jms, jme, kms, kme, &\r
923 its, ite, jts, jte, kts, kte\r
924 integer , intent(in) :: p_qi,p_first_scalar\r
925 real , dimension( ims:ime , kms:kme , jms:jme ), intent(out) :: &\r
926 rublten, &\r
927 rvblten, &\r
928 rthblten, &\r
929 rqvblten, &\r
930 rqcblten, &\r
931 rqiblten\r
932 real, dimension( ims:ime, kms:kme, jms:jme ),intent(out) :: pek, pep\r
933 integer :: i, j, k, itf, jtf, ktf\r
934 !\r
935 jtf = min0(jte,jde-1)\r
936 ktf = min0(kte,kde-1)\r
937 itf = min0(ite,ide-1)\r
938 !\r
939 if(.not.restart)then\r
940 do j = jts,jtf\r
941 do k = kts,ktf\r
942 do i = its,itf\r
943 rublten(i,k,j) = 0.\r
944 rvblten(i,k,j) = 0.\r
945 rthblten(i,k,j) = 0.\r
946 rqvblten(i,k,j) = 0.\r
947 rqcblten(i,k,j) = 0.\r
948 pek(i,k,j) = minpek\r
949 pep(i,k,j) = minpep\r
950 enddo\r
951 enddo\r
952 enddo\r
953 endif\r
954 !\r
955 if (p_qi .ge. p_first_scalar .and. .not.restart) then\r
956 do j = jts,jtf\r
957 do k = kts,ktf\r
958 do i = its,itf\r
959 rqiblten(i,k,j) = 0.\r
960 enddo\r
961 enddo\r
962 enddo\r
963 endif\r
964 !\r
965 end subroutine eepsinit\r
966 !------------------------------------------------------------------------------- \r
967 end module module_bl_eeps\r