| blob | 58b5f7b7393dad6ea7cc05d1d89456562c2040a2 |
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 ! \r
10 ! Updates since the original release\r
11 ! 03/2023\r
12 ! Added the cloud top cooling effect for stratocumulus, which is similar\r
13 ! to the one used in the MYNN PBL scheme. The critial step is to\r
14 ! idendity stratocumulus. We used the method by Marquet and Bechtold (2020)\r
15 ! Added addtinal source term for stable boundary layer in the epsilon \r
16 ! equation to represent the energy drain on the eddy scale\r
17 ! (Zeng et al. (2020))\r
18 ! Fixed a bug for the TKE global varaible. A new variable pek_pbl is defined \r
19 ! instead of using tke_pbl in Register.EM_COMMON. Now pek_pbl,pep_pbl,pek_adv\r
20 ! and pep_adv are all on the half levels. \r
21 ! Addded a few tweaks\r
22 !\r
23 !######################################################################\r
24 module module_bl_eeps\r
25 \r
26 use module_model_constants, only:rgas=>r_d,rv=>r_v, &\r
27 & t13=>Prandtl,ep1=>EP_1,ep2=>EP_2,grv=>g, &\r
28 & akv=>KARMAN,cp,rcp,xlv,xls,xlf,p1000mb\r
29 \r
30 real, private, parameter :: cs2=17.67,cs3=273.15,cs4=29.65\r
31 real, private, parameter :: ci2=22.514,ci3=273.16,ci4=0.0\r
32 real, private, parameter :: hgfr=233.15,t0=273.15\r
33 real, private, parameter :: epsl=1.e-20\r
34 real, private, parameter :: us2min=0.1\r
35 real, private, parameter :: akvmx = 1600.\r
36 real, private, parameter :: ricr=0.25,bt=5.0\r
37 real, private, parameter :: minpek = 1.0e-4\r
38 real, private, parameter :: minpep = 1.0e-7\r
39 real, private, parameter :: maxpek = 160.\r
40 real, private, parameter :: maxpep = 25.0\r
41 real, private, parameter :: zfmin = 0.01\r
42 integer,private,parameter:: bl_eeps_tkeadv = 1\r
43 !Adding top-down diffusion driven by cloud-top radiative cooling\r
44 integer,private,parameter:: bl_eeps_topdown = 1\r
45 !---------------------------------------------------------------\r
46 ! ci--- five constants used in TKE tubrbulence closure scheme\r
47 !---------------------------------------------------------------\r
48 ! Detering and Etling (1986); Langland and Liou (1996)\r
49 ! real,parameter :: c1=1.35,c2=0.026,c3=1.13,c4=1.90,c5=0.77\r
50 ! Duynkerke and Driedonks (1987); Stubley and Rooney (1986)\r
51 real,parameter :: c1=1.35,c2=0.09,c3=1.44,c4=1.92,c5=0.77\r
52 ! Beljaars et al. (1987)\r
53 ! real,parameter :: c1=1.35,c2=0.032,c3=1.44,c4=1.92,c5=0.54\r
54 !\r
55 contains\r
56 !-------------------------------------------------------------------------------\r
57 !\r
58 subroutine eeps(u3d,v3d,t3d,qv3d,qc3d,qi3d,qr3d,qs3d,qg3d,p3d,pi3d, &\r
59 rho3d,rthraten,rublten,rvblten,rthblten, &\r
60 rqvblten,rqcblten,rqiblten, &\r
61 pek,pep, &\r
62 dz8w,psfc,qsfc,tsk,ust,rmol,wspd, &\r
63 xland,hfx,qfx, &\r
64 dt,dx,itimestep,exch_h,exch_m,pblh,kpbl, &\r
65 pek_adv,pep_adv, &\r
66 ids,ide, jds,jde, kds,kde, &\r
67 ims,ime, jms,jme, kms,kme, &\r
68 its,ite, jts,jte, kts,kte &\r
69 )\r
70 !-------------------------------------------------------------------------------\r
71 implicit none\r
72 !-------------------------------------------------------------------------------\r
73 !-- u3d 3d u-velocity interpolated to theta points (m/s)\r
74 !-- v3d 3d v-velocity interpolated to theta points (m/s)\r
75 !-- t3d temperature (k)\r
76 !-- qv3d 3d water vapor mixing ratio (kg/kg)\r
77 !-- qc3d 3d cloud mixing ratio (kg/kg)\r
78 !-- qi3d 3d ice mixing ratio (kg/kg)\r
79 !-- qr3d 3d rain mixing ratio (kg/kg)\r
80 !-- qs3d 3d snow mixing ratio (kg/kg)\r
81 !-- qg3d 3d graupel mixing ratio (kg/kg)\r
82 ! (note: if P_QI<PARAM_FIRST_SCALAR this should be zero filled)\r
83 !-- p3d 3d pressure (pa)\r
84 !-- pi3d 3d exner function (dimensionless)\r
85 !-- rho3d 3d air density (kg/m^3)\r
86 !-- rthraten theta tendency due to\r
87 ! radiation parameterization (K/s)\r
88 !-- rublten u tendency due to\r
89 ! pbl parameterization (m/s/s)\r
90 !-- rvblten v tendency due to\r
91 ! pbl parameterization (m/s/s)\r
92 !-- rthblten theta tendency due to\r
93 ! pbl parameterization (K/s)\r
94 !-- rqvblten qv tendency due to\r
95 ! pbl parameterization (kg/kg/s)\r
96 !-- rqcblten qc tendency due to\r
97 ! pbl parameterization (kg/kg/s)\r
98 !-- rqiblten qi tendency due to\r
99 ! pbl parameterization (kg/kg/s)\r
100 !-- dz8w dz between full levels (m)\r
101 !-- pek 3d turbulent kinetic energy (TKE, m2/s2)\r
102 !-- pep 3d dissipation rate of TKE \r
103 !-- psfc pressure at the surface (pa)\r
104 !-- qsfc ground saturated mixing ratio\r
105 !-- ust u* in similarity theory (m/s)\r
106 !-- rmol inverse (1/L) of Monin-Obukhov length (m)\r
107 !-- xland land mask (1 for land, 2 for water)\r
108 !-- tsk surface skin temperature\r
109 !-- hfx upward heat flux at the surface (w/m^2)\r
110 !-- qfx upward moisture flux at the surface (kg/m^2/s)\r
111 !-- dt time step (s)\r
112 !-- ids start index for i in domain\r
113 !-- ide end index for i in domain\r
114 !-- jds start index for j in domain\r
115 !-- jde end index for j in domain\r
116 !-- kds start index for k in domain\r
117 !-- kde end index for k in domain\r
118 !-- ims start index for i in memory\r
119 !-- ime end index for i in memory\r
120 !-- jms start index for j in memory\r
121 !-- jme end index for j in memory\r
122 !-- kms start index for k in memory\r
123 !-- kme end index for k in memory\r
124 !-- its start index for i in tile\r
125 !-- ite end index for i in tile\r
126 !-- jts start index for j in tile\r
127 !-- jte end index for j in tile\r
128 !-- kts start index for k in tile\r
129 !-- kte end index for k in tile\r
130 !-------------------------------------------------------------------------------\r
131 !\r
132 integer, intent(in ) :: ids,ide, jds,jde, kds,kde, &\r
133 ims,ime, jms,jme, kms,kme, &\r
134 its,ite, jts,jte, kts,kte\r
135 \r
136 integer, intent(in ) :: itimestep\r
137 \r
138 real, intent(in ) :: dt, dx\r
139 !\r
140 !\r
141 real, dimension( ims:ime, kms:kme, jms:jme ) , &\r
142 intent(in ) :: qv3d, &\r
143 qc3d, &\r
144 qi3d, &\r
145 p3d, &\r
146 pi3d, &\r
147 t3d, &\r
148 dz8w, &\r
149 u3d, &\r
150 v3d, &\r
151 rho3d\r
152 \r
153 real, dimension( ims:ime, kms:kme, jms:jme ) , &\r
154 optional , &\r
155 intent(in ) :: qr3d, &\r
156 qs3d, &\r
157 qg3d\r
158 \r
159 real, dimension( ims:ime, kms:kme, jms:jme ) , &\r
160 intent(in) :: rthraten\r
161 \r
162 real, dimension( ims:ime, kms:kme, jms:jme ) , &\r
163 intent(inout) :: rublten, &\r
164 rvblten, &\r
165 rthblten, &\r
166 rqvblten, &\r
167 rqcblten, &\r
168 rqiblten\r
169 !\r
170 real, dimension( ims:ime, jms:jme ) , &\r
171 intent(in ) :: xland, &\r
172 hfx, &\r
173 qfx, &\r
174 psfc, &\r
175 qsfc, &\r
176 tsk, &\r
177 ust, &\r
178 rmol,&\r
179 wspd\r
180 !\r
181 real, dimension( ims:ime, kms:kme, jms:jme ) , &\r
182 intent(inout) :: pek, pep\r
183 \r
184 real, dimension( ims:ime, kms:kme, jms:jme ) , &\r
185 intent(inout) :: pek_adv,pep_adv\r
186 \r
187 real, dimension( ims:ime, kms:kme, jms:jme ) , &\r
188 intent(out) :: exch_h,exch_m\r
189 \r
190 real, dimension( ims:ime, jms:jme ) , &\r
191 intent(out ) :: pblh\r
192 \r
193 integer, dimension( ims:ime, jms:jme ) , &\r
194 intent(out ) :: kpbl\r
195 !\r
196 !local\r
197 real, dimension( its:ite, kts:kte+1):: zi,ghti \r
198 real, dimension( its:ite, kts:kte ) :: zl \r
199 \r
200 real, dimension( its:ite, kts:kte ) :: u1, &\r
201 v1, &\r
202 q1, &\r
203 q2, &\r
204 q3, &\r
205 q4, &\r
206 q5, &\r
207 q6, &\r
208 t1, &\r
209 ghtl, &\r
210 prsl, &\r
211 rho1, &\r
212 rthraten1\r
213 real, dimension( its:ite, kts:kte ) :: pek2d,pep2d\r
214 real, dimension( its:ite, kts:kte+1 ) :: exh2d,exm2d\r
215 real, dimension( its:ite) :: psfc1,ust1,rmol1, &\r
216 xland1,qsfc1,hfx1, &\r
217 qfx1,tsk1,pblh1,wspd1\r
218 integer, dimension( its:ite) :: kpbl1\r
219 \r
220 real :: rdelt\r
221 integer :: i,j,k,zz,im,kx,pp\r
222 \r
223 im=ite-its+1\r
224 kx=kte-kts+1\r
225 rdelt = 1./dt\r
226 \r
227 if(itimestep .eq. 1) then\r
228 pek_adv = pek\r
229 pep_adv = pep\r
230 end if\r
231 \r
232 do j=jts,jte\r
233 ! --------------- compute zi and zl -----------------------------------------\r
234 do i=its,ite\r
235 zi(i,kts)=0.0\r
236 enddo\r
237 !\r
238 do k=kts,kte\r
239 do i=its,ite\r
240 zi(i,k+1)=zi(i,k)+dz8w(i,k,j)\r
241 enddo\r
242 enddo\r
243 !\r
244 do k=kts,kte\r
245 do i=its,ite\r
246 zl(i,k)=0.5*(zi(i,k)+zi(i,k+1))\r
247 enddo\r
248 enddo\r
249 \r
250 ! reverse the vertical levels\r
251 pp = 0\r
252 do k=kts,kte\r
253 zz = kte-pp\r
254 do i=its,ite\r
255 u1(i,zz)=u3d(i,k,j)\r
256 v1(i,zz)=v3d(i,k,j)\r
257 t1(i,zz)=t3d(i,k,j)\r
258 q1(i,zz)=qv3d(i,k,j)/(1.+qv3d(i,k,j))\r
259 q2(i,zz)=qc3d(i,k,j)/(1.+qc3d(i,k,j))\r
260 q3(i,zz)=qi3d(i,k,j)/(1.+qi3d(i,k,j))\r
261 if(present(qr3d)) then\r
262 q4(i,zz)=qr3d(i,k,j)/(1.+qr3d(i,k,j))\r
263 else\r
264 q4(i,zz)=0.\r
265 end if\r
266 if(present(qs3d)) then\r
267 q5(i,zz)=qs3d(i,k,j)/(1.+qs3d(i,k,j))\r
268 else\r
269 q5(i,zz)=0.\r
270 end if\r
271 if(present(qg3d)) then\r
272 q6(i,zz)=qg3d(i,k,j)/(1.+qg3d(i,k,j))\r
273 else\r
274 q6(i,zz)=0.\r
275 end if\r
276 ghtl(i,zz)=zl(i,k)\r
277 prsl(i,zz) = p3d(i,k,j)\r
278 rho1(i,zz)=rho3d(i,k,j)\r
279 rthraten1(i,zz)=rthraten(i,k,j)\r
280 if(bl_eeps_tkeadv==1)then\r
281 pek2d(i,zz)= min(maxpek,max(minpek,pek_adv(i,k,j)))\r
282 pep2d(i,zz)= min(maxpep,max(minpep,pep_adv(i,k,j)))\r
283 else\r
284 pek2d(i,zz)= min(maxpek,max(minpek,pek(i,k,j)))\r
285 pep2d(i,zz)= min(maxpep,max(minpep,pep(i,k,j)))\r
286 end if\r
287 end do\r
288 pp = pp + 1\r
289 end do\r
290 \r
291 pp = 0\r
292 do k=kts,kte+1\r
293 zz = kte+1-pp\r
294 do i=its,ite\r
295 ghti(i,zz) = zi(i,k)\r
296 enddo\r
297 pp = pp + 1\r
298 enddo\r
299 \r
300 do i=its,ite\r
301 psfc1(i) = psfc(i,j)\r
302 ust1(i) = ust(i,j)\r
303 rmol1(i) = rmol(i,j)\r
304 wspd1(i) = wspd(i,j)\r
305 xland1(i)= xland(i,j)\r
306 qsfc1(i) = qsfc(i,j)\r
307 hfx1(i) = hfx(i,j)\r
308 qfx1(i) = qfx(i,j)\r
309 tsk1(i) = tsk(i,j)\r
310 end do\r
311 ! \r
312 call eeps2d(u1,v1,t1,q1,q2,q3,q4,q5,q6,prsl,ghtl,ghti,rho1,rthraten1,pek2d,pep2d &\r
313 ,psfc1,ust1,rmol1,wspd1,xland1,qsfc1,hfx1,qfx1,tsk1 &\r
314 ,dt,dx,itimestep,exh2d,exm2d,pblh1,kpbl1,im,kx)\r
315 !\r
316 pp = 0\r
317 do k=kts,kte\r
318 zz = kte-pp\r
319 do i=its,ite\r
320 rthblten(i,k,j)=(t1(i,zz)-t3d(i,k,j))/pi3d(i,k,j)*rdelt\r
321 rqvblten(i,k,j)=(q1(i,zz)/(1.-q1(i,zz))-qv3d(i,k,j))*rdelt\r
322 rqcblten(i,k,j)=(q2(i,zz)/(1.-q2(i,zz))-qc3d(i,k,j))*rdelt\r
323 rqiblten(i,k,j)=(q3(i,zz)/(1.-q3(i,zz))-qi3d(i,k,j))*rdelt\r
324 rublten(i,k,j) =(u1(i,zz)-u3d(i,k,j))*rdelt\r
325 rvblten(i,k,j) =(v1(i,zz)-v3d(i,k,j))*rdelt\r
326 \r
327 pek(i,k,j) = pek2d(i,zz)\r
328 pep(i,k,j) = pep2d(i,zz)\r
329 enddo\r
330 pp = pp + 1\r
331 enddo\r
332 \r
333 pp = 0\r
334 do k=kts,kte+1\r
335 zz = kte+1-pp\r
336 do i=its,ite\r
337 exch_h(i,k,j) = exh2d(i,zz)\r
338 exch_m(i,k,j) = exm2d(i,zz)\r
339 enddo\r
340 pp = pp + 1\r
341 enddo\r
342 \r
343 do i=its,ite\r
344 pblh(i,j) = pblh1(i)\r
345 kpbl(i,j) = kx-kpbl1(i)+1\r
346 end do\r
347 \r
348 if(bl_eeps_tkeadv==1)then\r
349 do k=kts,kte\r
350 do i=its,ite\r
351 pek_adv(i,k,j) = pek(i,k,j)\r
352 pep_adv(i,k,j) = pep(i,k,j)\r
353 end do\r
354 end do\r
355 end if\r
356 !\r
357 enddo\r
358 \r
359 !\r
360 end subroutine eeps\r
361 !\r
362 !-------------------------------------------------------------------------------\r
363 !\r
364 subroutine eeps2d(pu,pv,tz,pqv,pqc,pqi,pqr,pqs,pqg,prs,poz,zz,rho,rthraten,pek,pep, &\r
365 psfcpa,ust,rmol,wspd,xland,qsfc,hfx,qfx,tsk, &\r
366 dt,dx,itimestep,k_h,k_m,pblh,kpbl,lq,km)\r
367 !-------------------------------------------------------------------------------\r
368 implicit none\r
369 !-------------------------------------------------------------------------------\r
370 !\r
371 integer, intent(in ) :: itimestep,lq,km\r
372 real, intent(in ) :: dt, dx\r
373 !\r
374 real, dimension( lq ) , &\r
375 intent(in) :: ust, &\r
376 xland, &\r
377 hfx, &\r
378 qfx, &\r
379 rmol,&\r
380 wspd,&\r
381 psfcpa,qsfc,&\r
382 tsk\r
383 !\r
384 real, dimension( lq,km ),intent(inout) :: pu, &\r
385 pv, &\r
386 tz, &\r
387 pqv, &\r
388 pqc, &\r
389 pqi, &\r
390 pek, &\r
391 pep, &\r
392 rho, &\r
393 rthraten\r
394 \r
395 real, dimension( lq,km ),intent(in) :: prs, &\r
396 poz, &\r
397 pqr, &\r
398 pqs, &\r
399 pqg\r
400 real, dimension( lq,km+1 ),intent(in) :: zz\r
401 \r
402 real, dimension( lq,km+1 ),intent(out) :: k_h,k_m\r
403 real, dimension( lq ),intent(out) :: pblh \r
404 integer, dimension( lq ),intent(out) :: kpbl\r
405 \r
406 ! Local Vars\r
407 real, dimension( lq,km ) :: fac,psp,pt,ztv,szv,disht,up2,cpm\r
408 real, dimension( lq,km ) :: qvs,st,sh,rich,richb,dum2,dzz,doz,st0,rich0\r
409 real, dimension( lq,km ) :: am,src,src1,ax,bx,cx,ac,ab,ba,alt\r
410 real, dimension( lq,km-1 ) :: ax1,bx1,cx1,yy1,amt\r
411 real, dimension( lq,km+1 ) :: akm,akh,pnt\r
412 \r
413 real, dimension( lq) :: tkm,zkm,qkm,us2,us,dens,wstar2,sflux\r
414 real, dimension( lq) :: rm2,sfa,hs,hq\r
415 real, dimension( lq) :: gamat,gamaq\r
416 real, dimension( lq,km) :: sgk1,sgk2,edrain \r
417 integer :: j,k,mdt,lvl\r
418 real :: dtt,dum,dum0,alh,eqp,qvsw,qvsi,faf,rdz,du,dv,dss,bvf,rbdts\r
419 real :: tk,qvk,alk,af,ag,aks,aco\r
420 real :: duj,alv0,richf,rch,alv,ak\r
421 real :: govrth,bfx0,wstar3,coef\r
422 \r
423 ! Local Vars for top-down diffusion\r
424 integer, dimension( lq) :: kminrad,kcld,k700,k950\r
425 logical, dimension( lq) :: flg,flg1,scuflg\r
426 real, dimension( lq) :: zl1,wm3,ent_eff,radsum,minrad,zminrad\r
427 real, dimension( lq) :: eis\r
428 real, dimension( lq,km) :: pci,TKEprodTD,zfac,zfacent\r
429 real :: dthvx,tmp1,temps,templ\r
430 integer :: ktop,kk1,kk2\r
431 real :: radflux,rcldb,rvls,he0,he1,he2\r
432 !end top down\r
433 \r
434 mdt = 1 + int(sqrt(dt/9.0)+0.00001)\r
435 dtt = dt/mdt\r
436 \r
437 ! dzz for dz between full levels, doz for dz between half levels\r
438 do k=1,km\r
439 do j=1,lq\r
440 dzz(j,k) = zz(j,k)-zz(j,k+1)\r
441 alt(j,k) = dt/dzz(j,k)\r
442 up2(j,k)=grv/max(us2min,pu(j,k)*pu(j,k)+pv(j,k)*pv(j,k))\r
443 TKEprodTD(j,k) = 0.0\r
444 zfac(j,k) = 0.0\r
445 zfacent(j,k) = 0.0\r
446 enddo\r
447 enddo\r
448 \r
449 do k=1,km-1\r
450 do j=1,lq\r
451 doz(j,k) = poz(j,k)-poz(j,k+1)\r
452 end do\r
453 end do\r
454 !\r
455 do j=1,lq\r
456 doz(j,km)= 2.0*poz(j,km) ! not used\r
457 enddo\r
458 !\r
459 do k=1,km-1\r
460 do j=1,lq\r
461 amt(j,k) = dtt/doz(j,k)\r
462 enddo\r
463 enddo\r
464 \r
465 do k=1,km\r
466 do j=1,lq\r
467 sgk1(j,k) = 0.1\r
468 sgk2(j,k) = 0.01\r
469 end do\r
470 end do\r
471 !\r
472 do k=1,km\r
473 do j=1,lq\r
474 psp(j,k)=(prs(j,k)/p1000mb)**rcp\r
475 pt(j,k)=tz(j,k)/psp(j,k)\r
476 cpm(j,k) = cp*(1.+0.8*pqv(j,k))\r
477 dum=1.0+ep1*pqv(j,k)\r
478 ztv(j,k)=tz(j,k)*dum\r
479 szv(j,k)=pt(j,k)*dum\r
480 pci(j,k) = pqc(j,k) + pqi(j,k)\r
481 fac(j,k) = pqc(j,k) + pqi(j,k) + pqr(j,k) + pqs(j,k) + pqg(j,k)\r
482 cx(j,k)=1.0+pqv(j,k)+fac(j,k)\r
483 \r
484 alh=3.1484e6-2.37e3*tz(j,k)\r
485 eqp=611.2*exp(cs2*(tz(j,k)-cs3)/(tz(j,k)-cs4))\r
486 eqp=min(0.5*prs(j,k),eqp)\r
487 qvsw=0.622*eqp/(prs(j,k)-eqp)\r
488 eqp=611.0*exp(ci2*(tz(j,k)-ci3)/(tz(j,k)-ci4))\r
489 eqp=min(0.5*prs(j,k),eqp)\r
490 qvsi=0.622*eqp/(prs(j,k)-eqp)\r
491 if(tz(j,k).ge.t0) then\r
492 faf=0.0\r
493 else if(tz(j,k).le.hgfr) then\r
494 faf=1.0\r
495 else\r
496 if(pqi(j,k).le.epsl) then\r
497 faf=0.0\r
498 else if(pqc(j,k).le.epsl) then\r
499 faf=1.0\r
500 else\r
501 faf=pqi(j,k)/(pqc(j,k)+pqi(j,k))\r
502 endif\r
503 endif\r
504 qvs(j,k)=(1.0-faf)*qvsw+faf*qvsi\r
505 ax(j,k)=(1.0-faf)*alh+faf*xls\r
506 end do\r
507 end do\r
508 \r
509 !\r
510 ! to calculate the Richardson number and vertical shear of\r
511 ! the horizontal wind above the surface layer\r
512 !\r
513 do k=1,km-1\r
514 do j=1,lq\r
515 rdz=1.0/(poz(j,k)-poz(j,k+1))\r
516 du=pu(j,k)-pu(j,k+1)\r
517 dv=pv(j,k)-pv(j,k+1)\r
518 sh(j,k)=(du*du+dv*dv)*(rdz*rdz)\r
519 if((pqv(j,k).lt.qvs(j,k)).or.(pqv(j,k+1).lt.qvs(j,k+1))) then\r
520 dss=log(szv(j,k)/szv(j,k+1))\r
521 st(j,k)=grv*dss*rdz\r
522 else\r
523 dss=log(pt(j,k)/pt(j,k+1))\r
524 tk=0.5*(tz(j,k+1)+tz(j,k))\r
525 qvk=0.5*(pqv(j,k+1)+pqv(j,k))\r
526 alk=0.5*(ax(j,k+1)+ax(j,k))\r
527 af=alk*qvk/(rgas*tk)\r
528 ag=alk/(cp*tk)\r
529 aks=ep2*af*ag\r
530 aco=(1.0+af)/(1.0+aks)\r
531 st(j,k)=(aco*(dss+ag*(pqv(j,k)-pqv(j,k+1))) &\r
532 -(cx(j,k)-cx(j,k+1)))*grv*rdz\r
533 endif\r
534 rich(j,k)=st(j,k)/max(1.0e-7,sh(j,k))\r
535 end do\r
536 end do\r
537 \r
538 ! calculate first guess pblh\r
539 do k=1,km\r
540 do j=1,lq\r
541 richb(j,k)=poz(j,k)*(szv(j,k)/szv(j,km)-1.0)*up2(j,k)\r
542 enddo\r
543 enddo\r
544 !\r
545 do j=1,lq\r
546 pblh(j)=poz(j,km)\r
547 kpbl(j)=km\r
548 lvl=km\r
549 do k=km-1,5,-1\r
550 if(richb(j,k).ge.ricr) then\r
551 lvl=k\r
552 exit\r
553 endif\r
554 enddo\r
555 \r
556 if(lvl.lt.km) then\r
557 dum=(poz(j,lvl)-poz(j,lvl+1))/(richb(j,lvl)-richb(j,lvl+1))\r
558 pblh(j)=max(poz(j,lvl+1),poz(j,lvl)-(richb(j,lvl)-ricr)*dum)\r
559 pblh(j)=min(pblh(j),poz(j,lvl))\r
560 kpbl(j)=lvl\r
561 endif\r
562 end do\r
563 \r
564 ! preparing for tke/ep calculations\r
565 !\r
566 do j=1,lq\r
567 tkm(j) = tz(j,km)\r
568 zkm(j) = poz(j,km)\r
569 qkm(j) = pqv(j,km)\r
570 sfa(j) = tsk(j)*(p1000mb/psfcpa(j))**rcp\r
571 dens(j) = psfcpa(j)/(rgas*tkm(j)*(1.+ep1*qkm(j)))\r
572 \r
573 sh(j,km)=0.0\r
574 st(j,km)=0.0\r
575 rm2(j)=ust(j)*ust(j)/wspd(j)\r
576 hs(j)=hfx(j)/(dens(j)*cpm(j,km))\r
577 hq(j)=qfx(j)/dens(j)\r
578 !\r
579 duj=ust(j)*ust(j)\r
580 pep(j,km)=duj*ust(j)/(akv*zkm(j))\r
581 pep(j,km) = min(maxpep,max(minpep,pep(j,km)))\r
582 govrth = grv/pt(j,km)\r
583 sflux(j) = hs(j) + hq(j)*ep1*pt(j,km)\r
584 bfx0 = max(sflux(j),0.)\r
585 wstar3 = govrth*bfx0*pblh(j)\r
586 wstar2(j) = (wstar3**2)**t13\r
587 \r
588 rich(j,km) = grv*zkm(j)* &\r
589 log(szv(j,km)/(sfa(j)*(1.+ep1*qsfc(j))))/(wspd(j)**2.)\r
590 \r
591 if(rmol(j).lt.0.0 .and. sflux(j).gt.0.0) then\r
592 pek(j,km)=3.75*duj+0.2*wstar2(j)+duj*((-zkm(j)*rmol(j))**2.)**t13\r
593 else\r
594 pek(j,km)=3.75*duj\r
595 endif\r
596 pek(j,km) = min(maxpek,max(minpek,pek(j,km)))\r
597 \r
598 end do\r
599 \r
600 ! to recalculate pblh here\r
601 do j=1,lq\r
602 if(rmol(j).lt.0.0 .and. sflux(j).gt.0.0) then\r
603 dss=szv(j,km)+min(5.0,bt*sflux(j)/(max(0.1,sqrt(wstar2(j)))))\r
604 do k=1,km\r
605 richb(j,k)=poz(j,k)*(szv(j,k)/dss-1.0)*up2(j,k)\r
606 end do\r
607 \r
608 lvl=km\r
609 do k=km-1,5,-1\r
610 if(richb(j,k).ge.ricr) then\r
611 lvl=k\r
612 exit\r
613 end if\r
614 enddo\r
615 \r
616 if(lvl.lt.km) then\r
617 dum=(poz(j,lvl)-poz(j,lvl+1))/(richb(j,lvl)-richb(j,lvl+1))\r
618 pblh(j)=max(poz(j,lvl+1),poz(j,lvl)-(richb(j,lvl)-ricr)*dum)\r
619 pblh(j)=min(pblh(j),poz(j,lvl))\r
620 kpbl(j)=lvl\r
621 end if\r
622 end if\r
623 end do\r
624 \r
625 ! to add the countergradient term\r
626 do j=1,lq\r
627 if(rmol(j).lt.0.0 .and. sflux(j).gt.0.0) then\r
628 du=max(0.1,sqrt(wstar2(j)))\r
629 gamat(j)=max(bt*hs(j)/(pblh(j)*du),0.)\r
630 gamaq(j)=max(bt*hq(j)/(pblh(j)*du),0.)\r
631 else\r
632 gamat(j)=0.0\r
633 gamaq(j)=0.0\r
634 endif\r
635 end do\r
636 \r
637 ! to add TKE source driven by stratocumulus cloud top cooling\r
638 ! a switch bl_eeps_topdown controls if it is on or off\r
639 if (bl_eeps_topdown.eq.1)then\r
640 do j=1,lq\r
641 scuflg(j)=.true.\r
642 flg(j) =.true.\r
643 flg1(j) =.true.\r
644 minrad(j) = 100.\r
645 kminrad(j) = km\r
646 kcld(j) = km\r
647 k700(j) = km\r
648 k950(j) = km\r
649 zminrad(j) = poz(j,km)\r
650 zl1(j) = poz(j,km)\r
651 end do\r
652 \r
653 ! look for stratocumulus\r
654 ! look for the highest model level we will use\r
655 do k=km,2,-1\r
656 do j=1,lq\r
657 if(flg(j) .and. (prs(j,km)-prs(j,k))>3.0e4) then\r
658 k700(j)=k\r
659 flg(j)=.false.\r
660 end if\r
661 if(flg1(j) .and. (prs(j,km)-prs(j,k))>5.e3) then\r
662 k950(j)=k\r
663 flg1(j)=.false.\r
664 end if\r
665 end do\r
666 end do\r
667 ktop = minval(k700)\r
668 \r
669 ! calculate EIS based on Marquet and Bechtold (2020)\r
670 do j=1,lq\r
671 kk1 = k700(j)\r
672 kk2 = k950(j)\r
673 he1 = cp*tz(j,kk1)*(1.+5.87*(pqv(j,kk1)+pci(j,kk1))) &\r
674 -xlv*pqc(j,kk1)-xls*pqi(j,kk1)+grv*poz(j,kk1)\r
675 he2 = cp*tz(j,kk2)*(1.+5.87*(pqv(j,kk2)+pci(j,kk2))) &\r
676 -xlv*pqc(j,kk2)-xls*pqi(j,kk2)+grv*poz(j,kk2)\r
677 he0 = cp*tz(j,km)*(1.+5.87*(pqv(j,km)+pci(j,km))) &\r
678 -xlv*pqc(j,km)-xls*pqi(j,km)+grv*poz(j,km) \r
679 eis(j) = max(he1-he2,he2-he0)\r
680 eis(j) = eis(j)/cp\r
681 if(eis(j) < 5) scuflg(j)=.false.\r
682 end do\r
683 \r
684 ! look for the highest cloud level\r
685 do j=1,lq\r
686 flg(j)=scuflg(j)\r
687 end do\r
688 \r
689 do k=ktop,km\r
690 do j=1,lq\r
691 if(flg(j) .and. k >= k700(j) .and. pci(j,k) >= 1.e-6) then\r
692 kcld(j)=k\r
693 flg(j)=.false.\r
694 end if\r
695 end do\r
696 end do\r
697 \r
698 ! find the level of the minimum radiative cooling\r
699 do k=ktop,km\r
700 do j=1,lq\r
701 if(scuflg(j)) then\r
702 if(k >= kcld(j) .and. pci(j,k) >= 1.e-6) then\r
703 if(rthraten(j,k) < minrad(j)) then\r
704 minrad(j)=rthraten(j,k)\r
705 kminrad(j)=k\r
706 zminrad(j)=poz(j,k)\r
707 endif\r
708 endif\r
709 endif\r
710 end do\r
711 end do\r
712 \r
713 ! to exclude some scenarios\r
714 ! 1) cloud top should be below PBL height\r
715 ! 2) minimum radiative cooling level can't too close to surface\r
716 ! 3) minimum radiative cooling should be negative\r
717 do j=1,lq\r
718 if(scuflg(j) .and. minrad(j)>=0.) scuflg(j)=.false.\r
719 if(scuflg(j) .and. kcld(j)<kpbl(j)) scuflg(j)=.false.\r
720 if(scuflg(j) .and. kminrad(j)>=(km-1)) scuflg(j)=.false.\r
721 end do\r
722 \r
723 do j=1,lq\r
724 if (scuflg(j)) then\r
725 k = kminrad(j)\r
726 templ=tz(j,k)\r
727 !rvls is ws at full level\r
728 rvls=100.*6.112*exp(17.67*(templ-273.16)/(templ-29.65))*(ep2/prs(j,k))\r
729 temps=templ + (pqv(j,k)+pci(j,k)-rvls)/(cp/xlv + ep2*xlv*rvls/(rgas*templ**2))\r
730 rvls=100.*6.112*exp(17.67*(temps-273.15)/(temps-29.65))*(ep2/prs(j,k))\r
731 rcldb=max(pqv(j,k)+pci(j,k)-rvls,0.)\r
732 \r
733 !entrainment efficiency\r
734 dthvx = szv(j,k-2)-szv(j,k)\r
735 dthvx = max(dthvx,0.1)\r
736 tmp1 = xlv / cp * rcldb/(psp(j,k)*dthvx)\r
737 !Originally from Nichols and Turton (1986), where a2 = 60, but lowered\r
738 !here to 8, as in Grenier and Bretherton (2001).\r
739 ent_eff(j) = 0.2 + 0.2*8.*tmp1\r
740 end if\r
741 end do\r
742 \r
743 do k=km,ktop,-1\r
744 do j=1,lq\r
745 if(scuflg(j) .and. k<=(kminrad(j)+1) .and. k>=kcld(j)) then\r
746 radflux=rthraten(j,k)*psp(j,k) !converts theta/s to temp/s\r
747 radflux=radflux*cp/grv*(prs(j,k)-prs(j,k-1)) ! converts temp/s to W/m^2\r
748 if (radflux < 0.0 ) radsum(j)=abs(radflux)+radsum(j)\r
749 end if\r
750 end do\r
751 end do\r
752 \r
753 do j=1,lq\r
754 if(scuflg(j)) then\r
755 k = kminrad(j)\r
756 radsum(j)=min(radsum(j),120.0)\r
757 bfx0 = max(radsum(j)/rho(j,k)/cp,0.)\r
758 wm3(j) = grv/szv(j,k)*bfx0*min(pblh(j),1500.) ! this is wstar3(i)\r
759 end if\r
760 end do\r
761 \r
762 do k=km,ktop,-1\r
763 do j=1,lq\r
764 if(scuflg(j) .and. k>=kcld(j)) then\r
765 !Analytic vertical profile\r
766 zfac(j,k) = min(max((1.-(poz(j,k)-zl1(j))/(zminrad(j)-zl1(j))),zfmin),1.)\r
767 zfacent(j,k) = 10.*max((zminrad(j)-poz(j,k))/zminrad(j),0.0)*(1.-zfac(j,k))**3\r
768 !Calculate TKE production = 2(g/TH)(w'TH'), where w'TH' = A(TH/g)wstar^3/PBLH,\r
769 !A = ent_eff, and wstar is associated with the radiative cooling at top of PBL.\r
770 !An analytic profile controls the magnitude of this TKE prod in the vertical. \r
771 TKEprodTD(j,k)=2.*ent_eff(j)*wm3(j)/max(pblh(j),100.)*zfacent(j,k)\r
772 TKEprodTD(j,k)= max(TKEprodTD(j,k),0.0)\r
773 end if !end cloud check\r
774 end do ! end lq loop\r
775 end do ! end km loop\r
776 \r
777 end if !end top-down check\r
778 !\r
779 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
780 ! time integration of tke and ep equations\r
781 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
782 do k=1,km\r
783 do j=1,lq\r
784 richf=0.1620092\r
785 if(rich(j,k).lt.0.145) then\r
786 richf=0.6588*(rich(j,k)+0.1776-sqrt(rich(j,k)*rich(j,k)- &\r
787 0.3221*rich(j,k)+0.03156))\r
788 richf=max(-1000.0,min(0.1620092,richf))\r
789 endif\r
790 dum2(j,k)=1.116893\r
791 if(richf.lt.0.1620092) dum2(j,k)=1.318*(0.2231-richf) &\r
792 /(0.2341-richf)\r
793 end do\r
794 end do\r
795 !\r
796 lvl=1\r
797 do while(lvl .le. mdt)\r
798 \r
799 do k=1,km-1\r
800 do j=1,lq\r
801 akm(j,k)=c2*pek(j,k)*pek(j,k)/pep(j,k)\r
802 akh(j,k)=dum2(j,k)*akm(j,k)\r
803 akm(j,k)=sgk1(j,k)+akm(j,k)\r
804 akh(j,k)=sgk2(j,k)+akh(j,k)\r
805 akm(j,k)=min(akvmx,akm(j,k))\r
806 akh(j,k)=min(akvmx,akh(j,k))\r
807 if(st(j,k).gt.0.0) then\r
808 bvf=sqrt(st(j,k))\r
809 rbdts=min(1.,sqrt(rich(j,k)/0.145))\r
810 edrain(j,k)=0.1/c3*rbdts*bvf*pek(j,k)\r
811 else\r
812 edrain(j,k)=0.0 \r
813 endif\r
814 if(poz(j,k).gt.pblh(j))then\r
815 src(j,k)=akm(j,k)*sh(j,k)-akh(j,k)*st(j,k)+TKEprodTD(j,k)\r
816 src1(j,k)=akm(j,k)*sh(j,k)-akh(j,k)*st(j,k)+edrain(j,k)+TKEprodTD(j,k)\r
817 else\r
818 src(j,k)=akm(j,k)*sh(j,k)-akh(j,k)*(st(j,k) &\r
819 -2.0*grv*gamat(j)/(szv(j,k+1)+szv(j,k)))+TKEprodTD(j,k)\r
820 src1(j,k)=akm(j,k)*sh(j,k)-akh(j,k)*(st(j,k) &\r
821 -2.0*grv*gamat(j)/(szv(j,k+1)+szv(j,k)))+edrain(j,k)+TKEprodTD(j,k)\r
822 end if\r
823 end do\r
824 end do\r
825 !\r
826 do j=1,lq\r
827 akm(j,km)=c2*pek(j,km)*pek(j,km)/pep(j,km)\r
828 akh(j,km)=dum2(j,km)*akm(j,km)\r
829 end do\r
830 !\r
831 ! to calculate the vertical diffusion coefficients at u,v,levels\r
832 !\r
833 do k=2,km\r
834 do j=1,lq\r
835 am(j,k)=0.5*(akm(j,k-1)+akm(j,k))/dzz(j,k)\r
836 end do\r
837 end do\r
838 \r
839 do j=1,lq\r
840 am(j,1)=0.5*akm(j,1)/dzz(j,1)\r
841 end do\r
842 !\r
843 ! to calculate the cofficients for the triangle matrix\r
844 ! for tke dissipation and solve the matrix\r
845 !\r
846 do k=1,km-1\r
847 do j=1,lq\r
848 dum=(c3*src1(j,k)-c4*pep(j,k))*dtt/pek(j,k)\r
849 ax1(j,k)=-c5*amt(j,k)*am(j,k+1)\r
850 cx1(j,k)=-c5*amt(j,k)*am(j,k)\r
851 if(dum.le.0.5) then\r
852 bx1(j,k)=1.0+c5*amt(j,k)*(am(j,k)+am(j,k+1))-dum\r
853 yy1(j,k)=pep(j,k)\r
854 else\r
855 bx1(j,k)=1.0+c5*amt(j,k)*(am(j,k)+am(j,k+1))\r
856 yy1(j,k)=pep(j,k)+dum*pep(j,k)\r
857 endif\r
858 end do\r
859 end do\r
860 !\r
861 do j=1,lq\r
862 cx1(j,1)=0.0\r
863 ax1(j,km-1)=0.0\r
864 yy1(j,km-1)=yy1(j,km-1)+c5*amt(j,km-1)*am(j,km)*pep(j,km)\r
865 end do\r
866 !\r
867 call tridiag(ax1,bx1,cx1,yy1,lq,km-1)\r
868 !\r
869 ! to calculate the cofficients for the triangle matrix\r
870 ! for tke and solve the matrix\r
871 !\r
872 do k=1,km-1\r
873 do j=1,lq\r
874 pep(j,k)=min(maxpep,max(minpep, yy1(j,k)))\r
875 dum=(src(j,k)-pep(j,k))*dtt\r
876 ax1(j,k)=-c1*amt(j,k)*am(j,k+1)\r
877 bx1(j,k)=1.0+c1*amt(j,k)*(am(j,k)+am(j,k+1))\r
878 cx1(j,k)=-c1*amt(j,k)*am(j,k)\r
879 yy1(j,k)=pek(j,k)+dum\r
880 end do\r
881 end do\r
882 \r
883 do j=1,lq\r
884 cx1(j,1)=0.0\r
885 ax1(j,km-1)=0.0\r
886 yy1(j,km-1)=yy1(j,km-1)+c1*amt(j,km-1)*am(j,km)*pek(j,km)\r
887 end do\r
888 !\r
889 call tridiag(ax1,bx1,cx1,yy1,lq,km-1)\r
890 !\r
891 do k=1,km-1\r
892 do j=1,lq\r
893 pek(j,k)=min(maxpek,max(minpek, yy1(j,k)))\r
894 end do\r
895 end do\r
896 \r
897 lvl = lvl + 1\r
898 end do ! end while\r
899 !\r
900 !-------------------------------------------------------------\r
901 ! to calculate the vertical diffusion coeficient for momentum \r
902 ! heat, moisture\r
903 !-------------------------------------------------------------\r
904 ! note: akm and akh at the first level (km) will not be used \r
905 do k=2,km\r
906 do j=1,lq\r
907 akm(j,k)=c2*pek(j,k-1)*pek(j,k-1)/pep(j,k-1)\r
908 akh(j,k)=dum2(j,k)*akm(j,k)\r
909 akm(j,k)=sgk1(j,k)+akm(j,k)\r
910 akh(j,k)=sgk2(j,k)+akh(j,k)\r
911 akm(j,k)=min(akvmx,akm(j,k))\r
912 akh(j,k)=min(akvmx,akh(j,k))\r
913 pnt(j,k)=akh(j,k)\r
914 if(poz(j,k).gt.pblh(j)) pnt(j,k)=0.\r
915 akm(j,k)=akm(j,k)/doz(j,k-1)\r
916 akh(j,k)=akh(j,k)/doz(j,k-1)\r
917 end do\r
918 end do\r
919 !\r
920 do j=1,lq\r
921 akm(j,1)=0.0\r
922 akh(j,1)=0.0\r
923 pnt(j,1)=0.0\r
924 akm(j,km+1)=0.0\r
925 akh(j,km+1)=0.0\r
926 pnt(j,km+1)=0.0\r
927 end do\r
928 !\r
929 ! Momentum fluxes calculation\r
930 !\r
931 do k=1,km\r
932 do j=1,lq\r
933 ax(j,k)=-alt(j,k)*akm(j,k+1)\r
934 bx(j,k)=1.0+alt(j,k)*(akm(j,k)+akm(j,k+1))\r
935 cx(j,k)=-alt(j,k)*akm(j,k)\r
936 ba(j,k)=ax(j,k)\r
937 ab(j,k)=bx(j,k)\r
938 ac(j,k)=cx(j,k)\r
939 end do\r
940 end do\r
941 \r
942 do j=1,lq\r
943 bx(j,km) = bx(j,km)+alt(j,km)*rm2(j)\r
944 ab(j,km)=bx(j,km)\r
945 end do\r
946 \r
947 !\r
948 call tridiag(ba,ab,ac,pu,lq,km)\r
949 !\r
950 call tridiag(ax,bx,cx,pv,lq,km)\r
951 !\r
952 ! Heat and moisture fluxes calculation\r
953 !\r
954 do k=1,km\r
955 do j=1,lq\r
956 ax(j,k)=-alt(j,k)*akh(j,k+1)\r
957 bx(j,k)=1.0+alt(j,k)*(akh(j,k)+akh(j,k+1))\r
958 cx(j,k)=-alt(j,k)*akh(j,k)\r
959 ba(j,k)=ax(j,k)\r
960 ab(j,k)=bx(j,k)\r
961 ac(j,k)=cx(j,k)\r
962 dum0=0.0 \r
963 if(k.gt.1) dum0=pep(j,k-1)\r
964 dum=pep(j,k)\r
965 if(dum0.le.1.e-6) dum0=0.0\r
966 if(dum.le.1.e-6) dum=0.0\r
967 disht(j,k)=0.5*(dum0+dum)/cpm(j,k)\r
968 pt(j,k)=pt(j,k)+disht(j,k)*dt/psp(j,k)+alt(j,k)* &\r
969 gamat(j)*(pnt(j,k+1)-pnt(j,k))\r
970 end do\r
971 end do\r
972 !\r
973 do j=1,lq\r
974 pt(j,km)=pt(j,km)+alt(j,km)*hs(j)\r
975 end do\r
976 !\r
977 call tridiag(ba,ab,ac,pt,lq,km)\r
978 !\r
979 do k=1,km\r
980 do j=1,lq\r
981 ba(j,k)=ax(j,k)\r
982 ab(j,k)=bx(j,k)\r
983 ac(j,k)=cx(j,k)\r
984 end do\r
985 end do\r
986 \r
987 do k=1,km\r
988 do j=1,lq\r
989 pqv(j,k)=pqv(j,k)+alt(j,k)*gamaq(j)* &\r
990 (pnt(j,k+1)-pnt(j,k))\r
991 end do\r
992 end do\r
993 \r
994 do j=1,lq\r
995 pqv(j,km)=pqv(j,km)+alt(j,km)*hq(j)\r
996 end do\r
997 !\r
998 call tridiag(ba,ab,ac,pqv,lq,km)\r
999 !\r
1000 do k=1,km\r
1001 do j=1,lq\r
1002 ba(j,k)=ax(j,k)\r
1003 ab(j,k)=bx(j,k)\r
1004 ac(j,k)=cx(j,k)\r
1005 end do\r
1006 end do\r
1007 \r
1008 call tridiag(ba,ab,ac,pqc,lq,km)\r
1009 !\r
1010 call tridiag(ax,bx,cx,pqi,lq,km)\r
1011 !\r
1012 do k=1,km\r
1013 do j=1,lq\r
1014 tz(j,k)=pt(j,k)*psp(j,k)\r
1015 pqv(j,k)=max(epsl,pqv(j,k))\r
1016 pqc(j,k)=max(epsl,pqc(j,k))\r
1017 pqi(j,k)=max(epsl,pqi(j,k))\r
1018 !-----------------------------------------------------\r
1019 ! immediate melting of cloud ice below freezing level\r
1020 !-----------------------------------------------------\r
1021 if(tz(j,k).ge.t0.and.pqi(j,k).gt.epsl) then\r
1022 if(pqc(j,k).le.epsl) then\r
1023 pqc(j,k)=pqi(j,k)\r
1024 else\r
1025 pqc(j,k)=pqc(j,k)+pqi(j,k)\r
1026 endif\r
1027 alh=3.1484e6-2.37e3*tz(j,k)\r
1028 tz(j,k)=tz(j,k)-(xls-alh)*pqi(j,k)/cpm(j,k)\r
1029 pqi(j,k)=epsl\r
1030 endif\r
1031 end do\r
1032 end do\r
1033 \r
1034 do k=2,km\r
1035 do j=1,lq\r
1036 k_h(j,k)=akh(j,k)*doz(j,k-1)\r
1037 k_m(j,k)=akm(j,k)*doz(j,k-1)\r
1038 end do\r
1039 end do\r
1040 \r
1041 do j=1,lq\r
1042 k_h(j,1)=0.\r
1043 k_m(j,1)=0.\r
1044 k_m(j,km+1)= 0. !c2*pek(j,km)*pek(j,km)/pep(j,km) ! we don't need it\r
1045 k_h(j,km+1)= 0. !dum2(j,km)*k_m(j,km+1) ! we don't need it\r
1046 end do\r
1047 \r
1048 end subroutine eeps2d\r
1049 \r
1050 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
1051 !\r
1052 ! ==================================================================\r
1053 subroutine tridiag(a,b,c,d,lq,km)\r
1054 !! to solve system of linear eqs on tridiagonal matrix n times n\r
1055 !! after Peaceman and Rachford, 1955\r
1056 !! a,b,c,d - are vectors of order n \r
1057 !! a,b,c - are coefficients on the LHS\r
1058 !! d - is initially RHS on the output becomes a solution vector\r
1059 \r
1060 !-------------------------------------------------------------------\r
1061 implicit none\r
1062 INTEGER, INTENT(in):: lq,km\r
1063 REAL, DIMENSION(lq,km), INTENT(in) :: a,b\r
1064 REAL, DIMENSION(lq,km), INTENT(inout) :: c,d\r
1065 \r
1066 INTEGER :: j,k\r
1067 REAL :: p\r
1068 REAL, DIMENSION(lq,km) :: q\r
1069 \r
1070 do j=1,lq\r
1071 c(j,1)=0.\r
1072 q(j,km)=-c(j,km)/b(j,km)\r
1073 d(j,km)=d(j,km)/b(j,km)\r
1074 end do\r
1075 \r
1076 DO k=km-1,1,-1\r
1077 do j=1,lq\r
1078 p=1./(b(j,k)+a(j,k)*q(j,k+1))\r
1079 q(j,k)=-c(j,k)*p\r
1080 d(j,k)=(d(j,k)-a(j,k)*d(j,k+1))*p\r
1081 end do\r
1082 ENDDO\r
1083 \r
1084 DO k=2,km\r
1085 do j=1,lq\r
1086 d(j,k)=d(j,k)+q(j,k)*d(j,k-1)\r
1087 end do\r
1088 ENDDO\r
1089 \r
1090 end subroutine tridiag\r
1091 \r
1092 !-------------------------------------------------------------------------------\r
1093 subroutine eepsinit(rublten,rvblten,rthblten,rqvblten, &\r
1094 rqcblten,rqiblten,p_qi,p_first_scalar,pek,pep, &\r
1095 restart, allowed_to_read, &\r
1096 ids, ide, jds, jde, kds, kde, &\r
1097 ims, ime, jms, jme, kms, kme, &\r
1098 its, ite, jts, jte, kts, kte )\r
1099 !-------------------------------------------------------------------------------\r
1100 implicit none\r
1101 !-------------------------------------------------------------------------------\r
1102 !\r
1103 logical , intent(in) :: restart, allowed_to_read\r
1104 integer , intent(in) :: ids, ide, jds, jde, kds, kde, &\r
1105 ims, ime, jms, jme, kms, kme, &\r
1106 its, ite, jts, jte, kts, kte\r
1107 integer , intent(in) :: p_qi,p_first_scalar\r
1108 real , dimension( ims:ime , kms:kme , jms:jme ), intent(out) :: &\r
1109 rublten, &\r
1110 rvblten, &\r
1111 rthblten, &\r
1112 rqvblten, &\r
1113 rqcblten, &\r
1114 rqiblten\r
1115 real, dimension( ims:ime, kms:kme, jms:jme ),intent(out) :: pek, pep\r
1116 integer :: i, j, k, itf, jtf, ktf\r
1117 !\r
1118 jtf = min0(jte,jde-1)\r
1119 ktf = min0(kte,kde-1)\r
1120 itf = min0(ite,ide-1)\r
1121 !\r
1122 if(.not.restart)then\r
1123 do j = jts,jtf\r
1124 do k = kts,ktf\r
1125 do i = its,itf\r
1126 rublten(i,k,j) = 0.\r
1127 rvblten(i,k,j) = 0.\r
1128 rthblten(i,k,j) = 0.\r
1129 rqvblten(i,k,j) = 0.\r
1130 rqcblten(i,k,j) = 0.\r
1131 pek(i,k,j) = minpek\r
1132 pep(i,k,j) = minpep\r
1133 enddo\r
1134 enddo\r
1135 enddo\r
1136 endif\r
1137 !\r
1138 if (p_qi .ge. p_first_scalar .and. .not.restart) then\r
1139 do j = jts,jtf\r
1140 do k = kts,ktf\r
1141 do i = its,itf\r
1142 rqiblten(i,k,j) = 0.\r
1143 enddo\r
1144 enddo\r
1145 enddo\r
1146 endif\r
1147 !\r
1148 end subroutine eepsinit\r
1149 !------------------------------------------------------------------------------- \r
1150 end module module_bl_eeps\r