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