module_mp_ncloud5.F
References to this file elsewhere.
1 MODULE module_mp_ncloud5
2
3 REAL, PARAMETER, PRIVATE :: dtcldcr = 240.
4 INTEGER, PARAMETER, PRIVATE :: mstepmax = 100
5
6 REAL, PARAMETER, PRIVATE :: r0 = .8e-5 ! 8 microm in contrast to 10 micro m
7 REAL, PARAMETER, PRIVATE :: n0r = 8.e6
8 REAL, PARAMETER, PRIVATE :: avtr = 841.9
9 REAL, PARAMETER, PRIVATE :: bvtr = 0.8
10 REAL, PARAMETER, PRIVATE :: peaut = .55 ! collection efficiency
11 REAL, PARAMETER, PRIVATE :: xncr = 3.e8 ! maritime cloud in contrast to 3.e8 in tc80
12 REAL, PARAMETER, PRIVATE :: xmyu = 1.718e-5 ! the dynamic viscosity kgm-1s-1
13
14 REAL, PARAMETER, PRIVATE :: avts = 16.2
15 REAL, PARAMETER, PRIVATE :: bvts = .527
16 REAL, PARAMETER, PRIVATE :: xncmax = 1.e8
17 REAL, PARAMETER, PRIVATE :: n0smax = 1.e9
18 REAL, PARAMETER, PRIVATE :: betai = .6
19 REAL, PARAMETER, PRIVATE :: xn0 = 1.e-2
20 REAL, PARAMETER, PRIVATE :: dicon = 16.3
21 REAL, PARAMETER, PRIVATE :: di0 = 12.9e-6*.8
22 REAL, PARAMETER, PRIVATE :: dimax = 400.e-6
23 REAL, PARAMETER, PRIVATE :: n0s = 2.e6 ! temperature dependent n0s
24 REAL, PARAMETER, PRIVATE :: alpha = 1./8.18 ! .122 exponen factor for n0s
25 REAL, PARAMETER, PRIVATE :: pfrz1 = 100.
26 REAL, PARAMETER, PRIVATE :: pfrz2 = 0.66
27 ! REAL, PARAMETER, PRIVATE :: lamdarmax = 1.e15
28 REAL, PARAMETER, PRIVATE :: lamdarmax = 1.e5
29 REAL, PARAMETER, PRIVATE :: qcrmin = 1.e-6
30
31 REAL, PARAMETER, PRIVATE :: t40c = 233.16
32 REAL, PARAMETER, PRIVATE :: eacrc = 1.0
33
34 REAL, SAVE :: &
35 qc0, qck1,bvtr1,bvtr2,bvtr3,bvtr4,g1pbr,&
36 g3pbr,g4pbr,g5pbro2,pvtr,eacrr,pacrr, &
37 precr1,precr2,xm0,xmmax,bvts1, &
38 bvts2,bvts3,bvts4,g1pbs,g3pbs,g4pbs, &
39 g5pbso2,pvts,pacrs,precs1,precs2,pidn0r,&
40 pidn0s,xlv1,pacrc
41
42 CONTAINS
43
44 !===================================================================
45 !
46 SUBROUTINE ncloud5(th, q, qc, qr, qi, qs &
47 ,w, den, pii, p, delz &
48 ,delt,g, cpd, cpv, rd, rv, t0c &
49 ,ep1, ep2, qmin &
50 ,XLS, XLV0, XLF0, den0, denr &
51 ,cliq,cice,psat &
52 ,rain, rainncv &
53 ,ids,ide, jds,jde, kds,kde &
54 ,ims,ime, jms,jme, kms,kme &
55 ,its,ite, jts,jte, kts,kte &
56 )
57 !-------------------------------------------------------------------
58 IMPLICIT NONE
59 !-------------------------------------------------------------------
60 !
61 !Coded by Song-You Hong (NCEP) and implemented by Shuhua Chen (NCAR)
62 !
63 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
64 ims,ime, jms,jme, kms,kme , &
65 its,ite, jts,jte, kts,kte
66
67 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
68 INTENT(INOUT) :: &
69 th, &
70 q, &
71 qc, &
72 qi, &
73 qr, &
74 qs
75
76 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
77 INTENT(IN ) :: w, &
78 den, &
79 pii, &
80 p, &
81 delz
82
83
84 REAL, INTENT(IN ) :: delt, &
85 g, &
86 rd, &
87 rv, &
88 t0c, &
89 den0, &
90 cpd, &
91 cpv, &
92 ep1, &
93 ep2, &
94 qmin, &
95 xls, &
96 xlv0, &
97 xlf0, &
98 cliq, &
99 cice, &
100 psat, &
101 denr
102 REAL, DIMENSION( ims:ime , jms:jme ), &
103 INTENT(INOUT) :: rain, &
104 rainncv
105
106 ! LOCAL VAR
107
108 REAL, DIMENSION( its:ite , kts:kte ) :: t
109 REAL, DIMENSION( its:ite , kts:kte, 2 ) :: qci, qrs
110 INTEGER :: i,j,k
111
112 !-------------------------------------------------------------------
113 DO J=jts,jte
114
115 DO k=kts,kte
116 DO i=its,ite
117 t(i,k)=th(i,k,j)*pii(i,k,j)
118 qci(i,k,1) = qc(i,k,j)
119 qci(i,k,2) = qi(i,k,j)
120 qrs(i,k,1) = qr(i,k,j)
121 qrs(i,k,2) = qs(i,k,j)
122 ENDDO
123 ENDDO
124
125 CALL ncloud52D(t, q(ims,kms,j), qci, qrs &
126 ,w(ims,kms,j), den(ims,kms,j) &
127 ,p(ims,kms,j), delz(ims,kms,j) &
128 ,delt,g, cpd, cpv, rd, rv, t0c &
129 ,ep1, ep2, qmin &
130 ,XLS, XLV0, XLF0, den0, denr &
131 ,cliq,cice,psat &
132 ,j &
133 ,rain(ims,j),rainncv(ims,j) &
134 ,ids,ide, jds,jde, kds,kde &
135 ,ims,ime, jms,jme, kms,kme &
136 ,its,ite, jts,jte, kts,kte &
137 )
138
139 DO K=kts,kte
140 DO I=its,ite
141 th(i,k,j)=t(i,k)/pii(i,k,j)
142 qc(i,k,j) = qci(i,k,1)
143 qi(i,k,j) = qci(i,k,2)
144 qr(i,k,j) = qrs(i,k,1)
145 qs(i,k,j) = qrs(i,k,2)
146 ENDDO
147 ENDDO
148
149 ENDDO
150
151 END SUBROUTINE ncloud5
152
153 !===================================================================
154 !
155 SUBROUTINE ncloud52D(t, q, qci, qrs,w, den, p, delz &
156 ,delt,g, cpd, cpv, rd, rv, t0c &
157 ,ep1, ep2, qmin &
158 ,XLS, XLV0, XLF0, den0, denr &
159 ,cliq,cice,psat &
160 ,lat &
161 ,rain,rainncv &
162 ,ids,ide, jds,jde, kds,kde &
163 ,ims,ime, jms,jme, kms,kme &
164 ,its,ite, jts,jte, kts,kte &
165 )
166 !-------------------------------------------------------------------
167 IMPLICIT NONE
168 !-------------------------------------------------------------------
169 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
170 ims,ime, jms,jme, kms,kme , &
171 its,ite, jts,jte, kts,kte, &
172 lat
173
174 REAL, DIMENSION( its:ite , kts:kte ), &
175 INTENT(INOUT) :: &
176 t
177 REAL, DIMENSION( its:ite , kts:kte, 2 ), &
178 INTENT(INOUT) :: &
179 qci, qrs
180
181 REAL, DIMENSION( ims:ime , kms:kme ), &
182 INTENT(INOUT) :: &
183 q
184
185 REAL, DIMENSION( ims:ime , kms:kme ), &
186 INTENT(IN ) :: p, &
187 w, &
188 den, &
189 delz
190
191 REAL, INTENT(IN ) :: delt, &
192 g, &
193 cpd, &
194 cpv, &
195 t0c, &
196 den0, &
197 rd, &
198 rv, &
199 ep1, &
200 ep2, &
201 qmin, &
202 xls, &
203 xlv0, &
204 xlf0, &
205 cliq, &
206 cice, &
207 psat, &
208 denr
209 REAL, DIMENSION( ims:ime ), &
210 INTENT(INOUT) :: rain, &
211 rainncv
212
213 ! LOCAL VAR
214
215 INTEGER, PARAMETER :: iun = 84
216
217 REAL, DIMENSION( its:ite , kts:kte , 2) :: &
218 rh, qs, slope, slope2, slopeb, &
219 paut, pres, falk, fall, work1
220 REAL, DIMENSION( its:ite , kts:kte ) :: &
221 falkc, work1c, work2c, fallc
222 REAL, DIMENSION( its:ite , kts:kte, 3 ) :: &
223 pacr
224 REAL, DIMENSION( its:ite , kts:kte ) :: &
225 pgen, pisd, pcon, xl, cpm, work2, q1, t1, denfac, &
226 pgens, pauts, pacrss, pisds, press, pcons, psml, psev
227
228 INTEGER, DIMENSION( its:ite ) :: mstep
229 LOGICAL, DIMENSION( its:ite ) :: flgcld
230
231 REAL :: n0sfac, pi, &
232 cpmcal, xlcal, lamdar, lamdas, diffus, &
233 viscos, xka, venfac, conden, diffac, &
234 x, y, z, a, b, c, d, e, &
235 qdt, holdrr, holdrs, supcol, &
236 coeres, supsat, dtcld, xmi, eacrs, satdt, xnc, &
237 fallsum, xlwork2, factor, source, value, &
238 xlf,pfrzdtc,pfrzdtr
239
240 INTEGER :: i,j,k, &
241 iprt, latd, lond, loop, loops, ifsat, n, numdt
242
243 !
244 !=================================================================
245 ! compute internal functions
246 !
247 cpmcal(x) = cpd*(1.-max(x,qmin))+max(x,qmin)*cpv
248 xlcal(x) = xlv0-xlv1*(x-t0c)
249 ! tvcal(x,y) = x+x*ep1*max(y,qmin)
250 !----------------------------------------------------------------
251 ! size distributions: (x=mixing ratio, y=air density):
252 ! valid for mixing ratio > 1.e-30 kg/kg.
253 ! otherwise use uniform distribution value (1.e15)
254 !
255 lamdar(x,y)=(pidn0r/(x*y))**.25
256 lamdas(x,y,z)=(pidn0s*z/(x*y))**.25
257 !
258 !----------------------------------------------------------------
259 ! diffus: diffusion coefficient of the water vapor
260 ! viscos: kinematic viscosity(m2s-1)
261 !
262 diffus(x,y) = 8.794e-5*x**1.81/y
263 viscos(x,y) = 1.496e-6*x**1.5/(x+120.)/y
264 xka(x,y) = 1.414e3*viscos(x,y)*y
265 diffac(a,b,c,d,e) = d*a*a/(xka(c,d)*rv*c*c)+1./(e*diffus(c,b))
266 venfac(a,b,c) = (viscos(b,c)/diffus(b,a))**(.3333333) &
267 /viscos(b,c)**(.5)*(den0/c)**0.25
268 conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+d*d/(rv*e)*c/(a*a))
269 !
270 pi = 4. * atan(1.)
271 !
272 !=================================================================
273 ! set iprt = 0 for no unit fort.84 output
274 !
275 iprt = 1
276 if(iprt.eq.1) then
277 qdt = delt * 1000.
278 latd = (jte+jts)/2 + 1
279 lond = (ite+its)/2 + 1
280 else
281 latd = jts
282 lond = its
283 endif
284 !
285 !----------------------------------------------------------------
286 ! paddint 0 for negative values generated by dynamics
287 !
288 do k = kts, kte
289 do i = its, ite
290 qci(i,k,1) = max(qci(i,k,1),0.0)
291 qrs(i,k,1) = max(qrs(i,k,1),0.0)
292 qci(i,k,2) = max(qci(i,k,2),0.0)
293 qrs(i,k,2) = max(qrs(i,k,2),0.0)
294 enddo
295 enddo
296 !
297 !----------------------------------------------------------------
298 ! latent heat for phase changes and heat capacity. neglect the
299 ! changes during microphysical process calculation
300 ! emanuel(1994)
301 !
302 do k = kts, kte
303 do i = its, ite
304 cpm(i,k) = cpmcal(q(i,k))
305 xl(i,k) = xlcal(t(i,k))
306 enddo
307 enddo
308 if(lat.eq.latd) then
309 i = lond
310 do k = kts, kte
311 press(i,k) = 0.
312 pauts(i,k) = 0.
313 pacrss(i,k)= 0.
314 pgens(i,k) = 0.
315 pisds(i,k) = 0.
316 pcons(i,k) = 0.
317 t1(i,k) = t(i,k)
318 q1(i,k) = q(i,k)
319 enddo
320 endif
321 !
322 !----------------------------------------------------------------
323 ! compute the minor time steps.
324 !
325 loops = max(nint(delt/dtcldcr),1)
326 dtcld = delt/loops
327 if(delt.le.dtcldcr) dtcld = delt
328 !
329 do loop = 1,loops
330 !
331 do i = its, ite
332 mstep(i) = 1
333 flgcld(i) = .true.
334 enddo
335 !
336 do k = kts, kte
337 do i = its, ite
338 denfac(i,k) = sqrt(den0/den(i,k))
339 enddo
340 enddo
341 !
342 do k = kts, kte
343 do i = its, ite
344 qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
345 qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
346 qs(i,k,1) = max(qs(i,k,1),qmin)
347 rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
348 qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
349 qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
350 qs(i,k,2) = max(qs(i,k,2),qmin)
351 rh(i,k,2) = max(q(i,k) / qs(i,k,2),qmin)
352 ! if(lat.eq.latd.and.i.eq.lond) write(iun,700) qci(i,k,1)*1000., &
353 ! qrs(i,k,1)*1000.,qci(i,k,2)*1000.,qrs(i,k,2)*1000.
354 !700 format(1x,'before computation', 4f10.4)
355 enddo
356 enddo
357 !
358 !
359 !----------------------------------------------------------------
360 ! initialize the variables for microphysical physics
361 !
362 !
363 do k = kts, kte
364 do i = its, ite
365 pres(i,k,1) = 0.
366 pres(i,k,2) = 0.
367 paut(i,k,1) = 0.
368 paut(i,k,2) = 0.
369 pacr(i,k,1) = 0.
370 pacr(i,k,2) = 0.
371 pacr(i,k,3) = 0.
372 pgen(i,k) = 0.
373 pisd(i,k) = 0.
374 pcon(i,k) = 0.
375 psml(i,k) = 0.
376 psev(i,k) = 0.
377 falk(i,k,1) = 0.
378 falk(i,k,2) = 0.
379 fall(i,k,1) = 0.
380 fall(i,k,2) = 0.
381 fallc(i,k) = 0.
382 falkc(i,k) = 0.
383 enddo
384 enddo
385 !
386 !----------------------------------------------------------------
387 ! sloper: the slope parameter of the rain(m-1)
388 ! xka: thermal conductivity of air(jm-1s-1k-1)
389 ! work1: the thermodynamic term in the denominator associated with
390 ! heat conduction and vapor diffusion
391 ! (ry88, y93, h85)
392 ! work2: parameter associated with the ventilation effects(y93)
393 !
394 do k = kts, kte
395 do i = its, ite
396 if(qrs(i,k,1).le.qcrmin)then
397 slope(i,k,1) = lamdarmax
398 slopeb(i,k,1) = slope(i,k,1)**bvtr
399 else
400 slope(i,k,1) = lamdar(qrs(i,k,1),den(i,k))
401 slopeb(i,k,1) = slope(i,k,1)**bvtr
402 endif
403 if(qrs(i,k,2).le.qcrmin)then
404 slope(i,k,2) = lamdarmax
405 slopeb(i,k,2) = slope(i,k,2)**bvts
406 else
407 supcol = t0c-t(i,k)
408 n0sfac = min(exp(alpha*supcol),n0smax)
409 slope(i,k,2) = lamdas(qrs(i,k,2),den(i,k),n0sfac)
410 slopeb(i,k,2) = slope(i,k,2)**bvts
411 endif
412 slope2(i,k,1) = slope(i,k,1)*slope(i,k,1)
413 slope2(i,k,2) = slope(i,k,2)*slope(i,k,2)
414 enddo
415 enddo
416 !
417 do k = kts, kte
418 do i = its, ite
419 work1(i,k,1) = diffac(xl(i,k),p(i,k),t(i,k),den(i,k),qs(i,k,1))
420 work1(i,k,2) = diffac(xls,p(i,k),t(i,k),den(i,k),qs(i,k,2))
421 work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
422 enddo
423 enddo
424 !
425 ! instantaneous melting of cloud ice
426 !
427 do k = kts, kte
428 do i = its, ite
429 supcol = t0c-t(i,k)
430 xlf = xls-xl(i,k)
431 if(supcol.lt.0..and.qci(i,k,2).gt.qcrmin) then
432 qci(i,k,1) = qci(i,k,1) + qci(i,k,2)
433 t(i,k) = t(i,k) - xlf/cpm(i,k)*qci(i,k,2)
434 qci(i,k,2) = 0.
435 ! if(lat.eq.latd.and.i.eq.lond) write(iun,607) k,-supcol, &
436 ! qci(i,k,2)*1000.,xlf/cpm(i,k)*qci(i,k,2)
437 endif
438 !607 format(1x,'k = ',i3,' t = ',f5.1,' qi melting = ',f10.6, &
439 ! ' del t = ',f10.6)
440 !
441 ! homogeneous freezing of cloud water below -40c
442 !
443 if(supcol.gt.40..and.qci(i,k,1).gt.qcrmin) then
444 qci(i,k,2) = qci(i,k,2) + qci(i,k,1)
445 t(i,k) = t(i,k) + xlf/cpm(i,k)*qci(i,k,1)
446 qci(i,k,1) = 0.
447 ! if(lat.eq.latd.and.i.eq.lond) write(iun,608) k,-supcol, &
448 ! qci(i,k,1)*1000.,xlf/cpm(i,k)*qci(i,k,1)
449 endif
450 !608 format(1x,'k = ',i3,' t = ',f5.1,' qc home freezing = ',f10.6, &
451 ! ' del t = ',f10.6)
452 !
453 ! heterogeneous freezing of cloud water
454 !
455 if(supcol.gt.0..and.qci(i,k,1).gt.qcrmin) then
456 pfrzdtc = min(pfrz1*exp(pfrz2*supcol-1.) &
457 *den(i,k)/denr/xncr*qci(i,k,1)**2*dtcld,qci(i,k,1))
458 qci(i,k,2) = qci(i,k,2) + pfrzdtc
459 t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtc
460 qci(i,k,1) = qci(i,k,1)-pfrzdtc
461 ! if(lat.eq.latd.and.i.eq.lond) write(iun,609) k,-supcol, &
462 ! pfrzdtc*1000.,xlf/cpm(i,k)*pfrzdtc
463 endif
464 !609 format(1x,'k = ',i3,' t = ',f5.1,' qc hete freezing = ',f10.6, &
465 ! ' del t = ',f10.6)
466 !
467 ! freezing of rain water
468 !
469 if(supcol.gt.0..and.qrs(i,k,1).gt.qcrmin) then
470 pfrzdtr = min(20.*pi**2*pfrz1*n0r*denr/den(i,k) &
471 *exp(pfrz2*supcol-1.)*slope(i,k,1)**(-7)*dtcld, &
472 qrs(i,k,1))
473 qrs(i,k,2) = qrs(i,k,2) + pfrzdtr
474 t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtr
475 qrs(i,k,1) = qrs(i,k,1)-pfrzdtr
476 ! if(lat.eq.latd.and.i.eq.lond) write(iun,610) k,-supcol, &
477 ! pfrzdtr*1000.,xlf/cpm(i,k)*pfrzdtr
478 !610 format(1x,'k = ',i3,' t = ',f5.1,' qr freezing = ',f10.6, &
479 ! ' del t = ',f10.6)
480 endif
481 enddo
482 enddo
483 !
484 ! warm rain process
485 ! paut: auto conversion rate from cloud to rain (kgkg-1s-1)
486 ! pacr: accretion rate of rain by cloud(lin83)
487 ! pres: evaporation/condensation rate of rain(rh83)
488 !
489 do k = kts, kte
490 do i = its, ite
491 supsat = max(q(i,k),qmin)-qs(i,k,1)
492 satdt = supsat/dtcld
493 if(qci(i,k,1).gt.qc0) then
494 paut(i,k,1) = qck1*qci(i,k,1)**(7./3.)
495 paut(i,k,1) = min(paut(i,k,1),qci(i,k,1)/dtcld)
496 endif
497 if(qrs(i,k,1).gt.qcrmin) then
498 if(qci(i,k,1).gt.qcrmin) pacr(i,k,1) &
499 = min(pacrr/slope2(i,k,1)/slope(i,k,1)/slopeb(i,k,1) &
500 *qci(i,k,1)*denfac(i,k),qci(i,k,1)/dtcld)
501 coeres = slope2(i,k,1)*sqrt(slope(i,k,1)*slopeb(i,k,1))
502 pres(i,k,1) = (rh(i,k,1)-1.)*(precr1/slope2(i,k,1) &
503 +precr2*work2(i,k)/coeres)/work1(i,k,1)
504 if(pres(i,k,1).lt.0.) then
505 pres(i,k,1) = max(pres(i,k,1),-qrs(i,k,1)/dtcld)
506 pres(i,k,1) = max(pres(i,k,1),satdt/2)
507 else
508 pres(i,k,1) = min(pres(i,k,1),satdt/2)
509 pres(i,k,1) = min(pres(i,k,1),qrs(i,k,1)/dtcld)
510 endif
511 endif
512 enddo
513 enddo
514 !
515 !----------------------------------------------------------------
516 ! cold rain process
517 ! paut: conversion(aggregation) of ice to snow(kgkg-1s-1)(rh83)
518 ! pgen: generation(nucleation) of ice from vapor(kgkg-1s-1)(rh83)
519 ! pacr: accretion rate of snow by ice(lin83)
520 ! pisd: deposition/sublimation rate of ice(rh83)
521 ! pres: deposition/sublimation rate of snow(lin83)
522 !
523 do k = kts, kte
524 do i = its, ite
525 supcol = t0c-t(i,k)
526 supsat = max(q(i,k),qmin)-qs(i,k,2)
527 satdt = supsat/dtcld
528 ifsat = 0
529 n0sfac = min(exp(alpha*supcol),n0smax)
530 xnc = min(xn0 * exp(betai*supcol)/den(i,k),xncmax)
531 !
532 if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,2).gt.qcrmin) then
533 eacrs = exp(0.025*(-supcol))
534 pacr(i,k,2) = min(pacrs*n0sfac*eacrs/slope2(i,k,2)/slope(i,k,2) &
535 /slopeb(i,k,2)*qci(i,k,2)*denfac(i,k),qci(i,k,2)/dtcld)
536 endif
537 if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,1).gt.qcrmin) then
538 pacr(i,k,3) = min(pacrc/slope2(i,k,2)/slope(i,k,2) &
539 /slopeb(i,k,2)*qci(i,k,1)*denfac(i,k),qci(i,k,1)/dtcld)
540 endif
541 !
542 if(qci(i,k,2).gt.qcrmin) then
543 xmi = qci(i,k,2)*xnc
544 if(ifsat.ne.1) pisd(i,k) = 4.*dicon*sqrt(xmi)*den(i,k) &
545 *(rh(i,k,2)-1.)/work1(i,k,2)
546 if(pisd(i,k).lt.0.) then
547 pisd(i,k) = max(pisd(i,k),satdt)
548 pisd(i,k) = max(pisd(i,k),-qci(i,k,2)/dtcld)
549 else
550 pisd(i,k) = min(pisd(i,k),satdt)
551 endif
552 if(abs(pisd(i,k)).gt.abs(satdt)) ifsat = 1
553 endif
554 !
555 if(qrs(i,k,2).gt.qcrmin.and.ifsat.ne.1) then
556 coeres = slope2(i,k,2)*sqrt(slope(i,k,2)*slopeb(i,k,2))
557 if(ifsat.ne.1) pres(i,k,2) = (rh(i,k,2)-1.)*(precs1 &
558 /slope2(i,k,2)+precs2*work2(i,k) &
559 /coeres)/work1(i,k,2)
560 if(pres(i,k,2).lt.0.) then
561 pres(i,k,2) = max(pres(i,k,2),-qrs(i,k,2)/dtcld)
562 pres(i,k,2) = max(pres(i,k,2),satdt/2)
563 else
564 pres(i,k,2) = min(pres(i,k,2),satdt/2)
565 pres(i,k,2) = min(pres(i,k,2),qrs(i,k,2)/dtcld)
566 endif
567 if(abs(pisd(i,k)+pres(i,k,2)).ge.abs(satdt)) ifsat = 1
568 endif
569 !
570 if(supsat.gt.0.and.ifsat.ne.1) then
571 pgen(i,k) = max(0.,(xm0*xnc-max(qci(i,k,2),0.))/dtcld)
572 pgen(i,k) = min(pgen(i,k),satdt)
573 endif
574 !
575 if(qci(i,k,2).gt.qcrmin) paut(i,k,2) &
576 = max(0.,(qci(i,k,2)-xmmax*xnc)/dtcld)
577 !
578 if(t(i,k).gt.t0c) then
579 xlf = xls - xl(i,k)
580 if(qrs(i,k,2).gt.qcrmin) then
581 psml(i,k) = xka(t(i,k),den(i,k))/xlf*(t0c-t(i,k))*pi/2. &
582 *(precs1/slope2(i,k,2)+precs2*work2(i,k)/coeres)
583 psml(i,k) = min(max(psml(i,k),-qrs(i,k,2)/dtcld),0.)
584 endif
585 if(qrs(i,k,2).gt.qcrmin.and.rh(i,k,1).lt.1.) &
586 psev(i,k) = pres(i,k,2)*work1(i,k,2)/work1(i,k,1)
587 psev(i,k) = min(max(psev(i,k),-qrs(i,k,2)/dtcld),0.)
588 endif
589 enddo
590 enddo
591 !
592 !----------------------------------------------------------------
593 ! check mass conservation of generation terms and feedback to the
594 ! large scale
595 !
596 do k = kts, kte
597 do i = its, ite
598 if(t(i,k).le.t0c) then
599 !
600 ! cloud water
601 !
602 value = max(qcrmin,qci(i,k,1))
603 source = (paut(i,k,1)+pacr(i,k,1)+pacr(i,k,3))*dtcld
604 if (source.gt.value) then
605 factor = value/source
606 paut(i,k,1) = paut(i,k,1)*factor
607 pacr(i,k,1) = pacr(i,k,1)*factor
608 pacr(i,k,3) = pacr(i,k,3)*factor
609 endif
610 !
611 ! cloud ice
612 !
613 value = max(qcrmin,qci(i,k,2))
614 source = (paut(i,k,2)+pacr(i,k,2)-pgen(i,k)-pisd(i,k))*dtcld
615 if (source.gt.value) then
616 factor = value/source
617 paut(i,k,2) = paut(i,k,2)*factor
618 pacr(i,k,2) = pacr(i,k,2)*factor
619 pgen(i,k) = pgen(i,k)*factor
620 pisd(i,k) = pisd(i,k)*factor
621 endif
622 work2(i,k)=-(pres(i,k,1)+pres(i,k,2)+pgen(i,k)+pisd(i,k))
623 ! update
624 q(i,k) = q(i,k)+work2(i,k)*dtcld
625 qci(i,k,1) = max(qci(i,k,1)-(paut(i,k,1)+pacr(i,k,1) &
626 +pacr(i,k,3))*dtcld,0.)
627 qrs(i,k,1) = max(qrs(i,k,1)+(paut(i,k,1)+pacr(i,k,1) &
628 +pres(i,k,1))*dtcld,0.)
629 qci(i,k,2) = max(qci(i,k,2)-(paut(i,k,2)+pacr(i,k,2) &
630 -pgen(i,k)-pisd(i,k))*dtcld,0.)
631 qrs(i,k,2) = max(qrs(i,k,2)+(pres(i,k,2)+paut(i,k,2) &
632 +pacr(i,k,2)+pacr(i,k,3))*dtcld,0.)
633 xlf = xls-xl(i,k)
634 xlwork2 = -xls*(pres(i,k,2)+pisd(i,k)+pgen(i,k)) &
635 -xl(i,k)*pres(i,k,1) &
636 -xlf*pacr(i,k,3)
637 t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
638 else
639 !
640 ! cloud water
641 !
642 value = max(qcrmin,qci(i,k,1))
643 source=(paut(i,k,1)+pacr(i,k,1)+pacr(i,k,3))*dtcld
644 if (source.gt.value) then
645 factor = value/source
646 paut(i,k,1) = paut(i,k,1)*factor
647 pacr(i,k,1) = pacr(i,k,1)*factor
648 pacr(i,k,3) = pacr(i,k,3)*factor
649 endif
650 !
651 ! snow
652 !
653 value = max(qcrmin,qrs(i,k,2))
654 source=(-psev(i,k)-psml(i,k))*dtcld
655 if (source.gt.value) then
656 factor = value/source
657 psev(i,k) = psev(i,k)*factor
658 psml(i,k) = psml(i,k)*factor
659 endif
660 work2(i,k)=-(pres(i,k,1)+psev(i,k))
661 ! update
662 q(i,k) = q(i,k)+work2(i,k)*dtcld
663 qci(i,k,1) = max(qci(i,k,1)-(paut(i,k,1)+pacr(i,k,1) &
664 +pacr(i,k,3))*dtcld,0.)
665 qrs(i,k,1) = max(qrs(i,k,1)+(paut(i,k,1)+pacr(i,k,1) &
666 +pres(i,k,1) -psml(i,k)+pacr(i,k,3))*dtcld,0.)
667 qrs(i,k,2) = max(qrs(i,k,2)+(psml(i,k)+psev(i,k))*dtcld,0.)
668 xlf = xls-xl(i,k)
669 xlwork2 = -xlf*(psml(i,k))-xl(i,k)*(pres(i,k,1)+psev(i,k))
670 t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
671 endif
672 enddo
673 enddo
674 !
675 do k = kts, kte
676 do i = its, ite
677 qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
678 qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
679 qs(i,k,1) = max(qs(i,k,1),qmin)
680 qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
681 qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
682 qs(i,k,2) = max(qs(i,k,2),qmin)
683 enddo
684 enddo
685 !
686 !----------------------------------------------------------------
687 ! condensational/evaporational rate of cloud water if there exists
688 ! additional water vapor condensated/if evaporation of cloud water
689 ! is not enough to remove subsaturation.
690 ! use fall bariable for this process(pcon)
691 !
692 ! if(lat.eq.latd) write(iun,603)
693 do k = kts, kte
694 do i = its, ite
695 work1(i,k,1) = conden(t(i,k),q(i,k),qs(i,k,1),xl(i,k),cpm(i,k))
696 work2(i,k) = qci(i,k,1)+work1(i,k,1)
697 pcon(i,k) = min(max(work1(i,k,1)/dtcld,0.),max(q(i,k),0.)/dtcld)
698 if(qci(i,k,1).gt.qcrmin.and.work1(i,k,1).lt.0.and.t(i,k).gt.t0c) &
699 pcon(i,k) = max(work1(i,k,1),-qci(i,k,1))/dtcld
700 q(i,k) = q(i,k)-pcon(i,k)*dtcld
701 qci(i,k,1) = max(qci(i,k,1)+pcon(i,k)*dtcld,0.)
702 t(i,k) = t(i,k)+pcon(i,k)*xl(i,k)/cpm(i,k)*dtcld
703 !
704 if(lat.eq.latd.and.i.eq.lond) then
705 pgens(i,k) = pgens(i,k)+pgen(i,k)
706 pcons(i,k) = pcons(i,k)+pcon(i,k)
707 pisds(i,k) = pisds(i,k)+pisd(i,k)
708 pacrss(i,k) = pacrss(i,k)+pacr(i,k,1)+pacr(i,k,2)+pacr(i,k,3)
709 press(i,k) = press(i,k)+pres(i,k,1)+pres(i,k,2)
710 pauts(i,k) = pauts(i,k)+paut(i,k,1)+paut(i,k,2)
711 ! write(iun,604) k,p(i,k)/100., &
712 ! t(i,k)-t0c,t(i,k)-t1(i,k),q(i,k)*1000., &
713 ! (q(i,k)-q1(i,k))*1000.,rh(i,k,2)*100.,pgens(i,k)*qdt, &
714 ! pcons(i,k)*qdt,pisds(i,k)*qdt,pauts(i,k)*qdt,pacrss(i,k)*qdt, &
715 ! press(i,k)*qdt,qci(i,k,1)*1000.,qrs(i,k,1)*1000., &
716 ! qci(i,k,2)*1000.,qrs(i,k,2)*1000.
717 endif
718 enddo
719 enddo
720 !603 format(1x,' k',' p', &
721 ! ' t',' delt',' q',' delq',' rh', &
722 ! ' pgen',' pcon',' pisd',' paut',' pacr',' pres', &
723 ! ' qc',' qr',' qi',' qs')
724 !604 format(1x,i3,f6.0,4f5.1,f5.0,10f5.2)
725 !
726 !----------------------------------------------------------------
727 ! compute the fallout term:
728 ! first, vertical terminal velosity for minor loops
729 !
730 do k = kts, kte
731 do i = its, ite
732 if(qrs(i,k,1).le.qcrmin)then
733 slope(i,k,1) = lamdarmax
734 slopeb(i,k,1) = slope(i,k,1)**bvtr
735 else
736 slope(i,k,1) = lamdar(qrs(i,k,1),den(i,k))
737 slopeb(i,k,1) = slope(i,k,1)**bvtr
738 endif
739 if(qrs(i,k,2).le.qcrmin)then
740 slope(i,k,2) = lamdarmax
741 slopeb(i,k,2) = slope(i,k,2)**bvts
742 else
743 supcol = t0c-t(i,k)
744 n0sfac = min(exp(alpha*supcol),n0smax)
745 slope(i,k,2) = lamdas(qrs(i,k,2),den(i,k),n0sfac)
746 slopeb(i,k,2) = slope(i,k,2)**bvts
747 endif
748 slope2(i,k,1) = slope(i,k,1)*slope(i,k,1)
749 slope2(i,k,2) = slope(i,k,2)*slope(i,k,2)
750 enddo
751 enddo
752 !
753 do i = its, ite
754 do k = kte, kts, -1
755 work1(i,k,1) = pvtr/slopeb(i,k,1)*denfac(i,k)/delz(i,k)
756 work1(i,k,2) = pvts/slopeb(i,k,2)*denfac(i,k)/delz(i,k)
757 if(qrs(i,k,1).le.qcrmin) work1(i,k,1) = 0.
758 if(qrs(i,k,2).le.qcrmin) work1(i,k,2) = 0.
759 numdt = max(nint(max(work1(i,k,1),work1(i,k,2))*dtcld+.5),1)
760 if(qci(i,k,2).le.qmin) then
761 work2c(i,k) = 0.
762 else
763 work1c(i,k) = 3.29*(den(i,k)*qci(i,k,2))**0.16
764 work2c(i,k) = work1c(i,k)/delz(i,k)
765 endif
766 numdt = max(nint(work2c(i,k)*dtcld+.5),numdt)
767 if(numdt.ge.mstep(i)) mstep(i) = numdt
768 enddo
769 mstep(i) = min(mstep(i),mstepmax)
770 enddo
771 !
772 ! if(lat.eq.latd) write(iun,605)
773 do n = 1,mstepmax
774 k = kte
775 do i = its, ite
776 if(n.le.mstep(i)) then
777 falk(i,k,1) = den(i,k)*qrs(i,k,1)*work1(i,k,1)/mstep(i)
778 falk(i,k,2) = den(i,k)*qrs(i,k,2)*work1(i,k,2)/mstep(i)
779 fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
780 fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
781 holdrr = qrs(i,k,1)
782 holdrs = qrs(i,k,2)
783 qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcld/den(i,k),0.)
784 qrs(i,k,2) = max(qrs(i,k,2)-falk(i,k,2)*dtcld/den(i,k),0.)
785 falkc(i,k) = den(i,k)*qci(i,k,2)*work2c(i,k)/mstep(i)
786 fallc(i,k) = fallc(i,k)+falkc(i,k)
787 qci(i,k,2) = max(qci(i,k,2)-falkc(i,k)*dtcld/den(i,k),0.)
788 !
789 ! if(lat.eq.latd.and.i.eq.lond) &
790 ! write(iun,606) k,p(i,k)/100., &
791 ! t(i,k)-t0c,q(i,k)*1000.,rh(i,k,2)*100.,w(i,k),work1(i,k,1) &
792 ! *delz(i,k), work1(i,k,2)*delz(i,k),holdrr*1000.,holdrs*1000., &
793 ! qrs(i,k,1)*1000.,qrs(i,k,2)*1000.,n
794 endif
795 enddo
796 do k = kte-1, kts, -1
797 do i = its, ite
798 if(n.le.mstep(i)) then
799 falk(i,k,1) = den(i,k)*qrs(i,k,1)*work1(i,k,1)/mstep(i)
800 falk(i,k,2) = den(i,k)*qrs(i,k,2)*work1(i,k,2)/mstep(i)
801 fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
802 fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
803 holdrr = qrs(i,k,1)
804 holdrs = qrs(i,k,2)
805 qrs(i,k,1) = max(qrs(i,k,1)-(falk(i,k,1) &
806 -falk(i,k+1,1)*delz(i,k+1)/delz(i,k))*dtcld/den(i,k),0.)
807 qrs(i,k,2) = max(qrs(i,k,2)-(falk(i,k,2) &
808 -falk(i,k+1,2)*delz(i,k+1)/delz(i,k))*dtcld/den(i,k),0.)
809 falkc(i,k) = den(i,k)*qci(i,k,2)*work2c(i,k)/mstep(i)
810 fallc(i,k) = fallc(i,k)+falkc(i,k)
811 qci(i,k,2) = max(qci(i,k,2)-(falkc(i,k) &
812 -falkc(i,k+1)*delz(i,k+1)/delz(i,k))*dtcld/den(i,k),0.)
813 !
814 ! if(lat.eq.latd.and.i.eq.lond) &
815 ! write(iun,606) k,p(i,k)/100., &
816 ! t(i,k)-t0c,q(i,k)*1000.,rh(i,k,2)*100.,w(i,k),work1(i,k,1) &
817 ! *delz(i,k), work1(i,k,2)*delz(i,k),holdrr*1000.,holdrs*1000., &
818 ! qrs(i,k,1)*1000.,qrs(i,k,2)*1000.,n
819 endif
820 enddo
821 enddo
822 enddo
823 !605 format(1x,' k',' p',' t',' q',' rh',' w', &
824 ! ' vtr',' vts',' qri',' qsi',' qrf',' qsf',' mstep')
825 !606 format(1x,i3,f6.0,2f5.1,f5.0,f6.2,6f6.2,i5)
826 !
827 !----------------------------------------------------------------
828 ! rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
829 !
830 do i = its, ite
831 fallsum = fall(i,1,1)+fall(i,1,2)
832 if(fallsum.gt.0.) then
833 rainncv(i) = fallsum*delz(i,1)/denr*dtcld*1000.
834 rain(i) = fallsum*delz(i,1)/denr*dtcld*1000. &
835 + rain(i)
836 endif
837 enddo
838 !
839 ! if(lat.eq.latd) write(iun,601) latd,lond,loop,rain(lond)
840 ! 601 format(1x,' ncloud5 lat lon loop : rain(mm) ',3i6,f20.2)
841 !
842 enddo ! big loops
843
844 END SUBROUTINE ncloud52d
845 ! ...................................................................
846 REAL FUNCTION rgmma(x)
847 !-------------------------------------------------------------------
848 IMPLICIT NONE
849 !-------------------------------------------------------------------
850 ! rgmma function: use infinite product form
851 REAL :: euler
852 PARAMETER (euler=0.577215664901532)
853 REAL :: x, y
854 INTEGER :: i
855
856 if(x.eq.1.)then
857 rgmma=0.
858 else
859 rgmma=x*exp(euler*x)
860 do i=1,10000
861 y=float(i)
862 rgmma=rgmma*(1.000+x/y)*exp(-x/y)
863 enddo
864 rgmma=1./rgmma
865 endif
866 end function rgmma
867 !
868 !--------------------------------------------------------------------------
869 REAL FUNCTION fpvs(t,ice,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c)
870 !--------------------------------------------------------------------------
871 IMPLICIT NONE
872 !--------------------------------------------------------------------------
873 REAL t,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c,dldt,xa,xb,dldti, &
874 xai,xbi,ttp,tr
875 INTEGER ice
876 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
877 ttp=t0c+0.01
878 dldt=cvap-cliq
879 xa=-dldt/rv
880 xb=xa+hvap/(rv*ttp)
881 dldti=cvap-cice
882 xai=-dldti/rv
883 xbi=xai+hsub/(rv*ttp)
884 tr=ttp/t
885 if(t.lt.ttp.and.ice.eq.1) then
886 fpvs=psat*(tr**xai)*exp(xbi*(1.-tr))
887 else
888 fpvs=psat*(tr**xa)*exp(xb*(1.-tr))
889 endif
890 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
891 END FUNCTION fpvs
892
893 !-------------------------------------------------------------------
894 SUBROUTINE ncloud5init(den0,denr,dens,cl,cpv,allowed_to_read)
895 !-------------------------------------------------------------------
896 IMPLICIT NONE
897 !-------------------------------------------------------------------
898 !.... constants which may not be tunable
899
900 real, intent(in) :: den0,denr,dens,cl,cpv
901 logical, intent(in) :: allowed_to_read
902 real :: pi
903
904 pi = 4.*atan(1.)
905 xlv1 = cl-cpv
906
907 qc0 = 4./3.*pi*denr*r0**3*xncr/den0 ! 0.419e-3 -- .61e-3
908 qck1 = .104*9.8*peaut/(xncr*denr)**(1./3.)/xmyu ! 7.03
909 bvtr1 = 1.+bvtr
910 bvtr2 = 2.5+.5*bvtr
911 bvtr3 = 3.+bvtr
912 bvtr4 = 4.+bvtr
913 g1pbr = rgmma(bvtr1)
914 g3pbr = rgmma(bvtr3)
915 g4pbr = rgmma(bvtr4) ! 17.837825
916 g5pbro2 = rgmma(bvtr2) ! 1.8273
917 pvtr = avtr*g4pbr/6.
918 eacrr = 1.0
919 pacrr = pi*n0r*avtr*g3pbr*.25*eacrr
920 precr1 = 2.*pi*n0r*.78
921 precr2 = 2.*pi*n0r*.31*avtr**.5*g5pbro2
922 xm0 = (di0/dicon)**2
923 xmmax = (dimax/dicon)**2
924 !
925 bvts1 = 1.+bvts
926 bvts2 = 2.5+.5*bvts
927 bvts3 = 3.+bvts
928 bvts4 = 4.+bvts
929 g1pbs = rgmma(bvts1) !.8875
930 g3pbs = rgmma(bvts3)
931 g4pbs = rgmma(bvts4) ! 12.0786
932 g5pbso2 = rgmma(bvts2)
933 pvts = avts*g4pbs/6.
934 pacrs = pi*n0s*avts*g3pbs*.25
935 precs1 = 4.*n0s*.65
936 precs2 = 4.*n0s*.44*avts**.5*g5pbso2
937 pidn0r = pi*denr*n0r
938 pidn0s = pi*dens*n0s
939 !
940 pacrc = pi*n0s*avts*g3pbs*.25*eacrc
941 !
942 END SUBROUTINE ncloud5init
943
944 END MODULE module_mp_ncloud5
945