| blob | aa99b004b0cca18f90559d38205c06fe7667dd7b |
1 !
2 !
3 !
4 !
5 MODULE module_shcu_nscv
6 CONTAINS
7 !-------------------------------------------------------------------------------
8 subroutine shcu_nscv(dt,p3di,p3d,pi3d,qc3d,qi3d,rho3d, &
9 qv3d,t3d,raincv,xland,dz8w,w,u3d,v3d, &
10 hpbl,hfx,qfx, &
11 mp_physics, &
12 pgcon, &
13 cp,cliq,cpv,g,xlv,r_d,r_v,ep_1,ep_2, &
14 cice,xls,psat,f_qi,f_qc, &
15 rthshten,rqvshten,rqcshten,rqishten, &
16 rushten,rvshten, &
17 pratesh,hbot,htop, &
18 ids,ide, jds,jde, kds,kde, &
19 ims,ime, jms,jme, kms,kme, &
20 its,ite, jts,jte, kts,kte)
21 !-------------------------------------------------------------------------------
22 implicit none
23 !-------------------------------------------------------------------------------
24 !-- dt time step (s)
25 !-- p3di 3d pressure (pa) at interface level
26 !-- p3d 3d pressure (pa)
27 !-- pi3d 3d exner function (dimensionless)
28 !-- z height above sea level (m)
29 !-- qc3d cloud water mixing ratio (kg/kg)
30 !-- qi3d cloud ice mixing ratio (kg/kg)
31 !-- qv3d 3d water vapor mixing ratio (kg/kg)
32 !-- t3d temperature (k)
33 !-- w vertical velocity (m/s)
34 !-- dz8w dz between full levels (m)
35 !-- u3d 3d u-velocity interpolated to theta points (m/s)
36 !-- v3d 3d v-velocity interpolated to theta points (m/s)
37 !-- ids start index for i in domain
38 !-- ide end index for i in domain
39 !-- jds start index for j in domain
40 !-- jde end index for j in domain
41 !-- kds start index for k in domain
42 !-- kde end index for k in domain
43 !-- ims start index for i in memory
44 !-- ime end index for i in memory
45 !-- jms start index for j in memory
46 !-- jme end index for j in memory
47 !-- kms start index for k in memory
48 !-- kme end index for k in memory
49 !-- its start index for i in tile
50 !-- ite end index for i in tile
51 !-- jts start index for j in tile
52 !-- jte end index for j in tile
53 !-- kts start index for k in tile
54 !-- kte end index for k in tile
55 !-------------------------------------------------------------------------------
56 integer, intent(in ) :: ids,ide, jds,jde, kds,kde, &
57 ims,ime, jms,jme, kms,kme, &
58 its,ite, jts,jte, kts,kte
59 real, intent(in ) :: cp,cliq,cpv,g,xlv,r_d,r_v,ep_1,ep_2, &
60 cice,xls,psat
61 real, intent(in ) :: dt
62 real, optional, intent(in ) :: pgcon
63 real, dimension( ims:ime, kms:kme, jms:jme ),optional ,&
64 intent(inout) :: rthshten,&
65 rushten,&
66 rvshten,&
67 rqcshten,&
68 rqishten,&
69 rqvshten
70 logical, optional :: F_QC,F_QI
71 real, dimension( ims:ime, kms:kme, jms:jme ) ,&
72 intent(in ) :: qv3d,&
73 qc3d,&
74 qi3d,&
75 rho3d,&
76 p3d,&
77 pi3d,&
78 t3d
79 real, dimension( ims:ime, kms:kme, jms:jme ) ,&
80 intent(in ) :: p3di
81 real, dimension( ims:ime, kms:kme, jms:jme ) ,&
82 intent(in ) :: dz8w,&
83 w
84 real, dimension( ims:ime, jms:jme ) , &
85 intent(in ) :: raincv
86 real, dimension( ims:ime, jms:jme ) ,&
87 intent(inout) :: pratesh
88 real, dimension( ims:ime, jms:jme ) ,&
89 intent(out) :: hbot,&
90 htop
91 !
92 real, dimension( ims:ime, jms:jme ) ,&
93 intent(in ) :: xland
94 !
95 real, dimension( ims:ime, kms:kme, jms:jme ) ,&
96 intent(in ) :: u3d,&
97 v3d
98 !
99 real, dimension( ims:ime, jms:jme ) ,&
100 intent(in ) :: hpbl,&
101 hfx,&
102 qfx
103 integer, intent(in ) :: mp_physics
104 integer :: ncloud
105 !
106 ! local
107 !
108 real, dimension( its:ite, kts:kte ) :: del,&
109 prsll,&
110 dot,&
111 u1,&
112 v1,&
113 t1,&
114 q1, &
115 qc2,&
116 qi2
117 real, dimension( its:ite, kts:kte+1 ) :: prsii,&
118 zii
119 real, dimension( its:ite, kts:kte ) :: zll
120 real, dimension( its:ite) :: rain
121 real :: delt,rdelt
122 integer, dimension (its:ite) :: kbot,&
123 ktop,&
124 icps
125 real :: pgcon_use
126 integer :: i,j,k,kp
127 !
128 ! microphysics scheme --> ncloud
129 !
130 if (mp_physics .eq. 0) then
131 ncloud = 0
132 elseif ( mp_physics .eq. 1 .or. mp_physics .eq. 3 ) then
133 ncloud = 1
134 else
135 ncloud = 2
136 endif
137 !
138 if(present(pgcon)) then
139 pgcon_use = pgcon
140 else
141 ! pgcon_use = 0.7 ! Gregory et al. (1997, QJRMS)
142 pgcon_use = 0.55 ! Zhang & Wu (2003,JAS)
143 ! 0.55 is a physically-based value used by GFS
144 ! HWRF uses 0.2, for model tuning purposes
145 endif
146 !
147 delt=dt
148 rdelt=1./delt
149 !
150 ! outer most J_loop
151 !
152 do j = jts,jte
153 do k = kts,kte
154 kp = k+1
155 do i = its,ite
156 dot(i,k) = -5.0e-4*g*rho3d(i,k,j)*(w(i,k,j)+w(i,kp,j))
157 prsll(i,k)=p3d(i,k,j)*0.001
158 prsii(i,k)=p3di(i,k,j)*0.001
159 enddo
160 enddo
161 !
162 do i = its,ite
163 prsii(i,kte+1)=p3di(i,kte+1,j)*0.001
164 enddo
165 !
166 do i = its,ite
167 zii(i,1)=0.0
168 enddo
169 !
170 do k = kts,kte
171 do i = its,ite
172 zii(i,k+1)=zii(i,k)+dz8w(i,k,j)
173 enddo
174 enddo
175 !
176 do k = kts,kte
177 do i = its,ite
178 zll(i,k)=0.5*(zii(i,k)+zii(i,k+1))
179 enddo
180 enddo
181 !
182 do k = kts,kte
183 do i = its,ite
184 del(i,k)=prsll(i,k)*g/r_d*dz8w(i,k,j)/t3d(i,k,j)
185 u1(i,k)=u3d(i,k,j)
186 v1(i,k)=v3d(i,k,j)
187 q1(i,k)=qv3d(i,k,j)
188 ! q1(i,k)=qv3d(i,k,j)/(1.+qv3d(i,k,j))
189 t1(i,k)=t3d(i,k,j)
190 qi2(i,k) = qi3d(i,k,j)
191 qc2(i,k) = qc3d(i,k,j)
192 enddo
193 enddo
194 !
195 icps(:) = 0
196 do i = its,ite
197 if(raincv(i,j) .gt. 1.e-30) icps(i)=1
198 enddo
199 !
200 ! NCEP SCV
201 !
202 call nscv2d(delt=delt,del=del(its,kts),prsl=prsll(its,kts), &
203 prsi=prsii(its,kts),prslk=pi3d(ims,kms,j),zl=zll(its,kts), &
204 ncloud=ncloud,qc2=qc2(its,kts),qi2=qi2(its,kts), &
205 q1=q1(its,kts),t1=t1(its,kts),rain=rain(its), &
206 kbot=kbot(its),ktop=ktop(its), &
207 icps=icps(its), &
208 slimsk=xland(ims,j),dot=dot(its,kts), &
209 u1=u1(its,kts), v1=v1(its,kts), &
210 cp_=cp,cliq_=cliq,cvap_=cpv,g_=g,hvap_=xlv, &
211 rd_=r_d,rv_=r_v,fv_=ep_1,ep2=ep_2, &
212 cice=cice,xls=xls,psat=psat, &
213 hpbl=hpbl(ims,j),hfx=hfx(ims,j),qfx=qfx(ims,j), &
214 pgcon=pgcon_use, &
215 ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde, &
216 ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme, &
217 its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte )
218 !
219 do i = its,ite
220 pratesh(i,j) = rain(i)*1000./dt
221 hbot(i,j) = kbot(i)
222 htop(i,j) = ktop(i)
223 enddo
224 !
225 IF(PRESENT(rthshten).AND.PRESENT(rqvshten)) THEN
226 do k = kts,kte
227 do i = its,ite
228 rthshten(i,k,j)=(t1(i,k)-t3d(i,k,j))/pi3d(i,k,j)*rdelt
229 rqvshten(i,k,j)=(q1(i,k)-qv3d(i,k,j))*rdelt
230 enddo
231 enddo
232 ENDIF
233 !
234 IF(PRESENT(rushten).AND.PRESENT(rvshten)) THEN
235 do k = kts,kte
236 do i = its,ite
237 rushten(i,k,j)=(u1(i,k)-u3d(i,k,j))*rdelt
238 rvshten(i,k,j)=(v1(i,k)-v3d(i,k,j))*rdelt
239 enddo
240 enddo
241 ENDIF
242 !
243 IF(PRESENT( rqishten )) THEN
244 IF ( F_QI ) THEN
245 do k = kts,kte
246 do i = its,ite
247 rqishten(i,k,j)=(qi2(i,k)-qi3d(i,k,j))*rdelt
248 enddo
249 enddo
250 ENDIF
251 ENDIF
252 !
253 IF(PRESENT( rqcshten )) THEN
254 IF ( F_QC ) THEN
255 do k = kts,kte
256 do i = its,ite
257 rqcshten(i,k,j)=(qc2(i,k)-qc3d(i,k,j))*rdelt
258 enddo
259 enddo
260 ENDIF
261 ENDIF
262 !
263 enddo ! outer most J_loop
264 !
265 return
266 end subroutine shcu_nscv
267 !-------------------------------------------------------------------------------
268 !
269 !-------------------------------------------------------------------------------
270 ! NCEP SCV (Shallow Convection Scheme)
271 !-------------------------------------------------------------------------------
272 subroutine nscv2d(delt,del,prsl,prsi,prslk,zl, &
273 ncloud,qc2,qi2,q1,t1,rain,kbot,ktop, &
274 icps, &
275 slimsk,dot,u1,v1, &
276 cp_,cliq_,cvap_,g_,hvap_,rd_,rv_,fv_,ep2, &
277 cice,xls,psat, &
278 hpbl,hfx,qfx, &
279 pgcon, &
280 ids,ide, jds,jde, kds,kde, &
281 ims,ime, jms,jme, kms,kme, &
282 its,ite, jts,jte, kts,kte)
283 !-------------------------------------------------------------------------------
284 !
285 ! subprogram: nscv2d computes shallow-convective heating and moistening
286 !
287 ! abstract: computes non-precipitating convective heating and moistening
288 ! using a one cloud type arakawa-schubert convection scheme as in the ncep
289 ! sas scheme. the scheme has been operational at ncep gfs model since july 2010
290 ! the scheme includes updraft and downdraft effects, but the cloud depth is
291 ! limited less than 150 hpa.
292 !
293 ! developed by jong-il han and hua-lu pan
294 ! implemented into wrf by jihyeon jang and songyou hong
295 ! module with cpp-based options is available in grims
296 !
297 ! program history log:
298 ! 10-07-01 jong-il han initial operational implementation at ncep gfs
299 ! 10-12-01 jihyeon jang implemented into wrf
300 ! 18-05-04 jihyeon jang seperated the shallow convection module from nsas
301 !
302 ! subprograms called:
303 ! fpvs - function to compute saturation vapor pressure
304 !
305 ! references:
306 ! han and pan (2011, wea. forecasting)
307 !
308 !-------------------------------------------------------------------------------
309 implicit none
310 !-------------------------------------------------------------------------------
311 !
312 ! in/out variables
313 !
314 integer :: ids,ide, jds,jde, kds,kde, &
315 ims,ime, jms,jme, kms,kme, &
316 its,ite, jts,jte, kts,kte
317 real :: cp_,cliq_,cvap_,g_,hvap_,rd_,rv_,fv_,ep2
318 real :: pi_,qmin_,t0c_
319 real :: cice,xlv0,xls,psat
320 !
321 real :: delt
322 real :: del(its:ite,kts:kte), &
323 prsl(its:ite,kts:kte),prslk(ims:ime,kms:kme), &
324 prsi(its:ite,kts:kte+1),zl(its:ite,kts:kte)
325 integer :: ncloud
326 real :: slimsk(ims:ime)
327 real :: dot(its:ite,kts:kte)
328 real :: hpbl(ims:ime)
329 real :: rcs
330 real :: hfx(ims:ime),qfx(ims:ime)
331 !
332 real :: qi2(its:ite,kts:kte),qc2(its:ite,kts:kte)
333 real :: q1(its:ite,kts:kte), &
334 t1(its:ite,kts:kte), &
335 u1(its:ite,kts:kte), &
336 v1(its:ite,kts:kte)
337 integer :: icps(its:ite)
338 !
339 real :: rain(its:ite)
340 integer :: kbot(its:ite),ktop(its:ite)
341 !
342 ! local variables and arrays
343 !
344 integer :: i,j,indx, jmn, k, kk, km1
345 integer :: kpbl(its:ite)
346 !
347 real :: dellat, &
348 desdt, deta, detad, dg, &
349 dh, dhh, dlnsig, dp, &
350 dq, dqsdp, dqsdt, dt, &
351 dt2, dtmax, dtmin, &
352 dv1h, dv2h, dv3h, &
353 dv1q, dv2q, dv3q, &
354 dv1u, dv2u, dv3u, &
355 dv1v, dv2v, dv3v, &
356 dz, dz1, e1, clam, &
357 aafac, &
358 es, etah, &
359 evef, evfact, evfactl, &
360 factor, fjcap, &
361 gamma, pprime, betaw, &
362 qlk, qrch, qs, &
363 rfact, shear, tem1, &
364 tem2, val, val1, &
365 val2, w1, w1l, w1s, &
366 w2, w2l, w2s, w3, &
367 w3l, w3s, w4, w4l, &
368 w4s, tem, ptem, ptem1, &
369 pgcon
370 !
371 integer :: kb(its:ite), kbcon(its:ite), kbcon1(its:ite), &
372 ktcon(its:ite), ktcon1(its:ite), &
373 kbm(its:ite), kmax(its:ite)
374 !
375 real :: aa1(its:ite), &
376 delhbar(its:ite), delq(its:ite), &
377 delq2(its:ite), delqev(its:ite), rntot(its:ite), &
378 delqbar(its:ite), deltbar(its:ite), &
379 deltv(its:ite), edt(its:ite), &
380 wstar(its:ite), sflx(its:ite), &
381 pdot(its:ite), po(its:ite,kts:kte), &
382 qcond(its:ite), qevap(its:ite), hmax(its:ite), &
383 vshear(its:ite), &
384 xlamud(its:ite), xmb(its:ite), xmbmax(its:ite)
385 real :: delubar(its:ite), delvbar(its:ite)
386 !
387 real :: cincr
388 !
389 real :: thx(its:ite, kts:kte)
390 real :: rhox(its:ite)
391 real :: tvcon
392 !
393 real :: p(its:ite,kts:kte), to(its:ite,kts:kte), &
394 qo(its:ite,kts:kte), qeso(its:ite,kts:kte), &
395 uo(its:ite,kts:kte), vo(its:ite,kts:kte)
396 !
397 ! cloud water
398 !
399 real :: qlko_ktcon(its:ite), dellal(its:ite,kts:kte), &
400 dbyo(its:ite,kts:kte), &
401 xlamue(its:ite,kts:kte), &
402 heo(its:ite,kts:kte), heso(its:ite,kts:kte), &
403 dellah(its:ite,kts:kte), dellaq(its:ite,kts:kte), &
404 dellau(its:ite,kts:kte), dellav(its:ite,kts:kte), &
405 ucko(its:ite,kts:kte), vcko(its:ite,kts:kte), &
406 hcko(its:ite,kts:kte), qcko(its:ite,kts:kte), &
407 eta(its:ite,kts:kte), zi(its:ite,kts:kte+1), &
408 pwo(its:ite,kts:kte)
409 !
410 logical :: totflg, cnvflg(its:ite), flg(its:ite)
411 !
412 ! physical parameters
413 !
414 real,parameter :: c0=.002,c1=5.e-4
415 real,parameter :: cincrmax=180.,cincrmin=120.,dthk=25.
416 real :: el2orc,fact1,fact2,eps
417 real,parameter :: h1=0.33333333
418 real,parameter :: tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf)
419 !-------------------------------------------------------------------------------
420 pi_ = 3.14159
421 qmin_ = 1.0e-30
422 t0c_ = 273.15
423 xlv0 = hvap_
424 km1 = kte - 1
425 !
426 ! compute surface buoyancy flux
427 !
428 do k = kts,kte
429 do i = its,ite
430 thx(i,k) = t1(i,k)/prslk(i,k)
431 enddo
432 enddo
433 !
434 do i = its,ite
435 tvcon = (1.+fv_*q1(i,1))
436 rhox(i) = prsl(i,1)*1.e3/(rd_*t1(i,1)*tvcon)
437 enddo
438 !
439 do i = its,ite
440 ! sflx(i) = heat(i)+fv_*t1(i,1)*evap(i)
441 sflx(i) = hfx(i)/rhox(i)/cp_ + qfx(i)/rhox(i)*fv_*thx(i,1)
442 enddo
443 !
444 ! initialize arrays
445 !
446 do i = its,ite
447 cnvflg(i) = .true.
448 if(icps(i).eq.1) cnvflg(i) = .false.
449 if(sflx(i).le.0.) cnvflg(i) = .false.
450 if(cnvflg(i)) then
451 kbot(i)=kte+1
452 ktop(i)=0
453 endif
454 rain(i)=0.
455 kbcon(i)=kte
456 ktcon(i)=1
457 kb(i)=kte
458 pdot(i) = 0.
459 qlko_ktcon(i) = 0.
460 edt(i) = 0.
461 aa1(i) = 0.
462 vshear(i) = 0.
463 enddo
464 !
465 totflg = .true.
466 do i = its,ite
467 totflg = totflg .and. (.not. cnvflg(i))
468 enddo
469 if(totflg) return
470 !
471 dt2 = delt
472 val = 1200.
473 dtmin = max(dt2, val )
474 val = 3600.
475 dtmax = max(dt2, val )
476 !
477 ! model tunable parameters are all here
478 !
479 clam = .3
480 aafac = .1
481 betaw = .03
482 evfact = 0.3
483 evfactl = 0.3
484 val = 1.
485 !
486 ! define miscellaneous values
487 !
488 el2orc = hvap_*hvap_/(rv_*cp_)
489 eps = rd_/rv_
490 fact1 = (cvap_-cliq_)/rv_
491 fact2 = hvap_/rv_-fact1*t0c_
492 !
493 w1l = -8.e-3
494 w2l = -4.e-2
495 w3l = -5.e-3
496 w4l = -5.e-4
497 w1s = -2.e-4
498 w2s = -2.e-3
499 w3s = -1.e-3
500 w4s = -2.e-5
501 !
502 ! define top layer for search of the downdraft originating layer
503 ! and the maximum thetae for updraft
504 !
505 do i = its,ite
506 kbm(i) = kte
507 kmax(i) = kte
508 enddo
509 !
510 do k = kts,kte
511 do i = its,ite
512 if (prsl(i,k).gt.prsi(i,1)*0.70) kbm(i) = k + 1
513 if (prsl(i,k).gt.prsi(i,1)*0.60) kmax(i) = k + 1
514 enddo
515 enddo
516 !
517 do i = its,ite
518 kbm(i) = min(kbm(i),kmax(i))
519 enddo
520 !
521 ! hydrostatic height assume zero terr and compute
522 ! updraft entrainment rate as an inverse function of height
523 !
524 do k = kts+1,kte
525 do i = its,ite
526 zi(i,k) = 0.5*(zl(i,k-1)+zl(i,k))
527 enddo
528 enddo
529 !
530 do k = kts,km1
531 do i = its,ite
532 xlamue(i,k) = clam / zi(i,k+1)
533 enddo
534 enddo
535 !
536 do i = its,ite
537 xlamue(i,kte) = xlamue(i,km1)
538 enddo
539 !
540 ! pbl height
541 !
542 do i = its,ite
543 flg(i) = cnvflg(i)
544 kpbl(i)= 1
545 enddo
546 !
547 do k = kts+1,km1
548 do i = its,ite
549 if (flg(i).and.zl(i,k).le.hpbl(i)) then
550 kpbl(i) = k
551 else
552 flg(i) = .false.
553 endif
554 enddo
555 enddo
556 !
557 do i = its,ite
558 kpbl(i)= min(kpbl(i),kbm(i))
559 enddo
560 !
561 ! convert surface pressure to mb from cb
562 !
563 rcs = 1.
564 do k = kts,kte
565 do i = its,ite
566 if (cnvflg(i) .and. k .le. kmax(i)) then
567 p(i,k) = prsl(i,k) * 10.0
568 eta(i,k) = 1.
569 hcko(i,k) = 0.
570 qcko(i,k) = 0.
571 ucko(i,k) = 0.
572 vcko(i,k) = 0.
573 dbyo(i,k) = 0.
574 pwo(i,k) = 0.
575 dellal(i,k) = 0.
576 to(i,k) = t1(i,k)
577 qo(i,k) = q1(i,k)
578 uo(i,k) = u1(i,k) * rcs
579 vo(i,k) = v1(i,k) * rcs
580 endif
581 enddo
582 enddo
583 !
584 !
585 ! column variables
586 ! p is pressure of the layer (mb)
587 ! t is temperature at t-dt (k)..tn
588 ! q is mixing ratio at t-dt (kg/kg)..qn
589 ! to is temperature at t+dt (k)... this is after advection and turbulan
590 ! qo is mixing ratio at t+dt (kg/kg)..q1
591 !
592 do k = kts, kte
593 do i=its,ite
594 if (cnvflg(i) .and. k .le. kmax(i)) then
595 qeso(i,k) = 0.01 * fpvs(to(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
596 qeso(i,k) = eps * qeso(i,k) / (p(i,k) + (eps-1.)*qeso(i,k))
597 val1 = 1.e-8
598 qeso(i,k) = max(qeso(i,k), val1)
599 val2 = 1.e-10
600 qo(i,k) = max(qo(i,k), val2 )
601 endif
602 enddo
603 enddo
604 !
605 ! compute moist static energy
606 !
607 do k = kts,kte
608 do i=its,ite
609 if (cnvflg(i) .and. k .le. kmax(i)) then
610 tem = g_ * zl(i,k) + cp_ * to(i,k)
611 heo(i,k) = tem + hvap_ * qo(i,k)
612 heso(i,k) = tem + hvap_ * qeso(i,k)
613 endif
614 enddo
615 enddo
616 !
617 ! determine level with largest moist static energy within pbl
618 ! this is the level where updraft starts
619 !
620 do i=its,ite
621 if (cnvflg(i)) then
622 hmax(i) = heo(i,1)
623 kb(i) = 1
624 endif
625 enddo
626 !
627 do k = kts+1, kte
628 do i=its,ite
629 if (cnvflg(i).and.k.le.kpbl(i)) then
630 if(heo(i,k).gt.hmax(i)) then
631 kb(i) = k
632 hmax(i) = heo(i,k)
633 endif
634 endif
635 enddo
636 enddo
637 !
638 do k = kts, km1
639 do i=its,ite
640 if (cnvflg(i) .and. k .le. kmax(i)-1) then
641 dz = .5 * (zl(i,k+1) - zl(i,k))
642 dp = .5 * (p(i,k+1) - p(i,k))
643 es = 0.01*fpvs(to(i,k+1),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
644 pprime = p(i,k+1) + (eps-1.) * es
645 qs = eps * es / pprime
646 dqsdp = - qs / pprime
647 desdt = es * (fact1 / to(i,k+1) + fact2 / (to(i,k+1)**2))
648 dqsdt = qs * p(i,k+1) * desdt / (es * pprime)
649 gamma = el2orc * qeso(i,k+1) / (to(i,k+1)**2)
650 dt = (g_ * dz + hvap_ * dqsdp * dp) / (cp_ * (1. + gamma))
651 dq = dqsdt * dt + dqsdp * dp
652 to(i,k) = to(i,k+1) + dt
653 qo(i,k) = qo(i,k+1) + dq
654 po(i,k) = .5 * (p(i,k) + p(i,k+1))
655 endif
656 enddo
657 enddo
658 !
659 do k = kts, km1
660 do i=its,ite
661 if (cnvflg(i) .and. k .le. kmax(i)-1) then
662 qeso(i,k)=0.01*fpvs(to(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
663 qeso(i,k) = eps * qeso(i,k) / (po(i,k) + (eps-1.) * qeso(i,k))
664 val1 = 1.e-8
665 qeso(i,k) = max(qeso(i,k), val1)
666 val2 = 1.e-10
667 qo(i,k) = max(qo(i,k), val2 )
668 heo(i,k) = .5 * g_ * (zl(i,k) + zl(i,k+1)) + &
669 cp_ * to(i,k) + hvap_ * qo(i,k)
670 heso(i,k) = .5 * g_ * (zl(i,k) + zl(i,k+1)) + &
671 cp_ * to(i,k) + hvap_ * qeso(i,k)
672 uo(i,k) = .5 * (uo(i,k) + uo(i,k+1))
673 vo(i,k) = .5 * (vo(i,k) + vo(i,k+1))
674 endif
675 enddo
676 enddo
677 !
678 ! look for the level of free convection as cloud base
679 !
680 do i=its,ite
681 flg(i) = cnvflg(i)
682 if(flg(i)) kbcon(i) = kmax(i)
683 enddo
684 !
685 do k = kts+1, km1
686 do i=its,ite
687 if (flg(i).and.k.lt.kbm(i)) then
688 if(k.gt.kb(i).and.heo(i,kb(i)).gt.heso(i,k)) then
689 kbcon(i) = k
690 flg(i) = .false.
691 endif
692 endif
693 enddo
694 enddo
695 !
696 do i=its,ite
697 if(cnvflg(i)) then
698 if(kbcon(i).eq.kmax(i)) cnvflg(i) = .false.
699 endif
700 enddo
701 !
702 totflg = .true.
703 do i=its,ite
704 totflg = totflg .and. (.not. cnvflg(i))
705 enddo
706 if(totflg) return
707 !
708 ! determine critical convective inhibition
709 ! as a function of vertical velocity at cloud base.
710 !
711 do i=its,ite
712 if(cnvflg(i)) then
713 pdot(i) = 10.* dot(i,kbcon(i))
714 endif
715 enddo
716 !
717 do i=its,ite
718 if(cnvflg(i)) then
719 if(slimsk(i).eq.1.) then
720 w1 = w1l
721 w2 = w2l
722 w3 = w3l
723 w4 = w4l
724 else
725 w1 = w1s
726 w2 = w2s
727 w3 = w3s
728 w4 = w4s
729 endif
730 if(pdot(i).le.w4) then
731 ptem = (pdot(i) - w4) / (w3 - w4)
732 elseif(pdot(i).ge.-w4) then
733 ptem = - (pdot(i) + w4) / (w4 - w3)
734 else
735 ptem = 0.
736 endif
737 val1 = -1.
738 ptem = max(ptem,val1)
739 val2 = 1.
740 ptem = min(ptem,val2)
741 ptem = 1. - ptem
742 ptem1= .5*(cincrmax-cincrmin)
743 cincr = cincrmax - ptem * ptem1
744 tem1 = p(i,kb(i)) - p(i,kbcon(i))
745 if(tem1.gt.cincr) then
746 cnvflg(i) = .false.
747 endif
748 endif
749 enddo
750 !
751 totflg = .true.
752 do i=its,ite
753 totflg = totflg .and. (.not. cnvflg(i))
754 enddo
755 if(totflg) return
756 !
757 ! assume the detrainment rate for the updrafts to be same as
758 ! the entrainment rate at cloud base
759 !
760 do i = its,ite
761 if(cnvflg(i)) then
762 xlamud(i) = xlamue(i,kbcon(i))
763 endif
764 enddo
765 !
766 ! determine updraft mass flux for the subcloud layers
767 !
768 do k = km1, kts, -1
769 do i = its,ite
770 if (cnvflg(i)) then
771 if(k.lt.kbcon(i).and.k.ge.kb(i)) then
772 dz = zi(i,k+2) - zi(i,k+1)
773 ptem = 0.5*(xlamue(i,k)+xlamue(i,k+1))-xlamud(i)
774 eta(i,k) = eta(i,k+1) / (1. + ptem * dz)
775 endif
776 endif
777 enddo
778 enddo
779 !
780 ! compute mass flux above cloud base
781 !
782 do k = kts+1, km1
783 do i = its,ite
784 if(cnvflg(i))then
785 if(k.gt.kbcon(i).and.k.lt.kmax(i)) then
786 dz = zi(i,k+1) - zi(i,k)
787 ptem = 0.5*(xlamue(i,k)+xlamue(i,k-1))-xlamud(i)
788 eta(i,k) = eta(i,k-1) * (1 + ptem * dz)
789 endif
790 endif
791 enddo
792 enddo
793 !
794 ! compute updraft cloud property
795 !
796 do i = its,ite
797 if(cnvflg(i)) then
798 indx = kb(i)
799 hcko(i,indx) = heo(i,indx)
800 ucko(i,indx) = uo(i,indx)
801 vcko(i,indx) = vo(i,indx)
802 endif
803 enddo
804 !
805 do k = kts+1, km1
806 do i = its,ite
807 if (cnvflg(i)) then
808 if(k.gt.kb(i).and.k.lt.kmax(i)) then
809 dz = zi(i,k+1) - zi(i,k)
810 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
811 tem1 = 0.5 * xlamud(i) * dz
812 factor = 1. + tem - tem1
813 ptem = 0.5 * tem + pgcon
814 ptem1= 0.5 * tem - pgcon
815 hcko(i,k) = ((1.-tem1)*hcko(i,k-1)+tem*0.5* &
816 (heo(i,k)+heo(i,k-1)))/factor
817 ucko(i,k) = ((1.-tem1)*ucko(i,k-1)+ptem*uo(i,k) &
818 +ptem1*uo(i,k-1))/factor
819 vcko(i,k) = ((1.-tem1)*vcko(i,k-1)+ptem*vo(i,k) &
820 +ptem1*vo(i,k-1))/factor
821 dbyo(i,k) = hcko(i,k) - heso(i,k)
822 endif
823 endif
824 enddo
825 enddo
826 !
827 ! taking account into convection inhibition due to existence of
828 ! dry layers below cloud base
829 !
830 do i=its,ite
831 flg(i) = cnvflg(i)
832 kbcon1(i) = kmax(i)
833 enddo
834 !
835 do k = kts+1, km1
836 do i=its,ite
837 if (flg(i).and.k.lt.kbm(i)) then
838 if(k.ge.kbcon(i).and.dbyo(i,k).gt.0.) then
839 kbcon1(i) = k
840 flg(i) = .false.
841 endif
842 endif
843 enddo
844 enddo
845 !
846 do i=its,ite
847 if(cnvflg(i)) then
848 if(kbcon1(i).eq.kmax(i)) cnvflg(i) = .false.
849 endif
850 enddo
851 !
852 do i=its,ite
853 if(cnvflg(i)) then
854 tem = p(i,kbcon(i)) - p(i,kbcon1(i))
855 if(tem.gt.dthk) then
856 cnvflg(i) = .false.
857 endif
858 endif
859 enddo
860 !
861 totflg = .true.
862 do i = its,ite
863 totflg = totflg .and. (.not. cnvflg(i))
864 enddo
865 if(totflg) return
866 !
867 ! determine first guess cloud top as the level of zero buoyancy
868 ! limited to the level of sigma=0.7
869 !
870 do i = its,ite
871 flg(i) = cnvflg(i)
872 if(flg(i)) ktcon(i) = kbm(i)
873 enddo
874 !
875 do k = kts+1, km1
876 do i=its,ite
877 if (flg(i).and.k .lt. kbm(i)) then
878 if(k.gt.kbcon1(i).and.dbyo(i,k).lt.0.) then
879 ktcon(i) = k
880 flg(i) = .false.
881 endif
882 endif
883 enddo
884 enddo
885 !
886 ! specify upper limit of mass flux at cloud base
887 !
888 do i = its,ite
889 if(cnvflg(i)) then
890 k = kbcon(i)
891 dp = 1000. * del(i,k)
892 xmbmax(i) = dp / (g_ * dt2)
893 endif
894 enddo
895 !
896 ! compute cloud moisture property and precipitation
897 !
898 do i = its,ite
899 if (cnvflg(i)) then
900 aa1(i) = 0.
901 qcko(i,kb(i)) = qo(i,kb(i))
902 endif
903 enddo
904 !
905 do k = kts+1, km1
906 do i = its,ite
907 if (cnvflg(i)) then
908 if(k.gt.kb(i).and.k.lt.ktcon(i)) then
909 dz = zi(i,k+1) - zi(i,k)
910 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
911 qrch = qeso(i,k) + gamma * dbyo(i,k) / (hvap_ * (1. + gamma))
912 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
913 tem1 = 0.5 * xlamud(i) * dz
914 factor = 1. + tem - tem1
915 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5* &
916 (qo(i,k)+qo(i,k-1)))/factor
917 dq = eta(i,k) * (qcko(i,k) - qrch)
918 !
919 ! rhbar(i) = rhbar(i) + qo(i,k) / qeso(i,k)
920 !
921 ! below lfc check if there is excess moisture to release latent heat
922 !
923 if(k.ge.kbcon(i).and.dq.gt.0.) then
924 etah = .5 * (eta(i,k) + eta(i,k-1))
925 if(ncloud.gt.0) then
926 dp = 1000. * del(i,k)
927 qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
928 dellal(i,k) = etah * c1 * dz * qlk * g_ / dp
929 else
930 qlk = dq / (eta(i,k) + etah * c0 * dz)
931 endif
932 aa1(i) = aa1(i) - dz * g_ * qlk
933 qcko(i,k)= qlk + qrch
934 pwo(i,k) = etah * c0 * dz * qlk
935 endif
936 endif
937 endif
938 enddo
939 enddo
940 !
941 ! calculate cloud work function
942 !
943 do k = kts+1, km1
944 do i = its,ite
945 if (cnvflg(i)) then
946 if(k.ge.kbcon(i).and.k.lt.ktcon(i)) then
947 dz1 = zl(i,k+1) - zl(i,k)
948 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
949 rfact = 1. + fv_ * cp_ * gamma * to(i,k) / hvap_
950 aa1(i) = aa1(i) + dz1 * (g_ / (cp_ * to(i,k))) &
951 * dbyo(i,k) / (1. + gamma) * rfact
952 val = 0.
953 aa1(i)=aa1(i)+ dz1 * g_ * fv_ * max(val,(qeso(i,k) - qo(i,k)))
954 endif
955 endif
956 enddo
957 enddo
958 !
959 do i = its,ite
960 if(cnvflg(i).and.aa1(i).le.0.) cnvflg(i) = .false.
961 enddo
962 !
963 totflg = .true.
964 do i=its,ite
965 totflg = totflg .and. (.not. cnvflg(i))
966 enddo
967 if(totflg) return
968 !
969 ! estimate the convective overshooting as the level
970 ! where the [aafac * cloud work function] becomes zero,
971 ! which is the final cloud top limited to the level of sigma=0.7
972 !
973 do i = its,ite
974 if (cnvflg(i)) then
975 aa1(i) = aafac * aa1(i)
976 endif
977 enddo
978 !
979 do i = its,ite
980 flg(i) = cnvflg(i)
981 ktcon1(i) = kbm(i)
982 enddo
983 !
984 do k = kts+1,km1
985 do i = its,ite
986 if (flg(i)) then
987 if(k.ge.ktcon(i).and.k.lt.kbm(i)) then
988 dz1 = zl(i,k+1) - zl(i,k)
989 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
990 rfact = 1. + fv_ * cp_ * gamma &
991 * to(i,k) / hvap_
992 aa1(i) = aa1(i) + &
993 dz1 * (g_ / (cp_ * to(i,k))) &
994 * dbyo(i,k) / (1. + gamma) * rfact
995 if(aa1(i).lt.0.) then
996 ktcon1(i) = k
997 flg(i) = .false.
998 endif
999 endif
1000 endif
1001 enddo
1002 enddo
1003 !
1004 ! compute cloud moisture property, detraining cloud water
1005 ! and precipitation in overshooting layers
1006 !
1007 do k = kts+1,km1
1008 do i = its,ite
1009 if (cnvflg(i)) then
1010 if(k.ge.ktcon(i).and.k.lt.ktcon1(i)) then
1011 dz = zi(i,k+1) - zi(i,k)
1012 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1013 qrch = qeso(i,k) &
1014 + gamma * dbyo(i,k) / (hvap_ * (1. + gamma))
1015 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1)) * dz
1016 tem1 = 0.5 * xlamud(i) * dz
1017 factor = 1. + tem - tem1
1018 qcko(i,k) = ((1.-tem1)*qcko(i,k-1)+tem*0.5* &
1019 (qo(i,k)+qo(i,k-1)))/factor
1020 dq = eta(i,k) * (qcko(i,k) - qrch)
1021 !
1022 ! check if there is excess moisture to release latent heat
1023 !
1024 if(dq.gt.0.) then
1025 etah = .5 * (eta(i,k) + eta(i,k-1))
1026 if(ncloud.gt.0) then
1027 dp = 1000. * del(i,k)
1028 qlk = dq / (eta(i,k) + etah * (c0 + c1) * dz)
1029 dellal(i,k) = etah * c1 * dz * qlk * g_ / dp
1030 else
1031 qlk = dq / (eta(i,k) + etah * c0 * dz)
1032 endif
1033 qcko(i,k) = qlk + qrch
1034 pwo(i,k) = etah * c0 * dz * qlk
1035 endif
1036 endif
1037 endif
1038 enddo
1039 enddo
1040 !
1041 ! exchange ktcon with ktcon1
1042 !
1043 do i = its,ite
1044 if(cnvflg(i)) then
1045 kk = ktcon(i)
1046 ktcon(i) = ktcon1(i)
1047 ktcon1(i) = kk
1048 endif
1049 enddo
1050 !
1051 ! this section is ready for cloud water
1052 !
1053 if(ncloud.gt.0) then
1054 !
1055 ! compute liquid and vapor separation at cloud top
1056 !
1057 do i = its,ite
1058 if(cnvflg(i)) then
1059 k = ktcon(i) - 1
1060 gamma = el2orc * qeso(i,k) / (to(i,k)**2)
1061 qrch = qeso(i,k) &
1062 + gamma * dbyo(i,k) / (hvap_ * (1. + gamma))
1063 dq = qcko(i,k) - qrch
1064 !
1065 ! check if there is excess moisture to release latent heat
1066 !
1067 if(dq.gt.0.) then
1068 qlko_ktcon(i) = dq
1069 qcko(i,k) = qrch
1070 endif
1071 endif
1072 enddo
1073 !
1074 endif
1075 !
1076 !--- compute precipitation efficiency in terms of windshear
1077 !
1078 do i = its,ite
1079 if(cnvflg(i)) then
1080 vshear(i) = 0.
1081 endif
1082 enddo
1083 !
1084 do k = kts+1,kte
1085 do i = its,ite
1086 if (cnvflg(i)) then
1087 if(k.gt.kb(i).and.k.le.ktcon(i)) then
1088 shear= sqrt((uo(i,k)-uo(i,k-1)) ** 2 + (vo(i,k)-vo(i,k-1)) ** 2)
1089 vshear(i) = vshear(i) + shear
1090 endif
1091 endif
1092 enddo
1093 enddo
1094 !
1095 do i = its,ite
1096 if(cnvflg(i)) then
1097 vshear(i) = 1.e3 * vshear(i) / (zi(i,ktcon(i)+1)-zi(i,kb(i)+1))
1098 e1=1.591-.639*vshear(i) &
1099 +.0953*(vshear(i)**2)-.00496*(vshear(i)**3)
1100 edt(i)=1.-e1
1101 val = .9
1102 edt(i) = min(edt(i),val)
1103 val = .0
1104 edt(i) = max(edt(i),val)
1105 endif
1106 enddo
1107 !
1108 !--- what would the change be, that a cloud with unit mass
1109 !--- will do to the environment?
1110 !
1111 do k = kts,kte
1112 do i = its,ite
1113 if(cnvflg(i) .and. k .le. kmax(i)) then
1114 dellah(i,k) = 0.
1115 dellaq(i,k) = 0.
1116 dellau(i,k) = 0.
1117 dellav(i,k) = 0.
1118 endif
1119 enddo
1120 enddo
1121 !
1122 !--- changed due to subsidence and entrainment
1123 !
1124 do k = kts+1,km1
1125 do i = its,ite
1126 if (cnvflg(i)) then
1127 if(k.gt.kb(i).and.k.lt.ktcon(i)) then
1128 dp = 1000. * del(i,k)
1129 dz = zi(i,k+1) - zi(i,k)
1130 !
1131 dv1h = heo(i,k)
1132 dv2h = .5 * (heo(i,k) + heo(i,k-1))
1133 dv3h = heo(i,k-1)
1134 dv1q = qo(i,k)
1135 dv2q = .5 * (qo(i,k) + qo(i,k-1))
1136 dv3q = qo(i,k-1)
1137 dv1u = uo(i,k)
1138 dv2u = .5 * (uo(i,k) + uo(i,k-1))
1139 dv3u = uo(i,k-1)
1140 dv1v = vo(i,k)
1141 dv2v = .5 * (vo(i,k) + vo(i,k-1))
1142 dv3v = vo(i,k-1)
1143 !
1144 tem = 0.5 * (xlamue(i,k)+xlamue(i,k-1))
1145 tem1 = xlamud(i)
1146 !
1147 dellah(i,k) = dellah(i,k) + &
1148 ( eta(i,k)*dv1h - eta(i,k-1)*dv3h &
1149 - tem*eta(i,k-1)*dv2h*dz &
1150 + tem1*eta(i,k-1)*.5*(hcko(i,k)+hcko(i,k-1))*dz ) *g_/dp
1151 !
1152 dellaq(i,k) = dellaq(i,k) + &
1153 ( eta(i,k)*dv1q - eta(i,k-1)*dv3q &
1154 - tem*eta(i,k-1)*dv2q*dz &
1155 + tem1*eta(i,k-1)*.5*(qcko(i,k)+qcko(i,k-1))*dz ) *g_/dp
1156 !
1157 dellau(i,k) = dellau(i,k) + &
1158 ( eta(i,k)*dv1u - eta(i,k-1)*dv3u &
1159 - tem*eta(i,k-1)*dv2u*dz &
1160 + tem1*eta(i,k-1)*.5*(ucko(i,k)+ucko(i,k-1))*dz &
1161 - pgcon*eta(i,k-1)*(dv1u-dv3u) ) *g_/dp
1162 !
1163 dellav(i,k) = dellav(i,k) + &
1164 ( eta(i,k)*dv1v - eta(i,k-1)*dv3v &
1165 - tem*eta(i,k-1)*dv2v*dz &
1166 + tem1*eta(i,k-1)*.5*(vcko(i,k)+vcko(i,k-1))*dz &
1167 - pgcon*eta(i,k-1)*(dv1v-dv3v) ) *g_/dp
1168 !
1169 endif
1170 endif
1171 enddo
1172 enddo
1173 !
1174 !------- cloud top
1175 !
1176 do i = its,ite
1177 if(cnvflg(i)) then
1178 indx = ktcon(i)
1179 dp = 1000. * del(i,indx)
1180 dv1h = heo(i,indx-1)
1181 dellah(i,indx) = eta(i,indx-1) * &
1182 (hcko(i,indx-1) - dv1h) * g_ / dp
1183 dv1q = qo(i,indx-1)
1184 dellaq(i,indx) = eta(i,indx-1) * &
1185 (qcko(i,indx-1) - dv1q) * g_ / dp
1186 dv1u = uo(i,indx-1)
1187 dellau(i,indx) = eta(i,indx-1) * &
1188 (ucko(i,indx-1) - dv1u) * g_ / dp
1189 dv1v = vo(i,indx-1)
1190 dellav(i,indx) = eta(i,indx-1) * &
1191 (vcko(i,indx-1) - dv1v) * g_ / dp
1192 !
1193 ! cloud water
1194 !
1195 dellal(i,indx) = eta(i,indx-1) * &
1196 qlko_ktcon(i) * g_ / dp
1197 endif
1198 enddo
1199 !
1200 ! mass flux at cloud base for shallow convection
1201 ! (Grant, 2001)
1202 !
1203 do i= its,ite
1204 if(cnvflg(i)) then
1205 k = kbcon(i)
1206 ptem = g_*sflx(i)*hpbl(i)/t1(i,1)
1207 wstar(i) = ptem**h1
1208 tem = po(i,k)*100. / (rd_*t1(i,k))
1209 xmb(i) = betaw*tem*wstar(i)
1210 xmb(i) = min(xmb(i),xmbmax(i))
1211 endif
1212 enddo
1213 !
1214 do k = kts,kte
1215 do i = its,ite
1216 if (cnvflg(i) .and. k .le. kmax(i)) then
1217 qeso(i,k)=0.01* fpvs(t1(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls,psat,t0c_)
1218 qeso(i,k) = eps * qeso(i,k) / (p(i,k) + (eps-1.)*qeso(i,k))
1219 val = 1.e-8
1220 qeso(i,k) = max(qeso(i,k), val )
1221 endif
1222 enddo
1223 enddo
1224 !
1225 do i = its,ite
1226 delhbar(i) = 0.
1227 delqbar(i) = 0.
1228 deltbar(i) = 0.
1229 delubar(i) = 0.
1230 delvbar(i) = 0.
1231 qcond(i) = 0.
1232 enddo
1233 !
1234 do k = kts,kte
1235 do i = its,ite
1236 if (cnvflg(i)) then
1237 if(k.gt.kb(i).and.k.le.ktcon(i)) then
1238 dellat = (dellah(i,k) - hvap_ * dellaq(i,k)) / cp_
1239 t1(i,k) = t1(i,k) + dellat * xmb(i) * dt2
1240 q1(i,k) = q1(i,k) + dellaq(i,k) * xmb(i) * dt2
1241 tem = 1./rcs
1242 u1(i,k) = u1(i,k) + dellau(i,k) * xmb(i) * dt2 * tem
1243 v1(i,k) = v1(i,k) + dellav(i,k) * xmb(i) * dt2 * tem
1244 dp = 1000. * del(i,k)
1245 delhbar(i) = delhbar(i) + dellah(i,k)*xmb(i)*dp/g_
1246 delqbar(i) = delqbar(i) + dellaq(i,k)*xmb(i)*dp/g_
1247 deltbar(i) = deltbar(i) + dellat*xmb(i)*dp/g_
1248 delubar(i) = delubar(i) + dellau(i,k)*xmb(i)*dp/g_
1249 delvbar(i) = delvbar(i) + dellav(i,k)*xmb(i)*dp/g_
1250 endif
1251 endif
1252 enddo
1253 enddo
1254 !
1255 do k = kts,kte
1256 do i = its,ite
1257 if (cnvflg(i)) then
1258 if(k.gt.kb(i).and.k.le.ktcon(i)) then
1259 qeso(i,k)=0.01* fpvs(t1(i,k),1,rd_,rv_,cvap_,cliq_,cice,xlv0,xls &
1260 ,psat,t0c_)
1261 qeso(i,k) = eps * qeso(i,k)/(p(i,k) + (eps-1.)*qeso(i,k))
1262 val = 1.e-8
1263 qeso(i,k) = max(qeso(i,k), val )
1264 endif
1265 endif
1266 enddo
1267 enddo
1268 !
1269 do i = its,ite
1270 rntot(i) = 0.
1271 delqev(i) = 0.
1272 delq2(i) = 0.
1273 flg(i) = cnvflg(i)
1274 enddo
1275 !
1276 do k = kte,kts,-1
1277 do i = its,ite
1278 if (cnvflg(i)) then
1279 if(k.lt.ktcon(i).and.k.gt.kb(i)) then
1280 rntot(i) = rntot(i) + pwo(i,k) * xmb(i) * .001 * dt2
1281 endif
1282 endif
1283 enddo
1284 enddo
1285 !
1286 ! evaporating rain
1287 !
1288 do k = kte,kts,-1
1289 do i = its,ite
1290 if (k .le. kmax(i)) then
1291 deltv(i) = 0.
1292 delq(i) = 0.
1293 qevap(i) = 0.
1294 if(cnvflg(i)) then
1295 if(k.lt.ktcon(i).and.k.gt.kb(i)) then
1296 rain(i) = rain(i) + pwo(i,k) * xmb(i) * .001 * dt2
1297 endif
1298 endif
1299 if(flg(i).and.k.lt.ktcon(i)) then
1300 evef = edt(i) * evfact
1301 if(slimsk(i).eq.1.) evef=edt(i) * evfactl
1302 qcond(i) = evef * (q1(i,k) - qeso(i,k)) &
1303 / (1. + el2orc * qeso(i,k) / t1(i,k)**2)
1304 dp = 1000. * del(i,k)
1305 if(rain(i).gt.0..and.qcond(i).lt.0.) then
1306 qevap(i) = -qcond(i) * (1.-exp(-.32*sqrt(dt2*rain(i))))
1307 qevap(i) = min(qevap(i), rain(i)*1000.*g_/dp)
1308 delq2(i) = delqev(i) + .001 * qevap(i) * dp / g_
1309 endif
1310 if(rain(i).gt.0..and.qcond(i).lt.0..and.delq2(i).gt.rntot(i)) then
1311 qevap(i) = 1000.* g_ * (rntot(i) - delqev(i)) / dp
1312 flg(i) = .false.
1313 endif
1314 if(rain(i).gt.0..and.qevap(i).gt.0.) then
1315 tem = .001 * dp / g_
1316 tem1 = qevap(i) * tem
1317 if(tem1.gt.rain(i)) then
1318 qevap(i) = rain(i) / tem
1319 rain(i) = 0.
1320 else
1321 rain(i) = rain(i) - tem1
1322 endif
1323 q1(i,k) = q1(i,k) + qevap(i)
1324 t1(i,k) = t1(i,k) - (hvap_/cp_) * qevap(i)
1325 deltv(i) = - (hvap_/cp_)*qevap(i)/dt2
1326 delq(i) = + qevap(i)/dt2
1327 delqev(i) = delqev(i) + .001*dp*qevap(i)/g_
1328 endif
1329 dellaq(i,k) = dellaq(i,k) + delq(i) / xmb(i)
1330 delqbar(i) = delqbar(i) + delq(i)*dp/g_
1331 deltbar(i) = deltbar(i) + deltv(i)*dp/g_
1332 endif
1333 endif
1334 enddo
1335 enddo
1336 !
1337 do i = its,ite
1338 if(cnvflg(i)) then
1339 if(rain(i).lt.0..or..not.flg(i)) rain(i) = 0.
1340 ktop(i) = ktcon(i)
1341 kbot(i) = kbcon(i)
1342 icps(i) = 0
1343 endif
1344 enddo
1345 !
1346 ! cloud water
1347 !
1348 if (ncloud.gt.0) then
1349 !
1350 do k = kts,km1
1351 do i = its,ite
1352 if (cnvflg(i)) then
1353 if (k.ge.kbcon(i).and.k.le.ktcon(i)) then
1354 tem = dellal(i,k) * xmb(i) * dt2
1355 tem1 = max(0.0, min(1.0, (tcr-t1(i,k))*tcrf))
1356 if (ncloud.ge.2) then
1357 qi2(i,k) = qi2(i,k) + tem * tem1 ! ice
1358 qc2(i,k) = qc2(i,k) + tem *(1.0-tem1) ! water
1359 else
1360 qc2(i,k) = qc2(i,k) + tem
1361 endif
1362 endif
1363 endif
1364 enddo
1365 enddo
1366 !
1367 endif
1368 !
1369 end subroutine nscv2d
1370 !-------------------------------------------------------------------------------
1371 !
1372 !-------------------------------------------------------------------------------
1373 REAL FUNCTION fpvs(t,ice,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c)
1374 !-------------------------------------------------------------------------------
1375 IMPLICIT NONE
1376 !-------------------------------------------------------------------------------
1377 REAL :: t,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c,dldt,xa,xb,dldti, &
1378 xai,xbi,ttp,tr
1379 INTEGER :: ice
1380 !
1381 ttp=t0c+0.01
1382 dldt=cvap-cliq
1383 xa=-dldt/rv
1384 xb=xa+hvap/(rv*ttp)
1385 dldti=cvap-cice
1386 xai=-dldti/rv
1387 xbi=xai+hsub/(rv*ttp)
1388 tr=ttp/t
1389 if(t.lt.ttp.and.ice.eq.1) then
1390 fpvs=psat*(tr**xai)*exp(xbi*(1.-tr))
1391 else
1392 fpvs=psat*(tr**xa)*exp(xb*(1.-tr))
1393 endif
1394 !
1395 if (t.lt.180.) then
1396 tr=ttp/180.
1397 if(t.lt.ttp.and.ice.eq.1) then
1398 fpvs=psat*(tr**xai)*exp(xbi*(1.-tr))
1399 else
1400 fpvs=psat*(tr**xa)*exp(xb*(1.-tr))
1401 endif
1402 endif
1403 !
1404 if (t.ge.330.) then
1405 tr=ttp/330
1406 if(t.lt.ttp.and.ice.eq.1) then
1407 fpvs=psat*(tr**xai)*exp(xbi*(1.-tr))
1408 else
1409 fpvs=psat*(tr**xa)*exp(xb*(1.-tr))
1410 endif
1411 endif
1412 !
1413 END FUNCTION fpvs
1414 !-------------------------------------------------------------------------------
1415 !
1416 !-------------------------------------------------------------------------------
1417 subroutine nscvinit(rthshten,rqvshten,rqcshten,rqishten, &
1418 rushten,rvshten, &
1419 restart,p_qc,p_qi,p_first_scalar, &
1420 allowed_to_read, &
1421 ids, ide, jds, jde, kds, kde, &
1422 ims, ime, jms, jme, kms, kme, &
1423 its, ite, jts, jte, kts, kte )
1424 !-------------------------------------------------------------------------------
1425 implicit none
1426 !-------------------------------------------------------------------------------
1427 logical , intent(in) :: allowed_to_read,restart
1428 integer , intent(in) :: ids, ide, jds, jde, kds, kde, &
1429 ims, ime, jms, jme, kms, kme, &
1430 its, ite, jts, jte, kts, kte
1431 integer , intent(in) :: p_first_scalar, p_qi, p_qc
1432 real, dimension( ims:ime , kms:kme , jms:jme ) , intent(out) :: &
1433 rthshten, &
1434 rqvshten, &
1435 rushten, &
1436 rvshten, &
1437 rqcshten, &
1438 rqishten
1439 integer :: i, j, k, itf, jtf, ktf
1440 !
1441 jtf=min0(jte,jde-1)
1442 ktf=min0(kte,kde-1)
1443 itf=min0(ite,ide-1)
1444 !
1445 if(.not.restart)then
1446 do j = jts,jtf
1447 do k = kts,ktf
1448 do i = its,itf
1449 rthshten(i,k,j)=0.
1450 rqvshten(i,k,j)=0.
1451 rushten(i,k,j)=0.
1452 rvshten(i,k,j)=0.
1453 enddo
1454 enddo
1455 enddo
1456 !
1457 if (p_qc .ge. p_first_scalar) then
1458 do j = jts,jtf
1459 do k = kts,ktf
1460 do i = its,itf
1461 rqcshten(i,k,j)=0.
1462 enddo
1463 enddo
1464 enddo
1465 endif
1466 !
1467 if (p_qi .ge. p_first_scalar) then
1468 do j = jts,jtf
1469 do k = kts,ktf
1470 do i = its,itf
1471 rqishten(i,k,j)=0.
1472 enddo
1473 enddo
1474 enddo
1475 endif
1476 endif
1477 !
1478 end subroutine nscvinit
1479 !-------------------------------------------------------------------------------
1480 !
1481 END MODULE module_shcu_nscv
1482 !