| blob | ac76ff980062cea5f3156a94620211d2ffcc8231 |
1 #if ( RWORDSIZE == 4 )
2 # define VREC vsrec
3 # define VSQRT vssqrt
4 #else
5 # define VREC vrec
6 # define VSQRT vsqrt
7 #endif
8 !
9 !Including inline expansion statistical function
10 MODULE module_mp_wdm5
11 !
12 !
13 REAL, PARAMETER, PRIVATE :: dtcldcr = 120. ! maximum time step for minor loops
14 REAL, PARAMETER, PRIVATE :: n0r = 8.e6 ! intercept parameter rain
15 REAL, PARAMETER, PRIVATE :: avtr = 841.9 ! a constant for terminal velocity of rain
16 REAL, PARAMETER, PRIVATE :: bvtr = 0.8 ! a constant for terminal velocity of rain
17 REAL, PARAMETER, PRIVATE :: r0 = .8e-5 ! 8 microm in contrast to 10 micro m
18 REAL, PARAMETER, PRIVATE :: peaut = .55 ! collection efficiency
19 REAL, PARAMETER, PRIVATE :: xncr = 3.e8 ! maritime cloud in contrast to 3.e8 in tc80
20 REAL, PARAMETER, PRIVATE :: xmyu = 1.718e-5 ! the dynamic viscosity kgm-1s-1
21 REAL, PARAMETER, PRIVATE :: avts = 11.72 ! a constant for terminal velocity of snow
22 REAL, PARAMETER, PRIVATE :: bvts = .41 ! a constant for terminal velocity of snow
23 REAL, PARAMETER, PRIVATE :: n0smax = 1.e11 ! maximum n0s (t=-90C unlimited)
24 REAL, PARAMETER, PRIVATE :: lamdacmax = 1.e10 ! limited maximum value for slope parameter of cloud water
25 REAL, PARAMETER, PRIVATE :: lamdarmax = 1.e8 ! limited maximum value for slope parameter of rain
26 REAL, PARAMETER, PRIVATE :: lamdasmax = 1.e5 ! limited maximum value for slope parameter of snow
27 REAL, PARAMETER, PRIVATE :: lamdagmax = 6.e4 ! limited maximum value for slope parameter of graupel
28 REAL, PARAMETER, PRIVATE :: dicon = 11.9 ! constant for the cloud-ice diamter
29 REAL, PARAMETER, PRIVATE :: dimax = 500.e-6 ! limited maximum value for the cloud-ice diamter
30 REAL, PARAMETER, PRIVATE :: n0s = 2.e6 ! temperature dependent intercept parameter snow
31 REAL, PARAMETER, PRIVATE :: alpha = .12 ! .122 exponen factor for n0s
32 REAL, PARAMETER, PRIVATE :: pfrz1 = 100. ! constant in Biggs freezing
33 REAL, PARAMETER, PRIVATE :: pfrz2 = 0.66 ! constant in Biggs freezing
34 REAL, PARAMETER, PRIVATE :: qcrmin = 1.e-9 ! minimun values for qr, qs, and qg
35 REAL, PARAMETER, PRIVATE :: ncmin = 1.e1 ! minimum value for Nc
36 REAL, PARAMETER, PRIVATE :: nrmin = 1.e-2 ! minimum value for Nr
37 REAL, PARAMETER, PRIVATE :: eacrc = 1.0 ! Snow/cloud-water collection efficiency
38 !
39 REAL, PARAMETER, PRIVATE :: satmax = 1.0048 ! maximum saturation value for CCN activation
40 ! 1.008 for maritime air mass /1.0048 for conti
41 REAL, PARAMETER, PRIVATE :: actk = 0.6 ! parameter for the CCN activation
42 REAL, PARAMETER, PRIVATE :: actr = 1.5 ! radius of activated CCN drops
43 REAL, PARAMETER, PRIVATE :: ncrk1 = 3.03e3 ! Long's collection kernel coefficient
44 REAL, PARAMETER, PRIVATE :: ncrk2 = 2.59e15 ! Long's collection kernel coefficient
45 REAL, PARAMETER, PRIVATE :: di100 = 1.e-4 ! parameter related with accretion and collection of cloud drops
46 REAL, PARAMETER, PRIVATE :: di600 = 6.e-4 ! parameter related with accretion and collection of cloud drops
47 REAL, PARAMETER, PRIVATE :: di2000 = 20.e-4 ! parameter related with accretion and collection of cloud drops
48 REAL, PARAMETER, PRIVATE :: di82 = 82.e-6 ! dimater related with raindrops evaporation
49 REAL, PARAMETER, PRIVATE :: di15 = 15.e-6 ! auto conversion takes place beyond this diameter
50 REAL, SAVE :: &
51 qc0, qck1,pidnc,bvtr1,bvtr2,bvtr3,bvtr4, &
52 bvtr5,bvtr7,bvtr2o5,bvtr3o5,g1pbr,g2pbr, &
53 g3pbr,g4pbr,g5pbr,g7pbr,g5pbro2,g7pbro2, &
54 pvtr,pvtrn,eacrr,pacrr, pi, &
55 precr1,precr2,xmmax,roqimax,bvts1, &
56 bvts2,bvts3,bvts4,g1pbs,g3pbs,g4pbs, &
57 g5pbso2,pvts,pacrs,precs1,precs2,pidn0r, &
58 pidn0s,pidnr,xlv1,pacrc, &
59 rslopecmax,rslopec2max,rslopec3max, &
60 rslopermax,rslopesmax,rslopegmax, &
61 rsloperbmax,rslopesbmax,rslopegbmax, &
62 rsloper2max,rslopes2max,rslopeg2max, &
63 rsloper3max,rslopes3max,rslopeg3max
64 !
65 ! Specifies code-inlining of fpvs function in WDM52D below. JM 20040507
66 !
67 CONTAINS
68 !===================================================================
69 !
70 SUBROUTINE wdm5(th, q, qc, qr, qi, qs &
71 ,nn, nc, nr &
72 ,den, pii, p, delz &
73 ,delt,g, cpd, cpv, ccn0, rd, rv, t0c &
74 ,ep1, ep2, qmin &
75 ,XLS, XLV0, XLF0, den0, denr &
76 ,cliq,cice,psat &
77 ,rain, rainncv &
78 ,snow, snowncv &
79 ,sr &
80 ,itimestep &
81 ,ids,ide, jds,jde, kds,kde &
82 ,ims,ime, jms,jme, kms,kme &
83 ,its,ite, jts,jte, kts,kte &
84 )
85 !-------------------------------------------------------------------
86 IMPLICIT NONE
87 !-------------------------------------------------------------------
88 !
89 ! This code is a WRF double-moment 5-class mixed ice
90 ! microphyiscs scheme (WDM5). The WDM microphysics scheme predicts
91 ! number concentrations for warm rain species including clouds and
92 ! rain. cloud condensation nuclei (CCN) is also predicted.
93 ! The cold rain species including ice, snow, graupel follow the
94 ! WRF single-moment 5-class microphysics (WSM5)
95 ! in which theoretical background for WSM ice phase microphysics is
96 ! based on Hong et al. (2004).
97 ! The WDM scheme is described in Lim and Hong (2009).
98 ! All units are in m.k.s. and source/sink terms in kgkg-1s-1.
99 !
100 ! WDM5 cloud scheme
101 !
102 ! Coded by Kyo-Sun Lim and Song-You Hong (Yonsei Univ.) Fall 2008
103 !
104 ! Implemented by Kyo-Sun Lim and Jimy Dudhia (NCAR) Winter 2008
105 !
106 ! Reference) Lim and Hong (LH, 2010) Mon. Wea. Rev.
107 ! Hong, Dudhia, Chen (HDC, 2004) Mon. Wea. Rev.
108 ! Hong and Lim (HL, 2006) J. Korean Meteor. Soc.
109 ! Cohard and Pinty (CP, 2000) Quart. J. Roy. Meteor. Soc.
110 ! Khairoutdinov and Kogan (KK, 2000) Mon. Wea. Rev.
111 ! Dudhia, Hong and Lim (DHL, 2008) J. Meteor. Soc. Japan
112 !
113 ! Lin, Farley, Orville (LFO, 1983) J. Appl. Meteor.
114 ! Rutledge, Hobbs (RH83, 1983) J. Atmos. Sci.
115 ! Rutledge, Hobbs (RH84, 1984) J. Atmos. Sci.
116 !
117 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
118 ims,ime, jms,jme, kms,kme , &
119 its,ite, jts,jte, kts,kte
120 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
121 INTENT(INOUT) :: &
122 th, &
123 q, &
124 qc, &
125 qi, &
126 qr, &
127 qs, &
128 nn, &
129 nc, &
130 nr
131 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
132 INTENT(IN ) :: &
133 den, &
134 pii, &
135 p, &
136 delz
137 REAL, INTENT(IN ) :: delt, &
138 g, &
139 rd, &
140 rv, &
141 t0c, &
142 den0, &
143 cpd, &
144 cpv, &
145 ccn0, &
146 ep1, &
147 ep2, &
148 qmin, &
149 XLS, &
150 XLV0, &
151 XLF0, &
152 cliq, &
153 cice, &
154 psat, &
155 denr
156 INTEGER, INTENT(IN ) :: itimestep
157 REAL, DIMENSION( ims:ime , jms:jme ), &
158 INTENT(INOUT) :: rain, &
159 rainncv, &
160 sr
161 REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
162 INTENT(INOUT) :: snow, &
163 snowncv
164 ! LOCAL VAR
165 REAL, DIMENSION( its:ite , kts:kte ) :: t
166 REAL, DIMENSION( its:ite , kts:kte, 2 ) :: qci, qrs
167 REAL, DIMENSION( its:ite , kts:kte, 3 ) :: ncr
168 CHARACTER*256 :: emess
169 INTEGER :: mkx_test
170 INTEGER :: i,j,k
171 !-------------------------------------------------------------------
172 #ifndef RUN_ON_GPU
173 IF (itimestep .eq. 1) THEN
174 DO j=jms,jme
175 DO k=kms,kme
176 DO i=ims,ime
177 nn(i,k,j) = ccn0
178 ENDDO
179 ENDDO
180 ENDDO
181 ENDIF
182 !
183 DO j=jts,jte
184 DO k=kts,kte
185 DO i=its,ite
186 t(i,k)=th(i,k,j)*pii(i,k,j)
187 qci(i,k,1) = qc(i,k,j)
188 qci(i,k,2) = qi(i,k,j)
189 qrs(i,k,1) = qr(i,k,j)
190 qrs(i,k,2) = qs(i,k,j)
191 ncr(i,k,1) = nn(i,k,j)
192 ncr(i,k,2) = nc(i,k,j)
193 ncr(i,k,3) = nr(i,k,j)
194 ENDDO
195 ENDDO
196 ! Sending array starting locations of optional variables may cause
197 ! troubles, so we explicitly change the call.
198 CALL wdm52D(t, q(ims,kms,j), qci, qrs, ncr &
199 ,den(ims,kms,j) &
200 ,p(ims,kms,j), delz(ims,kms,j) &
201 ,delt,g, cpd, cpv, ccn0, rd, rv, t0c &
202 ,ep1, ep2, qmin &
203 ,XLS, XLV0, XLF0, den0, denr &
204 ,cliq,cice,psat &
205 ,j &
206 ,rain(ims,j),rainncv(ims,j) &
207 ,sr(ims,j) &
208 ,ids,ide, jds,jde, kds,kde &
209 ,ims,ime, jms,jme, kms,kme &
210 ,its,ite, jts,jte, kts,kte &
211 ,snow(ims,j),snowncv(ims,j) &
212 )
213 DO K=kts,kte
214 DO I=its,ite
215 th(i,k,j)=t(i,k)/pii(i,k,j)
216 qc(i,k,j) = qci(i,k,1)
217 qi(i,k,j) = qci(i,k,2)
218 qr(i,k,j) = qrs(i,k,1)
219 qs(i,k,j) = qrs(i,k,2)
220 nn(i,k,j) = ncr(i,k,1)
221 nc(i,k,j) = ncr(i,k,2)
222 nr(i,k,j) = ncr(i,k,3)
223 ENDDO
224 ENDDO
225 ENDDO
226 #else
227 CALL get_wsm5_gpu_levels ( mkx_test )
228 IF ( mkx_test .LT. kte ) THEN
229 WRITE(emess,*)'Number of levels compiled for GPU WSM5 too small. ', &
230 mkx_test,' < ',kte
231 CALL wrf_error_fatal(emess)
232 ENDIF
233 CALL wsm5_host ( &
234 th(its:ite,kts:kte,jts:jte), pii(its:ite,kts:kte,jts:jte) &
235 ,q(its:ite,kts:kte,jts:jte), qc(its:ite,kts:kte,jts:jte) &
236 ,qi(its:ite,kts:kte,jts:jte), qr(its:ite,kts:kte,jts:jte) &
237 ,qs(its:ite,kts:kte,jts:jte), den(its:ite,kts:kte,jts:jte) &
238 ,p(its:ite,kts:kte,jts:jte), delz(its:ite,kts:kte,jts:jte) &
239 ,delt &
240 ,rain(its:ite,jts:jte),rainncv(its:ite,jts:jte) &
241 ,snow(its:ite,jts:jte),snowncv(its:ite,jts:jte) &
242 ,sr(its:ite,jts:jte) &
243 ,its, ite, jts, jte, kts, kte &
244 ,its, ite, jts, jte, kts, kte &
245 ,its, ite, jts, jte, kts, kte &
246 )
247 #endif
248 END SUBROUTINE wdm5
249 !===================================================================
250 !
251 SUBROUTINE wdm52D(t, q, qci, qrs, ncr, den, p, delz &
252 ,delt,g, cpd, cpv, ccn0, rd, rv, t0c &
253 ,ep1, ep2, qmin &
254 ,XLS, XLV0, XLF0, den0, denr &
255 ,cliq,cice,psat &
256 ,lat &
257 ,rain,rainncv &
258 ,sr &
259 ,ids,ide, jds,jde, kds,kde &
260 ,ims,ime, jms,jme, kms,kme &
261 ,its,ite, jts,jte, kts,kte &
262 ,snow,snowncv &
263 )
264 !-------------------------------------------------------------------
265 IMPLICIT NONE
266 !-------------------------------------------------------------------
267 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
268 ims,ime, jms,jme, kms,kme , &
269 its,ite, jts,jte, kts,kte, &
270 lat
271 REAL, DIMENSION( its:ite , kts:kte ), &
272 INTENT(INOUT) :: &
273 t
274 REAL, DIMENSION( its:ite , kts:kte, 2 ), &
275 INTENT(INOUT) :: &
276 qci, &
277 qrs
278 REAL, DIMENSION( its:ite , kts:kte, 3 ), &
279 INTENT(INOUT) :: &
280 ncr
281 REAL, DIMENSION( ims:ime , kms:kme ), &
282 INTENT(INOUT) :: &
283 q
284 REAL, DIMENSION( ims:ime , kms:kme ), &
285 INTENT(IN ) :: &
286 den, &
287 p, &
288 delz
289 REAL, INTENT(IN ) :: delt, &
290 g, &
291 cpd, &
292 cpv, &
293 ccn0, &
294 t0c, &
295 den0, &
296 rd, &
297 rv, &
298 ep1, &
299 ep2, &
300 qmin, &
301 XLS, &
302 XLV0, &
303 XLF0, &
304 cliq, &
305 cice, &
306 psat, &
307 denr
308 REAL, DIMENSION( ims:ime ), &
309 INTENT(INOUT) :: rain, &
310 rainncv, &
311 sr
312 REAL, DIMENSION( ims:ime ), OPTIONAL, &
313 INTENT(INOUT) :: snow, &
314 snowncv
315 ! LOCAL VAR
316 REAL, DIMENSION( its:ite , kts:kte , 2) :: &
317 rh, qs, rslope, rslope2, rslope3, rslopeb, &
318 falk, fall, work1, qrs_tmp
319 REAL, DIMENSION( its:ite , kts:kte ) :: &
320 rslopec, rslopec2,rslopec3
321 REAL, DIMENSION( its:ite , kts:kte, 2) :: &
322 avedia
323 REAL, DIMENSION( its:ite , kts:kte ) :: &
324 workn,falln,falkn
325 REAL, DIMENSION( its:ite , kts:kte ) :: &
326 works
327 REAL, DIMENSION( its:ite , kts:kte ) :: &
328 den_tmp, delz_tmp, ncr_tmp
329 REAL, DIMENSION( its:ite , kts:kte ) :: &
330 falkc, work1c, work2c, fallc
331 REAL, DIMENSION( its:ite , kts:kte ) :: &
332 pcact, praut, psaut, prevp, psdep, pracw, psaci, psacw, &
333 pigen, pidep, pcond, prevp_s, &
334 xl, cpm, work2, psmlt, psevp, denfac, xni, &
335 n0sfac, denqrs2, denqci
336 REAL, DIMENSION( its:ite ) :: &
337 delqrs2, delqi
338 REAL, DIMENSION( its:ite , kts:kte ) :: &
339 nraut, nracw, nrevp, ncevp, nccol, nrcol, &
340 nsacw, nseml, ncact
341 REAL :: ifac, sfac
342 !
343 #define WSM_NO_CONDITIONAL_IN_VECTOR
344 #ifdef WSM_NO_CONDITIONAL_IN_VECTOR
345 REAL, DIMENSION(its:ite) :: xal, xbl
346 #endif
347 ! variables for optimization
348 REAL, DIMENSION( its:ite ) :: tvec1
349 INTEGER, DIMENSION( its:ite ) :: mnstep, numndt
350 INTEGER, DIMENSION( its:ite ) :: mstep, numdt
351 REAL, DIMENSION(its:ite) :: rmstep
352 REAL dtcldden, rdelz, rdtcld
353 LOGICAL, DIMENSION( its:ite ) :: flgcld
354 REAL :: &
355 cpmcal, xlcal, lamdac, diffus, &
356 viscos, xka, venfac, conden, diffac, &
357 x, y, z, a, b, c, d, e, &
358 ndt, qdt, holdrr, holdrs, supcol, supcolt, pvt, &
359 coeres, supsat, dtcld, xmi, eacrs, satdt, &
360 vt2i,vt2s,acrfac, coecol, &
361 nfrzdtr, nfrzdtc, &
362 taucon, lencon, lenconcr, &
363 qimax, diameter, xni0, roqi0, &
364 fallsum, fallsum_qsi, xlwork2, factor, source, &
365 value, xlf, pfrzdtc, pfrzdtr, supice
366 REAL :: temp
367 REAL :: holdc, holdci
368 INTEGER :: i, j, k, mstepmax, &
369 iprt, latd, lond, loop, loops, ifsat, n, idim, kdim
370 ! Temporaries used for inlining fpvs function
371 REAL :: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
372 REAL :: logtr
373 !
374 !=================================================================
375 ! compute internal functions
376 !
377 cpmcal(x) = cpd*(1.-max(x,qmin))+max(x,qmin)*cpv
378 xlcal(x) = xlv0-xlv1*(x-t0c)
379 !----------------------------------------------------------------
380 ! size distributions: (x=mixing ratio, y=air density):
381 ! valid for mixing ratio > 1.e-9 kg/kg.
382 !
383 ! Optimizatin : A**B => exp(log(A)*(B))
384 lamdac(x,y,z)= exp(log(((pidnc*z)/(x*y)))*((.33333333)))
385 !
386 !----------------------------------------------------------------
387 ! diffus: diffusion coefficient of the water vapor
388 ! viscos: kinematic viscosity(m2s-1)
389 ! diffus(x,y) = 8.794e-5 * exp(log(x)*(1.81)) / y
390 ! viscos(x,y) = 1.496e-6 * (x*sqrt(x)) /(x+120.)/y
391 ! xka(x,y) = 1.414e3*viscos(x,y)*y
392 ! diffac(a,b,c,d,e) = d*a*a/(xka(c,d)*rv*c*c)+1./(e*diffus(c,b))
393 ! venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333))) &
394 ! /sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
395 ! conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+d*d/(rv*e)*c/(a*a))
396 !
397 !
398 idim = ite-its+1
399 kdim = kte-kts+1
400 !
401 !----------------------------------------------------------------
402 ! paddint 0 for negative values generated by dynamics
403 !
404 do k = kts, kte
405 do i = its, ite
406 qci(i,k,1) = max(qci(i,k,1),0.0)
407 qrs(i,k,1) = max(qrs(i,k,1),0.0)
408 qci(i,k,2) = max(qci(i,k,2),0.0)
409 qrs(i,k,2) = max(qrs(i,k,2),0.0)
410 ncr(i,k,1) = max(ncr(i,k,1),0.)
411 ncr(i,k,2) = max(ncr(i,k,2),0.)
412 ncr(i,k,3) = max(ncr(i,k,3),0.)
413 enddo
414 enddo
415 !
416 ! latent heat for phase changes and heat capacity. neglect the
417 ! changes during microphysical process calculation
418 ! emanuel(1994)
419 !
420 do k = kts, kte
421 do i = its, ite
422 cpm(i,k) = cpmcal(q(i,k))
423 xl(i,k) = xlcal(t(i,k))
424 delz_tmp(i,k) = delz(i,k)
425 den_tmp(i,k) = den(i,k)
426 enddo
427 enddo
428 !
429 !----------------------------------------------------------------
430 ! compute the minor time steps.
431 !
432 loops = max(nint(delt/dtcldcr),1)
433 dtcld = delt/loops
434 if(delt.le.dtcldcr) dtcld = delt
435 !
436 do loop = 1,loops
437 !
438 !----------------------------------------------------------------
439 ! initialize the large scale variables
440 !
441 do i = its, ite
442 mstep(i) = 1
443 mnstep(i) = 1
444 flgcld(i) = .true.
445 enddo
446 !
447 ! do k = kts, kte
448 ! do i = its, ite
449 ! denfac(i,k) = sqrt(den0/den(i,k))
450 ! enddo
451 ! enddo
452 do k = kts, kte
453 CALL VREC( tvec1(its), den(its,k), ite-its+1)
454 do i = its, ite
455 tvec1(i) = tvec1(i)*den0
456 enddo
457 CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
458 enddo
459 !
460 ! Inline expansion for fpvs
461 ! qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
462 ! qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
463 hsub = xls
464 hvap = xlv0
465 cvap = cpv
466 ttp=t0c+0.01
467 dldt=cvap-cliq
468 xa=-dldt/rv
469 xb=xa+hvap/(rv*ttp)
470 dldti=cvap-cice
471 xai=-dldti/rv
472 xbi=xai+hsub/(rv*ttp)
473 ! this is for compilers where the conditional inhibits vectorization
474 #ifdef WSM_NO_CONDITIONAL_IN_VECTOR
475 do k = kts, kte
476 do i = its, ite
477 if(t(i,k).lt.ttp) then
478 xal(i) = xai
479 xbl(i) = xbi
480 else
481 xal(i) = xa
482 xbl(i) = xb
483 endif
484 enddo
485 do i = its, ite
486 tr=ttp/t(i,k)
487 logtr=log(tr)
488 qs(i,k,1)=psat*exp(logtr*(xa)+xb*(1.-tr))
489 qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
490 qs(i,k,1) = max(qs(i,k,1),qmin)
491 rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
492 qs(i,k,2)=psat*exp(logtr*(xal(i))+xbl(i)*(1.-tr))
493 qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
494 qs(i,k,2) = max(qs(i,k,2),qmin)
495 rh(i,k,2) = max(q(i,k) / qs(i,k,2),qmin)
496 enddo
497 enddo
498 #else
499 do k = kts, kte
500 do i = its, ite
501 tr=ttp/t(i,k)
502 logtr=log(tr)
503 qs(i,k,1)=psat*exp(logtr*(xa)+xb*(1.-tr))
504 qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
505 qs(i,k,1) = max(qs(i,k,1),qmin)
506 rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
507 if(t(i,k).lt.ttp) then
508 qs(i,k,2)=psat*exp(logtr*(xai)+xbi*(1.-tr))
509 else
510 qs(i,k,2)=psat*exp(logtr*(xa)+xb*(1.-tr))
511 endif
512 qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
513 qs(i,k,2) = max(qs(i,k,2),qmin)
514 rh(i,k,2) = max(q(i,k) / qs(i,k,2),qmin)
515 enddo
516 enddo
517 #endif
518 !
519 !----------------------------------------------------------------
520 ! initialize the variables for microphysical physics
521 !
522 !
523 do k = kts, kte
524 do i = its, ite
525 prevp(i,k) = 0.
526 psdep(i,k) = 0.
527 praut(i,k) = 0.
528 psaut(i,k) = 0.
529 pracw(i,k) = 0.
530 psaci(i,k) = 0.
531 psacw(i,k) = 0.
532 pigen(i,k) = 0.
533 pidep(i,k) = 0.
534 pcond(i,k) = 0.
535 psmlt(i,k) = 0.
536 psevp(i,k) = 0.
537 pcact(i,k) = 0.
538 prevp_s(i,k) = 0.
539 falk(i,k,1) = 0.
540 falk(i,k,2) = 0.
541 fall(i,k,1) = 0.
542 fall(i,k,2) = 0.
543 fallc(i,k) = 0.
544 falkc(i,k) = 0.
545 falln(i,k) = 0.
546 falkn(i,k) = 0.
547 xni(i,k) = 1.e3
548 nsacw(i,k) = 0.
549 nseml(i,k) = 0.
550 nracw(i,k) = 0.
551 nccol(i,k) = 0.
552 nrcol(i,k) = 0.
553 ncact(i,k) = 0.
554 nraut(i,k) = 0.
555 nrevp(i,k) = 0.
556 ncevp(i,k) = 0.
557 enddo
558 enddo
559 !
560 !----------------------------------------------------------------
561 ! compute the fallout term:
562 ! first, vertical terminal velosity for minor loops
563 !
564 do k = kts, kte
565 do i = its, ite
566 if(qci(i,k,1).le.qmin .or. ncr(i,k,2).le.ncmin)then
567 rslopec(i,k) = rslopecmax
568 rslopec2(i,k) = rslopec2max
569 rslopec3(i,k) = rslopec3max
570 else
571 rslopec(i,k) = 1./lamdac(qci(i,k,1),den(i,k),ncr(i,k,2))
572 rslopec2(i,k) = rslopec(i,k)*rslopec(i,k)
573 rslopec3(i,k) = rslopec2(i,k)*rslopec(i,k)
574 endif
575 !-------------------------------------------------------------
576 ! Ni: ice crystal number concentraiton [HDC 5c]
577 !-------------------------------------------------------------
578 ! xni(i,k) = min(max(5.38e7*(den(i,k) &
579 ! *max(qci(i,k,2),qmin))**0.75,1.e3),1.e6)
580 temp = (den(i,k)*max(qci(i,k,2),qmin))
581 temp = sqrt(sqrt(temp*temp*temp))
582 xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
583 enddo
584 enddo
585 do k = kts, kte
586 do i = its, ite
587 qrs_tmp(i,k,1) = qrs(i,k,1)
588 qrs_tmp(i,k,2) = qrs(i,k,2)
589 ncr_tmp(i,k) = ncr(i,k,3)
590 enddo
591 enddo
592 call slope_wdm5(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
593 rslope3,work1,workn,its,ite,kts,kte)
594 !----------------------------------------------------------------
595 ! compute the fallout term:
596 ! first, vertical terminal velosity for minor loops
597 !----------------------------------------------------------------
598 !
599 ! vt update for qr and nr
600 mstepmax = 1
601 numdt = 1
602 do k = kte, kts, -1
603 do i = its, ite
604 work1(i,k,1) = work1(i,k,1)/delz(i,k)
605 workn(i,k) = workn(i,k)/delz(i,k)
606 numdt(i) = max(nint(max(work1(i,k,1),workn(i,k))*dtcld+.5),1)
607 if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
608 enddo
609 enddo
610 do i = its, ite
611 if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
612 enddo
613 !
614 do n = 1, mstepmax
615 k = kte
616 do i = its, ite
617 if(n.le.mstep(i)) then
618 falk(i,k,1) = den(i,k)*qrs(i,k,1)*work1(i,k,1)/mstep(i)
619 falkn(i,k) = ncr(i,k,3)*workn(i,k)/mstep(i)
620 fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
621 falln(i,k) = falln(i,k)+falkn(i,k)
622 qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcld/den(i,k),0.)
623 ncr(i,k,3) = max(ncr(i,k,3)-falkn(i,k)*dtcld,0.)
624 endif
625 enddo
626 do k = kte-1, kts, -1
627 do i = its, ite
628 if(n.le.mstep(i)) then
629 falk(i,k,1) = den(i,k)*qrs(i,k,1)*work1(i,k,1)/mstep(i)
630 falkn(i,k) = ncr(i,k,3)*workn(i,k)/mstep(i)
631 fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
632 falln(i,k) = falln(i,k)+falkn(i,k)
633 qrs(i,k,1) = max(qrs(i,k,1)-(falk(i,k,1)-falk(i,k+1,1) &
634 *delz(i,k+1)/delz(i,k))*dtcld/den(i,k),0.)
635 ncr(i,k,3) = max(ncr(i,k,3)-(falkn(i,k)-falkn(i,k+1)*delz(i,k+1) &
636 /delz(i,k))*dtcld,0.)
637 endif
638 enddo
639 enddo
640 do k = kts, kte
641 do i = its, ite
642 qrs_tmp(i,k,1) = qrs(i,k,1)
643 ncr_tmp(i,k) = ncr(i,k,3)
644 enddo
645 enddo
646 call slope_rain(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
647 rslope3,work1,workn,its,ite,kts,kte)
648 do k = kte, kts, -1
649 do i = its, ite
650 work1(i,k,1) = work1(i,k,1)/delz(i,k)
651 workn(i,k) = workn(i,k)/delz(i,k)
652 enddo
653 enddo
654 enddo
655 ! for semi
656 do k = kte, kts, -1
657 do i = its, ite
658 works(i,k) = work1(i,k,2)
659 denqrs2(i,k) = den(i,k)*qrs(i,k,2)
660 if(qrs(i,k,2).le.0.0) works(i,k) = 0.0
661 enddo
662 enddo
663 call nislfv_rain_plm(idim,kdim,den_tmp,denfac,t,delz_tmp,works,denqrs2, &
664 delqrs2,dtcld,2,1)
665 do k = kts, kte
666 do i = its, ite
667 qrs(i,k,2) = max(denqrs2(i,k)/den(i,k),0.)
668 fall(i,k,2) = denqrs2(i,k)*works(i,k)/delz(i,k)
669 enddo
670 enddo
671 do i = its, ite
672 fall(i,1,2) = delqrs2(i)/delz(i,1)/dtcld
673 enddo
674 do k = kts, kte
675 do i = its, ite
676 qrs_tmp(i,k,1) = qrs(i,k,1)
677 qrs_tmp(i,k,2) = qrs(i,k,2)
678 ncr_tmp(i,k) = ncr(i,k,3)
679 enddo
680 enddo
681 call slope_wdm5(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
682 rslope3,work1,workn,its,ite,kts,kte)
683 !
684 do k = kte, kts, -1
685 do i = its, ite
686 supcol = t0c-t(i,k)
687 n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
688 if(t(i,k).gt.t0c.and.qrs(i,k,2).gt.0.) then
689 !----------------------------------------------------------------
690 ! psmlt: melting of snow [HL A33] [RH83 A25]
691 ! (T>T0: QS->QR)
692 !----------------------------------------------------------------
693 xlf = xlf0
694 ! work2(i,k)= venfac(p(i,k),t(i,k),den(i,k))
695 work2(i,k)= (exp(log(((1.496e-6*((t(i,k))*sqrt(t(i,k))) &
696 /((t(i,k))+120.)/(den(i,k)))/(8.794e-5 &
697 *exp(log(t(i,k))*(1.81))/p(i,k)))) &
698 *((.3333333)))/sqrt((1.496e-6*((t(i,k)) &
699 *sqrt(t(i,k)))/((t(i,k))+120.)/(den(i,k)))) &
700 *sqrt(sqrt(den0/(den(i,k)))))
701 coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
702 ! psmlt(i,k) = xka(t(i,k),den(i,k))/xlf*(t0c-t(i,k))*pi/2. &
703 ! *n0sfac(i,k)*(precs1*rslope2(i,k,2)+precs2 &
704 ! *work2(i,k)*coeres)
705 psmlt(i,k) = (1.414e3*(1.496e-6 * ((t(i,k))*sqrt(t(i,k))) &
706 /((t(i,k))+120.)/(den(i,k)))*(den(i,k)))/xlf &
707 *(t0c-t(i,k))*pi/2.*n0sfac(i,k) &
708 *(precs1*rslope2(i,k,2)+precs2*work2(i,k)*coeres)
709 psmlt(i,k) = min(max(psmlt(i,k)*dtcld/mstep(i),-qrs(i,k,2) &
710 /mstep(i)),0.)
711 !-------------------------------------------------------------------
712 ! nsmlt: melgin of snow [LH A27]
713 ! (T>T0: ->NR)
714 !-------------------------------------------------------------------
715 if(qrs(i,k,2).gt.qcrmin) then
716 sfac = rslope(i,k,2)*n0s*n0sfac(i,k)*mstep(i)/qrs(i,k,2)
717 ncr(i,k,3) = ncr(i,k,3) - sfac*psmlt(i,k)
718 endif
719 qrs(i,k,2) = qrs(i,k,2) + psmlt(i,k)
720 qrs(i,k,1) = qrs(i,k,1) - psmlt(i,k)
721 t(i,k) = t(i,k) + xlf/cpm(i,k)*psmlt(i,k)
722 endif
723 enddo
724 enddo
725 !---------------------------------------------------------------
726 ! Vice [ms-1] : fallout of ice crystal [HDC 5a]
727 !---------------------------------------------------------------
728 do k = kte, kts, -1
729 do i = its, ite
730 if(qci(i,k,2).le.0.) then
731 work1c(i,k) = 0.
732 else
733 xmi = den(i,k)*qci(i,k,2)/xni(i,k)
734 diameter = max(min(dicon * sqrt(xmi),dimax), 1.e-25)
735 work1c(i,k) = 1.49e4*exp(log(diameter)*(1.31))
736 endif
737 enddo
738 enddo
739 !
740 ! forward semi-laglangian scheme (JH), PCM (piecewise constant), (linear)
741 !
742 do k = kte, kts, -1
743 do i = its, ite
744 denqci(i,k) = den(i,k)*qci(i,k,2)
745 enddo
746 enddo
747 call nislfv_rain_plm(idim,kdim,den_tmp,denfac,t,delz_tmp,work1c,denqci, &
748 delqi,dtcld,1,0)
749 do k = kts, kte
750 do i = its, ite
751 qci(i,k,2) = max(denqci(i,k)/den(i,k),0.)
752 enddo
753 enddo
754 do i = its, ite
755 fallc(i,1) = delqi(i)/delz(i,1)/dtcld
756 enddo
757 !
758 !----------------------------------------------------------------
759 ! rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
760 !
761 do i = its, ite
762 fallsum = fall(i,1,1)+fall(i,1,2)+fallc(i,1)
763 fallsum_qsi = fall(i,1,2)+fallc(i,1)
764 rainncv(i) = 0.
765 if(fallsum.gt.0.) then
766 rainncv(i) = fallsum*delz(i,1)/denr*dtcld*1000.
767 rain(i) = fallsum*delz(i,1)/denr*dtcld*1000.+rain(i)
768 endif
769 if (PRESENT (snowncv) .and. PRESENT (snow)) then
770 snowncv(i) = 0.
771 if(fallsum_qsi.gt.0.) then
772 snowncv(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000.
773 snow(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000.+snow(i)
774 endif
775 endif
776 sr(i) = 0.
777 if(fallsum.gt.0.)sr(i)=fallsum_qsi*delz(i,kts)/denr*dtcld*1000. &
778 /(rainncv(i)+1.e-12)
779 enddo
780 !
781 !---------------------------------------------------------------
782 ! pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
783 ! (T>T0: QI->QC)
784 !---------------------------------------------------------------
785 do k = kts, kte
786 do i = its, ite
787 supcol = t0c-t(i,k)
788 xlf = xls-xl(i,k)
789 if(supcol.lt.0.) xlf = xlf0
790 if(supcol.lt.0 .and. qci(i,k,2).gt.0.) then
791 qci(i,k,1) = qci(i,k,1)+qci(i,k,2)
792 !---------------------------------------------------------------
793 ! nimlt: instantaneous melting of cloud ice [LH A18]
794 ! (T>T0: ->NC)
795 !--------------------------------------------------------------
796 ncr(i,k,2) = ncr(i,k,2) + xni(i,k)
797 t(i,k) = t(i,k) - xlf/cpm(i,k)*qci(i,k,2)
798 qci(i,k,2) = 0.
799 endif
800 !---------------------------------------------------------------
801 ! pihmf: homogeneous freezing of cloud water below -40c [HL A45]
802 ! (T<-40C: QC->QI)
803 !---------------------------------------------------------------
804 if(supcol.gt.40. .and. qci(i,k,1).gt.0.) then
805 qci(i,k,2) = qci(i,k,2) + qci(i,k,1)
806 !---------------------------------------------------------------
807 ! nihmf: homogeneous of cloud water below -40c [LH A17]
808 ! (T<-40C: NC->)
809 !---------------------------------------------------------------
810 if(ncr(i,k,2).gt.0.) ncr(i,k,2) = 0.
811 t(i,k) = t(i,k) + xlf/cpm(i,k)*qci(i,k,1)
812 qci(i,k,1) = 0.
813 endif
814 !---------------------------------------------------------------
815 ! pihtf: heterogeneous freezing of cloud water [HL A44]
816 ! (T0>T>-40C: QC->QI)
817 !---------------------------------------------------------------
818 if(supcol.gt.0. .and. qci(i,k,1).gt.0.) then
819 supcolt=min(supcol,70.)
820 pfrzdtc = min(pi*pi*pfrz1*(exp(pfrz2*supcolt)-1.)*denr/den(i,k) &
821 *ncr(i,k,2)*rslopec3(i,k)*rslopec3(i,k)/18.*dtcld,qci(i,k,1))
822 !---------------------------------------------------------------
823 ! nihtf: heterogeneous of cloud water [LH A16]
824 ! (T0>T>-40C: NC->)
825 !---------------------------------------------------------------
826 if(ncr(i,k,2).gt.ncmin) then
827 nfrzdtc = min(pi*pfrz1*(exp(pfrz2*supcolt)-1.)*ncr(i,k,2) &
828 *rslopec3(i,k)/6.*dtcld,ncr(i,k,2))
829 ncr(i,k,2) = ncr(i,k,2) - nfrzdtc
830 endif
831 qci(i,k,2) = qci(i,k,2) + pfrzdtc
832 t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtc
833 qci(i,k,1) = qci(i,k,1)-pfrzdtc
834 endif
835 !---------------------------------------------------------------
836 ! psfrz: freezing of rain water [HL A20] [LFO 45]
837 ! (T<T0, QR->QS)
838 !---------------------------------------------------------------
839 if(supcol.gt.0. .and. qrs(i,k,1).gt.0.) then
840 supcolt=min(supcol,70.)
841 pfrzdtr = min(140.*(pi*pi)*pfrz1*ncr(i,k,3)*denr/den(i,k) &
842 *(exp(pfrz2*supcolt)-1.)*rslope3(i,k,1)*rslope3(i,k,1) &
843 *dtcld,qrs(i,k,1))
844 !---------------------------------------------------------------
845 ! nsfrz: freezing of rain water [LH A26]
846 ! (T<T0, NR-> )
847 !---------------------------------------------------------------
848 if(ncr(i,k,3).gt.nrmin) then
849 nfrzdtr = min(4.*pi*pfrz1*ncr(i,k,3)*(exp(pfrz2*supcolt)-1.) &
850 *rslope3(i,k,1)*dtcld,ncr(i,k,3))
851 ncr(i,k,3) = ncr(i,k,3)-nfrzdtr
852 endif
853 qrs(i,k,2) = qrs(i,k,2) + pfrzdtr
854 t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtr
855 qrs(i,k,1) = qrs(i,k,1)-pfrzdtr
856 endif
857 enddo
858 enddo
859 !
860 do k = kts, kte
861 do i = its, ite
862 ncr(i,k,2) = max(ncr(i,k,2),0.0)
863 ncr(i,k,3) = max(ncr(i,k,3),0.0)
864 enddo
865 enddo
866 !----------------------------------------------------------------
867 ! update the slope parameters for microphysics computation
868 !
869 do k = kts, kte
870 do i = its, ite
871 qrs_tmp(i,k,1) = qrs(i,k,1)
872 qrs_tmp(i,k,2) = qrs(i,k,2)
873 ncr_tmp(i,k) = ncr(i,k,3)
874 enddo
875 enddo
876 call slope_wdm5(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
877 rslope3,work1,workn,its,ite,kts,kte)
878 do k = kts, kte
879 do i = its, ite
880 !-----------------------------------------------------------------
881 ! compute the mean-volume drop diameter [LH A10]
882 ! for raindrop distribution
883 !-----------------------------------------------------------------
884 avedia(i,k,2) = rslope(i,k,1)*((24.)**(.3333333))
885 if(qci(i,k,1).le.qmin .or. ncr(i,k,2).le.ncmin) then
886 rslopec(i,k) = rslopecmax
887 rslopec2(i,k) = rslopec2max
888 rslopec3(i,k) = rslopec3max
889 else
890 rslopec(i,k) = 1./lamdac(qci(i,k,1),den(i,k),ncr(i,k,2))
891 rslopec2(i,k) = rslopec(i,k)*rslopec(i,k)
892 rslopec3(i,k) = rslopec2(i,k)*rslopec(i,k)
893 endif
894 !--------------------------------------------------------------------
895 ! compute the mean-volume drop diameter [LH A7]
896 ! for cloud-droplet distribution
897 !--------------------------------------------------------------------
898 avedia(i,k,1) = rslopec(i,k)
899 enddo
900 enddo
901 !----------------------------------------------------------------
902 ! work1: the thermodynamic term in the denominator associated with
903 ! heat conduction and vapor diffusion
904 ! (ry88, y93, h85)
905 ! work2: parameter associated with the ventilation effects(y93)
906 !
907 do k = kts, kte
908 do i = its, ite
909 ! work1(i,k,1) = diffac(xl(i,k),p(i,k),t(i,k),den(i,k),qs(i,k,1))
910 work1(i,k,1) = ((((den(i,k))*(xl(i,k))*(xl(i,k)))*((t(i,k))+120.) &
911 *(den(i,k)))/(1.414e3*(1.496e-6*((t(i,k))*sqrt(t(i,k))))&
912 *(den(i,k))*(rv*(t(i,k))*(t(i,k))))) &
913 + p(i,k)/((qs(i,k,1))*(8.794e-5*exp(log(t(i,k))*(1.81))))
914 ! work1(i,k,2) = diffac(xls,p(i,k),t(i,k),den(i,k),qs(i,k,2))
915 work1(i,k,2) = ((((den(i,k))*(xls)*(xls))*((t(i,k))+120.)*(den(i,k)))&
916 /(1.414e3*(1.496e-6*((t(i,k))*sqrt(t(i,k))))*(den(i,k)) &
917 *(rv*(t(i,k))*(t(i,k)))) &
918 + p(i,k)/(qs(i,k,2)*(8.794e-5*exp(log(t(i,k))*(1.81)))))
919 ! work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
920 work2(i,k) = (exp(.3333333*log(((1.496e-6 * ((t(i,k))*sqrt(t(i,k)))) &
921 *p(i,k))/(((t(i,k))+120.)*den(i,k)*(8.794e-5 &
922 *exp(log(t(i,k))*(1.81))))))*sqrt(sqrt(den0/(den(i,k))))) &
923 /sqrt((1.496e-6*((t(i,k))*sqrt(t(i,k)))) &
924 /(((t(i,k))+120.)*den(i,k)))
925 enddo
926 enddo
927 !
928 !===============================================================
929 !
930 ! warm rain processes
931 !
932 ! - follows the processes in RH83 and LFO except for autoconcersion
933 !
934 !===============================================================
935 !
936 do k = kts, kte
937 do i = its, ite
938 supsat = max(q(i,k),qmin)-qs(i,k,1)
939 satdt = supsat/dtcld
940 !---------------------------------------------------------------
941 ! praut: auto conversion rate from cloud to rain [LH 9] [CP 17]
942 ! (QC->QR)
943 !---------------------------------------------------------------
944 lencon = 2.7e-2*den(i,k)*qci(i,k,1)*(1.e20/16.*rslopec2(i,k) &
945 *rslopec2(i,k)-0.4)
946 lenconcr = max(1.2*lencon,qcrmin)
947 if(avedia(i,k,1).gt.di15) then
948 taucon = 3.7/den(i,k)/qci(i,k,1)/(0.5e6*rslopec(i,k)-7.5)
949 praut(i,k) = lencon/taucon
950 praut(i,k) = min(max(praut(i,k),0.),qci(i,k,1)/dtcld)
951 !---------------------------------------------------------------
952 ! nraut: auto conversion rate from cloud to rain [LH A6][CP 18 & 19]
953 ! (NC->NR)
954 !---------------------------------------------------------------
955 nraut(i,k) = 3.5e9*den(i,k)*praut(i,k)
956 if(qrs(i,k,1).gt.lenconcr) &
957 nraut(i,k) = ncr(i,k,3)/qrs(i,k,1)*praut(i,k)
958 nraut(i,k) = min(nraut(i,k),ncr(i,k,2)/dtcld)
959 endif
960 !---------------------------------------------------------------
961 ! pracw: accretion of cloud water by rain [LH 10][CP 22 & 23]
962 ! (QC->QR)
963 ! nracw: accretion of cloud water by rain [LH A9]
964 ! (NC->)
965 !---------------------------------------------------------------
966 if(qrs(i,k,1).ge.lenconcr) then
967 if(avedia(i,k,2).ge.di100) then
968 nracw(i,k) = min(ncrk1*ncr(i,k,2)*ncr(i,k,3)*(rslopec3(i,k) &
969 + 24.*rslope3(i,k,1)),ncr(i,k,2)/dtcld)
970 pracw(i,k) = min(pi/6.*(denr/den(i,k))*ncrk1*ncr(i,k,2) &
971 *ncr(i,k,3)*rslopec3(i,k)*(2.*rslopec3(i,k) &
972 + 24.*rslope3(i,k,1)),qci(i,k,1)/dtcld)
973 else
974 nracw(i,k) = min(ncrk2*ncr(i,k,2)*ncr(i,k,3)*(2.*rslopec3(i,k) &
975 *rslopec3(i,k)+5040.*rslope3(i,k,1) &
976 *rslope3(i,k,1)),ncr(i,k,2)/dtcld)
977 pracw(i,k) = min(pi/6.*(denr/den(i,k))*ncrk2*ncr(i,k,2) &
978 *ncr(i,k,3)*rslopec3(i,k)*(6.*rslopec3(i,k) &
979 *rslopec3(i,k)+5040.*rslope3(i,k,1) &
980 *rslope3(i,k,1)),qci(i,k,1)/dtcld)
981 endif
982 endif
983 !----------------------------------------------------------------
984 ! nccol: self collection of cloud water [LH A8][CP 24 & 25]
985 ! (NC->)
986 !----------------------------------------------------------------
987 if(avedia(i,k,1).ge.di100) then
988 nccol(i,k) = ncrk1*ncr(i,k,2)*ncr(i,k,2)*rslopec3(i,k)
989 else
990 nccol(i,k) = 2.*ncrk2*ncr(i,k,2)*ncr(i,k,2)*rslopec3(i,k) &
991 *rslopec3(i,k)
992 endif
993 !----------------------------------------------------------------
994 ! nrcol: self collection of rain-drops and break-up [LH A21][CP 24 & 25]
995 ! (NR->)
996 !----------------------------------------------------------------
997 if(qrs(i,k,1).ge.lenconcr) then
998 if(avedia(i,k,2).lt.di100) then
999 nrcol(i,k) = 5040.*ncrk2*ncr(i,k,3)*ncr(i,k,3)*rslope3(i,k,1) &
1000 *rslope3(i,k,1)
1001 elseif(avedia(i,k,2).ge.di100 .and. avedia(i,k,2).lt.di600) then
1002 nrcol(i,k) = 24.*ncrk1*ncr(i,k,3)*ncr(i,k,3)*rslope3(i,k,1)
1003 elseif(avedia(i,k,2).ge.di600 .and. avedia(i,k,2).lt.di2000) then
1004 coecol = -2.5e3*(avedia(i,k,2)-di600)
1005 nrcol(i,k) = 24.*exp(coecol)*ncrk1*ncr(i,k,3)*ncr(i,k,3) &
1006 *rslope3(i,k,1)
1007 else
1008 nrcol(i,k) = 0.
1009 endif
1010 endif
1011 !---------------------------------------------------------------
1012 ! prevp: evaporation/condensation rate of rain [HL A41]
1013 ! (QV->QR or QR->QV)
1014 !---------------------------------------------------------------
1015 if(qrs(i,k,1).gt.0.) then
1016 coeres = rslope(i,k,1)*sqrt(rslope(i,k,1)*rslopeb(i,k,1))
1017 prevp(i,k) = (rh(i,k,1)-1.)*ncr(i,k,3)*(precr1*rslope(i,k,1) &
1018 +precr2*work2(i,k)*coeres)/work1(i,k,1)
1019 if(prevp(i,k).lt.0.) then
1020 prevp(i,k) = max(prevp(i,k),-qrs(i,k,1)/dtcld)
1021 prevp(i,k) = max(prevp(i,k),satdt/2)
1022 !----------------------------------------------------------------
1023 ! Nrevp: evaporation/condensation rate of rain [LH A14]
1024 ! (NR->NC)
1025 !----------------------------------------------------------------
1026 if(avedia(i,k,2).le.di82) then
1027 nrevp(i,k) = ncr(i,k,3)/dtcld
1028 !----------------------------------------------------------------
1029 ! Prevp_s: evaporation/condensation rate of rain [LH A15] [KK 23]
1030 ! (QR->QC)
1031 !----------------------------------------------------------------
1032 prevp_s(i,k) = qrs(i,k,1)/dtcld
1033 endif
1034 else
1035 !
1036 prevp(i,k) = min(prevp(i,k),satdt/2)
1037 endif
1038 endif
1039 enddo
1040 enddo
1041 !
1042 !===============================================================
1043 !
1044 ! cold rain processes
1045 !
1046 ! - follows the revised ice microphysics processes in HDC
1047 ! - the processes same as in RH83 and RH84 and LFO behave
1048 ! following ice crystal hapits defined in HDC, inclduing
1049 ! intercept parameter for snow (n0s), ice crystal number
1050 ! concentration (ni), ice nuclei number concentration
1051 ! (n0i), ice diameter (d)
1052 !
1053 !===============================================================
1054 !
1055 rdtcld = 1./dtcld
1056 do k = kts, kte
1057 do i = its, ite
1058 supcol = t0c-t(i,k)
1059 n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
1060 supsat = max(q(i,k),qmin)-qs(i,k,2)
1061 satdt = supsat/dtcld
1062 ifsat = 0
1063 !-------------------------------------------------------------
1064 ! Ni: ice crystal number concentraiton [HDC 5c]
1065 !-------------------------------------------------------------
1066 ! xni(i,k) = min(max(5.38e7*(den(i,k) &
1067 ! *max(qci(i,k,2),qmin))**0.75,1.e3),1.e6)
1068 temp = (den(i,k)*max(qci(i,k,2),qmin))
1069 temp = sqrt(sqrt(temp*temp*temp))
1070 xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
1071 eacrs = exp(0.07*(-supcol))
1072 !
1073 if(supcol.gt.0) then
1074 if(qrs(i,k,2).gt.qcrmin .and. qci(i,k,2).gt.qmin) then
1075 xmi = den(i,k)*qci(i,k,2)/xni(i,k)
1076 diameter = min(dicon * sqrt(xmi),dimax)
1077 vt2i = 1.49e4*diameter**1.31
1078 vt2s = pvts*rslopeb(i,k,2)*denfac(i,k)
1079 !-------------------------------------------------------------
1080 ! psaci: Accretion of cloud ice by rain [HDC 10]
1081 ! (T<T0: QI->QS)
1082 !-------------------------------------------------------------
1083 acrfac = 2.*rslope3(i,k,2)+2.*diameter*rslope2(i,k,2) &
1084 + diameter**2*rslope(i,k,2)
1085 psaci(i,k) = pi*qci(i,k,2)*eacrs*n0s*n0sfac(i,k)*abs(vt2s-vt2i) &
1086 *acrfac/4.
1087 endif
1088 endif
1089 !-------------------------------------------------------------
1090 ! psacw: Accretion of cloud water by snow [HL A7] [LFO 24]
1091 ! (T<T0: QC->QS, and T>=T0: QC->QR)
1092 !-------------------------------------------------------------
1093 if(qrs(i,k,2).gt.qcrmin .and. qci(i,k,1).gt.qmin) then
1094 psacw(i,k) = min(pacrc*n0sfac(i,k)*rslope3(i,k,2)*rslopeb(i,k,2) &
1095 *qci(i,k,1)*denfac(i,k),qci(i,k,1)*rdtcld)
1096 endif
1097 !-------------------------------------------------------------
1098 ! nsacw: Accretion of cloud water by snow [LH A12]
1099 ! (NC ->)
1100 !-------------------------------------------------------------
1101 if(qrs(i,k,2).gt.qcrmin .and. ncr(i,k,2).gt.ncmin) then
1102 nsacw(i,k) = min(pacrc*n0sfac(i,k)*rslope3(i,k,2)*rslopeb(i,k,2) &
1103 *ncr(i,k,2)*denfac(i,k),ncr(i,k,2)/dtcld)
1104 endif
1105 if(supcol.le.0) then
1106 xlf = xlf0
1107 !--------------------------------------------------------------
1108 ! nseml: Enhanced melting of snow by accretion of water [LH A29]
1109 ! (T>=T0: ->NR)
1110 !--------------------------------------------------------------
1111 if (qrs(i,k,2).gt.qcrmin) then
1112 sfac = rslope(i,k,2)*n0s*n0sfac(i,k)/qrs(i,k,2)
1113 nseml(i,k) = -sfac*min(max(cliq*supcol*(psacw(i,k))/xlf &
1114 ,-qrs(i,k,2)/dtcld),0.)
1115 endif
1116 endif
1117 if(supcol.gt.0) then
1118 !-------------------------------------------------------------
1119 ! pidep: Deposition/Sublimation rate of ice [HDC 9]
1120 ! (T<T0: QV->QI or QI->QV)
1121 !-------------------------------------------------------------
1122 if(qci(i,k,2).gt.0 .and. ifsat.ne.1) then
1123 xmi = den(i,k)*qci(i,k,2)/xni(i,k)
1124 diameter = dicon * sqrt(xmi)
1125 pidep(i,k) = 4.*diameter*xni(i,k)*(rh(i,k,2)-1.)/work1(i,k,2)
1126 supice = satdt-prevp(i,k)
1127 if(pidep(i,k).lt.0.) then
1128 ! pidep(i,k) = max(max(pidep(i,k),satdt/2),supice)
1129 ! pidep(i,k) = max(pidep(i,k),-qci(i,k,2)/dtcld)
1130 pidep(i,k) = max(max(pidep(i,k),satdt*.5),supice)
1131 pidep(i,k) = max(pidep(i,k),-qci(i,k,2)*rdtcld)
1132 else
1133 ! pidep(i,k) = min(min(pidep(i,k),satdt/2),supice)
1134 pidep(i,k) = min(min(pidep(i,k),satdt*.5),supice)
1135 endif
1136 if(abs(prevp(i,k)+pidep(i,k)).ge.abs(satdt)) ifsat = 1
1137 endif
1138 !-------------------------------------------------------------
1139 ! psdep: deposition/sublimation rate of snow [HDC 14]
1140 ! (QV->QS or QS->QV)
1141 !-------------------------------------------------------------
1142 if(qrs(i,k,2).gt.0. .and. ifsat.ne.1) then
1143 coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
1144 psdep(i,k) = (rh(i,k,2)-1.)*n0sfac(i,k)*(precs1*rslope2(i,k,2) &
1145 + precs2*work2(i,k)*coeres)/work1(i,k,2)
1146 supice = satdt-prevp(i,k)-pidep(i,k)
1147 if(psdep(i,k).lt.0.) then
1148 ! psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)/dtcld)
1149 ! psdep(i,k) = max(max(psdep(i,k),satdt/2),supice)
1150 psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)*rdtcld)
1151 psdep(i,k) = max(max(psdep(i,k),satdt*.5),supice)
1152 else
1153 ! psdep(i,k) = min(min(psdep(i,k),satdt/2),supice)
1154 psdep(i,k) = min(min(psdep(i,k),satdt*.5),supice)
1155 endif
1156 if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)).ge.abs(satdt)) ifsat = 1
1157 endif
1158 !-------------------------------------------------------------
1159 ! pigen: generation(nucleation) of ice from vapor [HL A50] [HDC 7-8]
1160 ! (T<T0: QV->QI)
1161 !-------------------------------------------------------------
1162 if(supsat.gt.0 .and. ifsat.ne.1) then
1163 supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)
1164 xni0 = 1.e3*exp(0.1*supcol)
1165 roqi0 = 4.92e-11*exp(log(xni0)*(1.33))
1166 pigen(i,k) = max(0.,(roqi0/den(i,k)-max(qci(i,k,2),0.))*rdtcld)
1167 pigen(i,k) = min(min(pigen(i,k),satdt),supice)
1168 endif
1169 !
1170 !-------------------------------------------------------------
1171 ! psaut: conversion(aggregation) of ice to snow [HDC 12]
1172 ! (T<T0: QI->QS)
1173 !-------------------------------------------------------------
1174 if(qci(i,k,2).gt.0.) then
1175 qimax = roqimax/den(i,k)
1176 ! psaut(i,k) = max(0.,(qci(i,k,2)-qimax)/dtcld)
1177 psaut(i,k) = max(0.,(qci(i,k,2)-qimax)*rdtcld)
1178 endif
1179 endif
1180 !-------------------------------------------------------------
1181 ! psevp: Evaporation of melting snow [HL A35] [RH83 A27]
1182 ! (T>T0: QS->QV)
1183 !-------------------------------------------------------------
1184 if(supcol.lt.0.) then
1185 if(qrs(i,k,2).gt.0. .and. rh(i,k,1).lt.1.) &
1186 psevp(i,k) = psdep(i,k)*work1(i,k,2)/work1(i,k,1)
1187 ! psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)/dtcld),0.)
1188 psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)*rdtcld),0.)
1189 endif
1190 enddo
1191 enddo
1192 !
1193 !
1194 !----------------------------------------------------------------
1195 ! check mass conservation of generation terms and feedback to the
1196 ! large scale
1197 !
1198 do k = kts, kte
1199 do i = its, ite
1200 if(t(i,k).le.t0c) then
1201 !
1202 ! Q_cloud water
1203 !
1204 value = max(qmin,qci(i,k,1))
1205 source = (praut(i,k)+pracw(i,k)+psacw(i,k)-prevp_s(i,k))*dtcld
1206 if (source.gt.value) then
1207 factor = value/source
1208 praut(i,k) = praut(i,k)*factor
1209 pracw(i,k) = pracw(i,k)*factor
1210 psacw(i,k) = psacw(i,k)*factor
1211 prevp_s(i,k) = prevp_s(i,k)*factor
1212 endif
1213 !
1214 ! Q_cloud ice
1215 !
1216 value = max(qmin,qci(i,k,2))
1217 source = (psaut(i,k)+psaci(i,k)-pigen(i,k)-pidep(i,k))*dtcld
1218 if (source.gt.value) then
1219 factor = value/source
1220 psaut(i,k) = psaut(i,k)*factor
1221 psaci(i,k) = psaci(i,k)*factor
1222 pigen(i,k) = pigen(i,k)*factor
1223 pidep(i,k) = pidep(i,k)*factor
1224 endif
1225 !
1226 ! Q_rain
1227 !
1228 !
1229 value = max(qmin,qrs(i,k,1))
1230 source = (-praut(i,k)-pracw(i,k)-prevp(i,k)+prevp_s(i,k))*dtcld
1231 if (source.gt.value) then
1232 factor = value/source
1233 praut(i,k) = praut(i,k)*factor
1234 pracw(i,k) = pracw(i,k)*factor
1235 prevp(i,k) = prevp(i,k)*factor
1236 prevp_s(i,k) = prevp_s(i,k)*factor
1237 endif
1238 !
1239 ! Q_snow
1240 !
1241 value = max(qmin,qrs(i,k,2))
1242 source = (-psdep(i,k)-psaut(i,k)-psaci(i,k)-psacw(i,k))*dtcld
1243 if (source.gt.value) then
1244 factor = value/source
1245 psdep(i,k) = psdep(i,k)*factor
1246 psaut(i,k) = psaut(i,k)*factor
1247 psaci(i,k) = psaci(i,k)*factor
1248 psacw(i,k) = psacw(i,k)*factor
1249 endif
1250 !
1251 ! N_cloud
1252 !
1253 value = max(ncmin,ncr(i,k,2))
1254 source = (+nraut(i,k)+nccol(i,k)+nracw(i,k)+nsacw(i,k) &
1255 -nrevp(i,k))*dtcld
1256 if (source.gt.value) then
1257 factor = value/source
1258 nraut(i,k) = nraut(i,k)*factor
1259 nccol(i,k) = nccol(i,k)*factor
1260 nracw(i,k) = nracw(i,k)*factor
1261 nsacw(i,k) = nsacw(i,k)*factor
1262 nrevp(i,k) = nrevp(i,k)*factor
1263 endif
1264 !
1265 ! N_rain
1266 !
1267 value = max(nrmin,ncr(i,k,3))
1268 source = (-nraut(i,k)+nrcol(i,k)+nrevp(i,k))*dtcld
1269 if (source.gt.value) then
1270 factor = value/source
1271 nraut(i,k) = nraut(i,k)*factor
1272 nrcol(i,k) = nrcol(i,k)*factor
1273 nrevp(i,k) = nrevp(i,k)*factor
1274 endif
1275 !
1276 work2(i,k)=-(prevp(i,k)+psdep(i,k)+pigen(i,k)+pidep(i,k))
1277 ! update
1278 q(i,k) = q(i,k)+work2(i,k)*dtcld
1279 qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k)+psacw(i,k) &
1280 +prevp_s(i,k))*dtcld,0.)
1281 qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k)+prevp(i,k) &
1282 -prevp_s(i,k))*dtcld,0.)
1283 qci(i,k,2) = max(qci(i,k,2)-(psaut(i,k)+psaci(i,k)-pigen(i,k) &
1284 -pidep(i,k))*dtcld,0.)
1285 qrs(i,k,2) = max(qrs(i,k,2)+(psdep(i,k)+psaut(i,k)+psaci(i,k) &
1286 +psacw(i,k))*dtcld,0.)
1287 ncr(i,k,2) = max(ncr(i,k,2)+(-nraut(i,k)-nccol(i,k)-nracw(i,k) &
1288 -nsacw(i,k)+nrevp(i,k))*dtcld,0.)
1289 ncr(i,k,3) = max(ncr(i,k,3)+(nraut(i,k)-nrcol(i,k)-nrevp(i,k)) &
1290 *dtcld,0.)
1291 xlf = xls-xl(i,k)
1292 xlwork2 = -xls*(psdep(i,k)+pidep(i,k)+pigen(i,k)) &
1293 -xl(i,k)*prevp(i,k)-xlf*psacw(i,k)
1294 t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
1295 else
1296 !
1297 ! Q_cloud water
1298 !
1299 value = max(qmin,qci(i,k,1))
1300 source=(praut(i,k)+pracw(i,k)+psacw(i,k)-prevp_s(i,k))*dtcld
1301 if (source.gt.value) then
1302 factor = value/source
1303 praut(i,k) = praut(i,k)*factor
1304 pracw(i,k) = pracw(i,k)*factor
1305 psacw(i,k) = psacw(i,k)*factor
1306 prevp_s(i,k) = prevp_s(i,k)*factor
1307 endif
1308 !
1309 ! Q_rain
1310 !
1311 value = max(qmin,qrs(i,k,1))
1312 source = (-praut(i,k)-pracw(i,k)-prevp(i,k)+prevp_s(i,k) &
1313 -psacw(i,k))*dtcld
1314 if (source.gt.value) then
1315 factor = value/source
1316 praut(i,k) = praut(i,k)*factor
1317 pracw(i,k) = pracw(i,k)*factor
1318 prevp(i,k) = prevp(i,k)*factor
1319 psacw(i,k) = psacw(i,k)*factor
1320 prevp_s(i,k) = prevp_s(i,k)*factor
1321 endif
1322 !
1323 ! Q_snow
1324 !
1325 value = max(qcrmin,qrs(i,k,2))
1326 source=(-psevp(i,k))*dtcld
1327 if (source.gt.value) then
1328 factor = value/source
1329 psevp(i,k) = psevp(i,k)*factor
1330 endif
1331 !
1332 ! N_cloud
1333 !
1334 value = max(ncmin,ncr(i,k,2))
1335 source = (+nraut(i,k)+nccol(i,k)+nracw(i,k)+nsacw(i,k) &
1336 -nrevp(i,k))*dtcld
1337 if (source.gt.value) then
1338 factor = value/source
1339 nraut(i,k) = nraut(i,k)*factor
1340 nccol(i,k) = nccol(i,k)*factor
1341 nracw(i,k) = nracw(i,k)*factor
1342 nsacw(i,k) = nsacw(i,k)*factor
1343 nrevp(i,k) = nrevp(i,k)*factor
1344 endif
1345 !
1346 ! N_rain
1347 !
1348 value = max(nrmin,ncr(i,k,3))
1349 source = (-nraut(i,k)-nseml(i,k)+nrcol(i,k)+nrevp(i,k))*dtcld
1350 if (source.gt.value) then
1351 factor = value/source
1352 nraut(i,k) = nraut(i,k)*factor
1353 nseml(i,k) = nseml(i,k)*factor
1354 nrcol(i,k) = nrcol(i,k)*factor
1355 nrevp(i,k) = nrevp(i,k)*factor
1356 endif
1357 work2(i,k)=-(prevp(i,k)+psevp(i,k))
1358 ! update
1359 q(i,k) = q(i,k)+work2(i,k)*dtcld
1360 qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k)+psacw(i,k) &
1361 +prevp_s(i,k))*dtcld,0.)
1362 qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k)+prevp(i,k) &
1363 +psacw(i,k)-prevp_s(i,k))*dtcld,0.)
1364 qrs(i,k,2) = max(qrs(i,k,2)+psevp(i,k)*dtcld,0.)
1365 ncr(i,k,2) = max(ncr(i,k,2)+(-nraut(i,k)-nccol(i,k)-nracw(i,k) &
1366 -nsacw(i,k)+nrevp(i,k))*dtcld,0.)
1367 ncr(i,k,3) = max(ncr(i,k,3)+(nraut(i,k)+nseml(i,k)-nrcol(i,k) &
1368 -nrevp(i,k))*dtcld,0.)
1369 xlf = xls-xl(i,k)
1370 xlwork2 = -xl(i,k)*(prevp(i,k)+psevp(i,k))
1371 t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
1372 endif
1373 enddo
1374 enddo
1375 !
1376 ! Inline expansion for fpvs
1377 ! qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1378 ! qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1379 hsub = xls
1380 hvap = xlv0
1381 cvap = cpv
1382 ttp=t0c+0.01
1383 dldt=cvap-cliq
1384 xa=-dldt/rv
1385 xb=xa+hvap/(rv*ttp)
1386 dldti=cvap-cice
1387 xai=-dldti/rv
1388 xbi=xai+hsub/(rv*ttp)
1389 do k = kts, kte
1390 do i = its, ite
1391 tr=ttp/t(i,k)
1392 logtr = log(tr)
1393 qs(i,k,1)=psat*exp(logtr*(xa)+xb*(1.-tr))
1394 qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
1395 qs(i,k,1) = max(qs(i,k,1),qmin)
1396 rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
1397 enddo
1398 enddo
1399 !
1400 do k = kts, kte
1401 do i = its, ite
1402 !-------------------------------------------------------------------
1403 ! rate of change of cloud drop concentration due to CCN activation
1404 ! pcact: QV -> QC [LH 8] [KK 14]
1405 ! ncact: NCCN -> NC [LH A2] [KK 12]
1406 !-------------------------------------------------------------------
1407 if(rh(i,k,1).gt.1.) then
1408 ncact(i,k) = max(0.,((ncr(i,k,1)+ncr(i,k,2)) &
1409 *min(1.,(rh(i,k,1)/satmax)**actk) - ncr(i,k,2)))/dtcld
1410 ncact(i,k) =min(ncact(i,k),max(ncr(i,k,1),0.)/dtcld)
1411 pcact(i,k) = min(4.*pi*denr*(actr*1.E-6)**3*ncact(i,k)/ &
1412 (3.*den(i,k)),max(q(i,k),0.)/dtcld)
1413 q(i,k) = max(q(i,k)-pcact(i,k)*dtcld,0.)
1414 qci(i,k,1) = max(qci(i,k,1)+pcact(i,k)*dtcld,0.)
1415 ncr(i,k,1) = max(ncr(i,k,1)-ncact(i,k)*dtcld,0.)
1416 ncr(i,k,2) = max(ncr(i,k,2)+ncact(i,k)*dtcld,0.)
1417 t(i,k) = t(i,k)+pcact(i,k)*xl(i,k)/cpm(i,k)*dtcld
1418 endif
1419 !---------------------------------------------------------------------
1420 ! pcond:condensational/evaporational rate of cloud water [HL A46] [RH83 A6]
1421 ! if there exists additional water vapor condensated/if
1422 ! evaporation of cloud water is not enough to remove subsaturation
1423 ! (QV->QC or QC->QV)
1424 !---------------------------------------------------------------------
1425 tr=ttp/t(i,k)
1426 qs(i,k,1)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
1427 qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
1428 qs(i,k,1) = max(qs(i,k,1),qmin)
1429 work1(i,k,1) = ((max(q(i,k),qmin)-(qs(i,k,1)))/(1.+(xl(i,k)) &
1430 *(xl(i,k))/(rv*(cpm(i,k)))*(qs(i,k,1))/((t(i,k)) &
1431 *(t(i,k)))))
1432 work2(i,k) = qci(i,k,1)+work1(i,k,1)
1433 pcond(i,k) = min(max(work1(i,k,1)/dtcld,0.),max(q(i,k),0.)/dtcld)
1434 if(qci(i,k,1).gt.0. .and. work1(i,k,1).lt.0.) &
1435 pcond(i,k) = max(work1(i,k,1),-qci(i,k,1))/dtcld
1436 !----------------------------------------------------------------------
1437 ! ncevp: evpration of Cloud number concentration [LH A3]
1438 ! (NC->NCCN)
1439 !----------------------------------------------------------------------
1440 if(pcond(i,k).eq.-qci(i,k,1)/dtcld) then
1441 ncr(i,k,2) = 0.
1442 ncr(i,k,1) = ncr(i,k,1)+ncr(i,k,2)
1443 endif
1444 !
1445 q(i,k) = q(i,k)-pcond(i,k)*dtcld
1446 qci(i,k,1) = max(qci(i,k,1)+pcond(i,k)*dtcld,0.)
1447 t(i,k) = t(i,k)+pcond(i,k)*xl(i,k)/cpm(i,k)*dtcld
1448 enddo
1449 enddo
1450 !
1451 !----------------------------------------------------------------
1452 ! padding for small values
1453 !
1454 do k = kts, kte
1455 do i = its, ite
1456 if(qci(i,k,1).le.qmin) qci(i,k,1) = 0.0
1457 if(qci(i,k,2).le.qmin) qci(i,k,2) = 0.0
1458 enddo
1459 enddo
1460 enddo ! big loops
1461 END SUBROUTINE wdm52d
1462 ! ...................................................................
1463 REAL FUNCTION rgmma(x)
1464 !-------------------------------------------------------------------
1465 IMPLICIT NONE
1466 !-------------------------------------------------------------------
1467 ! rgmma function: use infinite product form
1468 REAL :: euler
1469 PARAMETER (euler=0.577215664901532)
1470 REAL :: x, y
1471 INTEGER :: i
1472 if(x.eq.1.)then
1473 rgmma=0.
1474 else
1475 rgmma=x*exp(euler*x)
1476 do i=1,10000
1477 y=float(i)
1478 rgmma=rgmma*(1.000+x/y)*exp(-x/y)
1479 enddo
1480 rgmma=1./rgmma
1481 endif
1482 END FUNCTION rgmma
1483 !
1484 !--------------------------------------------------------------------------
1485 REAL FUNCTION fpvs(t,ice,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c)
1486 !--------------------------------------------------------------------------
1487 IMPLICIT NONE
1488 !--------------------------------------------------------------------------
1489 REAL t,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c,dldt,xa,xb,dldti, &
1490 xai,xbi,ttp,tr
1491 INTEGER ice
1492 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1493 ttp=t0c+0.01
1494 dldt=cvap-cliq
1495 xa=-dldt/rv
1496 xb=xa+hvap/(rv*ttp)
1497 dldti=cvap-cice
1498 xai=-dldti/rv
1499 xbi=xai+hsub/(rv*ttp)
1500 tr=ttp/t
1501 if(t.lt.ttp .and. ice.eq.1) then
1502 fpvs=psat*exp(log(tr)*(xai))*exp(xbi*(1.-tr))
1503 else
1504 fpvs=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
1505 endif
1506 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1507 END FUNCTION fpvs
1508 !-------------------------------------------------------------------
1509 SUBROUTINE wdm5init(den0,denr,dens,cl,cpv,ccn0,allowed_to_read)
1510 !-------------------------------------------------------------------
1511 IMPLICIT NONE
1512 !-------------------------------------------------------------------
1513 !.... constants which may not be tunable
1514 REAL, INTENT(IN) :: den0,denr,dens,cl,cpv,ccn0
1515 LOGICAL, INTENT(IN) :: allowed_to_read
1516 !
1517 pi = 4.*atan(1.)
1518 xlv1 = cl-cpv
1519 !
1520 qc0 = 4./3.*pi*denr*r0**3*xncr/den0 ! 0.419e-3 -- .61e-3
1521 qck1 = .104*9.8*peaut/(xncr*denr)**(1./3.)/xmyu*den0**(4./3.) ! 7.03
1522 pidnc = pi*denr/6.
1523 !
1524 bvtr1 = 1.+bvtr
1525 bvtr2 = 2.+bvtr
1526 bvtr3 = 3.+bvtr
1527 bvtr4 = 4.+bvtr
1528 bvtr5 = 5.+bvtr
1529 bvtr7 = 7.+bvtr
1530 bvtr2o5 = 2.5+.5*bvtr
1531 bvtr3o5 = 3.5+.5*bvtr
1532 g1pbr = rgmma(bvtr1)
1533 g2pbr = rgmma(bvtr2)
1534 g3pbr = rgmma(bvtr3)
1535 g4pbr = rgmma(bvtr4) ! 17.837825
1536 g5pbr = rgmma(bvtr5)
1537 g7pbr = rgmma(bvtr7)
1538 g5pbro2 = rgmma(bvtr2o5)
1539 g7pbro2 = rgmma(bvtr3o5)
1540 pvtr = avtr*g5pbr/24.
1541 pvtrn = avtr*g2pbr
1542 eacrr = 1.0
1543 pacrr = pi*n0r*avtr*g3pbr*.25*eacrr
1544 precr1 = 2.*pi*1.56
1545 precr2 = 2.*pi*.31*avtr**.5*g7pbro2
1546 pidn0r = pi*denr*n0r
1547 pidnr = 4.*pi*denr
1548 xmmax = (dimax/dicon)**2
1549 roqimax = 2.08e22*dimax**8
1550 !
1551 bvts1 = 1.+bvts
1552 bvts2 = 2.5+.5*bvts
1553 bvts3 = 3.+bvts
1554 bvts4 = 4.+bvts
1555 g1pbs = rgmma(bvts1) !.8875
1556 g3pbs = rgmma(bvts3)
1557 g4pbs = rgmma(bvts4) ! 12.0786
1558 g5pbso2 = rgmma(bvts2)
1559 pvts = avts*g4pbs/6.
1560 pacrs = pi*n0s*avts*g3pbs*.25
1561 precs1 = 4.*n0s*.65
1562 precs2 = 4.*n0s*.44*avts**.5*g5pbso2
1563 pidn0s = pi*dens*n0s
1564 pacrc = pi*n0s*avts*g3pbs*.25*eacrc
1565 !
1566 rslopecmax = 1./lamdacmax
1567 rslopermax = 1./lamdarmax
1568 rslopesmax = 1./lamdasmax
1569 rsloperbmax = rslopermax ** bvtr
1570 rslopesbmax = rslopesmax ** bvts
1571 rslopec2max = rslopecmax * rslopecmax
1572 rsloper2max = rslopermax * rslopermax
1573 rslopes2max = rslopesmax * rslopesmax
1574 rslopec3max = rslopec2max * rslopecmax
1575 rsloper3max = rsloper2max * rslopermax
1576 rslopes3max = rslopes2max * rslopesmax
1577 !
1578 END SUBROUTINE wdm5init
1579 !------------------------------------------------------------------------------
1580 subroutine slope_wdm5(qrs,ncr,den,denfac,t,rslope,rslopeb,rslope2,rslope3, &
1581 vt,vtn,its,ite,kts,kte)
1582 IMPLICIT NONE
1583 INTEGER :: its,ite, jts,jte, kts,kte
1584 REAL, DIMENSION( its:ite , kts:kte,2) :: &
1585 qrs, &
1586 rslope, &
1587 rslopeb, &
1588 rslope2, &
1589 rslope3, &
1590 vt
1591 REAL, DIMENSION( its:ite , kts:kte) :: &
1592 ncr, &
1593 vtn, &
1594 den, &
1595 denfac, &
1596 t
1597 REAL, PARAMETER :: t0c = 273.15
1598 REAL, DIMENSION( its:ite , kts:kte ) :: &
1599 n0sfac
1600 REAL :: lamdar, lamdas, x, y, z, supcol
1601 integer :: i, j, k
1602 !----------------------------------------------------------------
1603 ! size distributions: (x=mixing ratio, y=air density):
1604 ! valid for mixing ratio > 1.e-9 kg/kg.
1605 !
1606 ! Optimizatin : A**B => exp(log(A)*(B))
1607 lamdar(x,y,z)= exp(log(((pidnr*z)/(x*y)))*((.33333333)))
1608 lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y))) ! (pidn0s*z/(x*y))**.25
1609 !
1610 do k = kts, kte
1611 do i = its, ite
1612 supcol = t0c-t(i,k)
1613 !---------------------------------------------------------------
1614 ! n0s: Intercept parameter for snow [m-4] [HDC 6]
1615 !---------------------------------------------------------------
1616 n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
1617 if(qrs(i,k,1).le.qcrmin .or. ncr(i,k).le.nrmin ) then
1618 rslope(i,k,1) = rslopermax
1619 rslopeb(i,k,1) = rsloperbmax
1620 rslope2(i,k,1) = rsloper2max
1621 rslope3(i,k,1) = rsloper3max
1622 else
1623 rslope(i,k,1) = min(1./lamdar(qrs(i,k,1),den(i,k),ncr(i,k)),1.e-3)
1624 rslopeb(i,k,1) = rslope(i,k,1)**bvtr
1625 rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
1626 rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
1627 endif
1628 if(qrs(i,k,2).le.qcrmin) then
1629 rslope(i,k,2) = rslopesmax
1630 rslopeb(i,k,2) = rslopesbmax
1631 rslope2(i,k,2) = rslopes2max
1632 rslope3(i,k,2) = rslopes3max
1633 else
1634 rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k))
1635 rslopeb(i,k,2) = rslope(i,k,2)**bvts
1636 rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
1637 rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
1638 endif
1639 vt(i,k,1) = pvtr*rslopeb(i,k,1)*denfac(i,k)
1640 vt(i,k,2) = pvts*rslopeb(i,k,2)*denfac(i,k)
1641 vtn(i,k) = pvtrn*rslopeb(i,k,1)*denfac(i,k)
1642 if(qrs(i,k,1).le.0.0) vt(i,k,1) = 0.0
1643 if(qrs(i,k,2).le.0.0) vt(i,k,2) = 0.0
1644 if(ncr(i,k).le.0.0) vtn(i,k) = 0.0
1645 enddo
1646 enddo
1647 END subroutine slope_wdm5
1648 !-----------------------------------------------------------------------------
1649 subroutine slope_rain(qrs,ncr,den,denfac,t,rslope,rslopeb,rslope2,rslope3, &
1650 vt,vtn,its,ite,kts,kte)
1651 IMPLICIT NONE
1652 INTEGER :: its,ite, jts,jte, kts,kte
1653 REAL, DIMENSION( its:ite , kts:kte) :: &
1654 qrs, &
1655 ncr, &
1656 rslope, &
1657 rslopeb, &
1658 rslope2, &
1659 rslope3, &
1660 vt, &
1661 vtn, &
1662 den, &
1663 denfac, &
1664 t
1665 REAL, PARAMETER :: t0c = 273.15
1666 REAL, DIMENSION( its:ite , kts:kte ) :: &
1667 n0sfac
1668 REAL :: lamdar, x, y, z, supcol
1669 integer :: i, j, k
1670 !----------------------------------------------------------------
1671 ! size distributions: (x=mixing ratio, y=air density):
1672 ! valid for mixing ratio > 1.e-9 kg/kg.
1673 lamdar(x,y,z)= exp(log(((pidnr*z)/(x*y)))*((.33333333)))
1674 !
1675 do k = kts, kte
1676 do i = its, ite
1677 if(qrs(i,k).le.qcrmin .or. ncr(i,k).le.nrmin) then
1678 rslope(i,k) = rslopermax
1679 rslopeb(i,k) = rsloperbmax
1680 rslope2(i,k) = rsloper2max
1681 rslope3(i,k) = rsloper3max
1682 else
1683 rslope(i,k) = min(1./lamdar(qrs(i,k),den(i,k),ncr(i,k)),1.e-3)
1684 rslopeb(i,k) = rslope(i,k)**bvtr
1685 rslope2(i,k) = rslope(i,k)*rslope(i,k)
1686 rslope3(i,k) = rslope2(i,k)*rslope(i,k)
1687 endif
1688 vt(i,k) = pvtr*rslopeb(i,k)*denfac(i,k)
1689 vtn(i,k) = pvtrn*rslopeb(i,k)*denfac(i,k)
1690 if(qrs(i,k).le.0.0) vt(i,k) = 0.0
1691 if(ncr(i,k).le.0.0) vtn(i,k) = 0.0
1692 enddo
1693 enddo
1694 END subroutine slope_rain
1695 !------------------------------------------------------------------------------
1696 subroutine slope_snow(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3, &
1697 vt,its,ite,kts,kte)
1698 IMPLICIT NONE
1699 INTEGER :: its,ite, jts,jte, kts,kte
1700 REAL, DIMENSION( its:ite , kts:kte) :: &
1701 qrs, &
1702 rslope, &
1703 rslopeb, &
1704 rslope2, &
1705 rslope3, &
1706 vt, &
1707 den, &
1708 denfac, &
1709 t
1710 REAL, PARAMETER :: t0c = 273.15
1711 REAL, DIMENSION( its:ite , kts:kte ) :: &
1712 n0sfac
1713 REAL :: lamdas, x, y, z, supcol
1714 integer :: i, j, k
1715 !----------------------------------------------------------------
1716 ! size distributions: (x=mixing ratio, y=air density):
1717 ! valid for mixing ratio > 1.e-9 kg/kg.
1718 lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y))) ! (pidn0s*z/(x*y))**.25
1719 !
1720 do k = kts, kte
1721 do i = its, ite
1722 supcol = t0c-t(i,k)
1723 !---------------------------------------------------------------
1724 ! n0s: Intercept parameter for snow [m-4] [HDC 6]
1725 !---------------------------------------------------------------
1726 n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
1727 if(qrs(i,k).le.qcrmin)then
1728 rslope(i,k) = rslopesmax
1729 rslopeb(i,k) = rslopesbmax
1730 rslope2(i,k) = rslopes2max
1731 rslope3(i,k) = rslopes3max
1732 else
1733 rslope(i,k) = 1./lamdas(qrs(i,k),den(i,k),n0sfac(i,k))
1734 rslopeb(i,k) = rslope(i,k)**bvts
1735 rslope2(i,k) = rslope(i,k)*rslope(i,k)
1736 rslope3(i,k) = rslope2(i,k)*rslope(i,k)
1737 endif
1738 vt(i,k) = pvts*rslopeb(i,k)*denfac(i,k)
1739 if(qrs(i,k).le.0.0) vt(i,k) = 0.0
1740 enddo
1741 enddo
1742 END subroutine slope_snow
1743 !-------------------------------------------------------------------
1744 SUBROUTINE nislfv_rain_plm(im,km,denl,denfacl,tkl,dzl,wwl,rql,precip,dt,id,iter)
1745 !-------------------------------------------------------------------
1746 !
1747 ! for non-iteration semi-Lagrangain forward advection for cloud
1748 ! with mass conservation and positive definite advection
1749 ! 2nd order interpolation with monotonic piecewise linear method
1750 ! this routine is under assumption of decfl < 1 for semi_Lagrangian
1751 !
1752 ! dzl depth of model layer in meter
1753 ! wwl terminal velocity at model layer m/s
1754 ! rql cloud density*mixing ration
1755 ! precip precipitation
1756 ! dt time step
1757 ! id kind of precip: 0 test case; 1 raindrop 2: snow
1758 ! iter how many time to guess mean terminal velocity: 0 pure forward.
1759 ! 0 : use departure wind for advection
1760 ! 1 : use mean wind for advection
1761 ! > 1 : use mean wind after iter-1 iterations
1762 !
1763 ! author: hann-ming henry juang <henry.juang@noaa.gov>
1764 ! implemented by song-you hong
1765 !
1766 implicit none
1767 integer im,km,id
1768 real dt
1769 real dzl(im,km),wwl(im,km),rql(im,km),precip(im)
1770 real denl(im,km),denfacl(im,km),tkl(im,km)
1771 !
1772 integer i,k,n,m,kk,kb,kt,iter
1773 real tl,tl2,qql,dql,qqd
1774 real th,th2,qqh,dqh
1775 real zsum,qsum,dim,dip,c1,con1,fa1,fa2
1776 real allold, allnew, zz, dzamin, cflmax, decfl
1777 real dz(km), ww(km), qq(km), wd(km), wa(km), was(km)
1778 real den(km), denfac(km), tk(km)
1779 real wi(km+1), zi(km+1), za(km+1)
1780 real qn(km), qr(km),tmp(km),tmp1(km),tmp2(km),tmp3(km)
1781 real dza(km+1), qa(km+1), qmi(km+1), qpi(km+1)
1782 !
1783 precip(:) = 0.0
1784 !
1785 i_loop : do i=1,im
1786 ! -----------------------------------
1787 dz(:) = dzl(i,:)
1788 qq(:) = rql(i,:)
1789 ww(:) = wwl(i,:)
1790 den(:) = denl(i,:)
1791 denfac(:) = denfacl(i,:)
1792 tk(:) = tkl(i,:)
1793 ! skip for no precipitation for all layers
1794 allold = 0.0
1795 do k=1,km
1796 allold = allold + qq(k)
1797 enddo
1798 if(allold.le.0.0) then
1799 cycle i_loop
1800 endif
1801 !
1802 ! compute interface values
1803 zi(1)=0.0
1804 do k=1,km
1805 zi(k+1) = zi(k)+dz(k)
1806 enddo
1807 !
1808 ! save departure wind
1809 wd(:) = ww(:)
1810 n=1
1811 100 continue
1812 ! plm is 2nd order, we can use 2nd order wi or 3rd order wi
1813 ! 2nd order interpolation to get wi
1814 wi(1) = ww(1)
1815 wi(km+1) = ww(km)
1816 do k=2,km
1817 wi(k) = (ww(k)*dz(k-1)+ww(k-1)*dz(k))/(dz(k-1)+dz(k))
1818 enddo
1819 ! 3rd order interpolation to get wi
1820 fa1 = 9./16.
1821 fa2 = 1./16.
1822 wi(1) = ww(1)
1823 wi(2) = 0.5*(ww(2)+ww(1))
1824 do k=3,km-1
1825 wi(k) = fa1*(ww(k)+ww(k-1))-fa2*(ww(k+1)+ww(k-2))
1826 enddo
1827 wi(km) = 0.5*(ww(km)+ww(km-1))
1828 wi(km+1) = ww(km)
1829 !
1830 ! terminate of top of raingroup
1831 do k=2,km
1832 if( ww(k).eq.0.0 ) wi(k)=ww(k-1)
1833 enddo
1834 !
1835 ! diffusivity of wi
1836 con1 = 0.05
1837 do k=km,1,-1
1838 decfl = (wi(k+1)-wi(k))*dt/dz(k)
1839 if( decfl .gt. con1 ) then
1840 wi(k) = wi(k+1) - con1*dz(k)/dt
1841 endif
1842 enddo
1843 ! compute arrival point
1844 do k=1,km+1
1845 za(k) = zi(k) - wi(k)*dt
1846 enddo
1847 !
1848 do k=1,km
1849 dza(k) = za(k+1)-za(k)
1850 enddo
1851 dza(km+1) = zi(km+1) - za(km+1)
1852 !
1853 ! computer deformation at arrival point
1854 do k=1,km
1855 qa(k) = qq(k)*dz(k)/dza(k)
1856 qr(k) = qa(k)/den(k)
1857 enddo
1858 qa(km+1) = 0.0
1859 ! call maxmin(km,1,qa,' arrival points ')
1860 !
1861 ! compute arrival terminal velocity, and estimate mean terminal velocity
1862 ! then back to use mean terminal velocity
1863 if( n.le.iter ) then
1864 ! if (id.eq.1) then
1865 !
1866 ! call slope_rain(qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,1,1,1,km)
1867 ! else
1868 call slope_snow(qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,1,1,1,km)
1869 ! endif
1870 if( n.ge.2 ) wa(1:km)=0.5*(wa(1:km)+was(1:km))
1871 do k=1,km
1872 !#ifdef DEBUG
1873 ! print*,' slope_wsm3 ',qr(k)*1000.,den(k),denfac(k),tk(k),tmp(k),tmp1(k),tmp2(k),ww(k),wa(k)
1874 !#endif
1875 ! mean wind is average of departure and new arrival winds
1876 ww(k) = 0.5* ( wd(k)+wa(k) )
1877 enddo
1878 was(:) = wa(:)
1879 n=n+1
1880 go to 100
1881 endif
1882 !
1883 ! estimate values at arrival cell interface with monotone
1884 do k=2,km
1885 dip=(qa(k+1)-qa(k))/(dza(k+1)+dza(k))
1886 dim=(qa(k)-qa(k-1))/(dza(k-1)+dza(k))
1887 if( dip*dim.le.0.0 ) then
1888 qmi(k)=qa(k)
1889 qpi(k)=qa(k)
1890 else
1891 qpi(k)=qa(k)+0.5*(dip+dim)*dza(k)
1892 qmi(k)=2.0*qa(k)-qpi(k)
1893 if( qpi(k).lt.0.0 .or. qmi(k).lt.0.0 ) then
1894 qpi(k) = qa(k)
1895 qmi(k) = qa(k)
1896 endif
1897 endif
1898 enddo
1899 qpi(1)=qa(1)
1900 qmi(1)=qa(1)
1901 qmi(km+1)=qa(km+1)
1902 qpi(km+1)=qa(km+1)
1903 !
1904 ! interpolation to regular point
1905 qn = 0.0
1906 kb=1
1907 kt=1
1908 intp : do k=1,km
1909 kb=max(kb-1,1)
1910 kt=max(kt-1,1)
1911 ! find kb and kt
1912 if( zi(k).ge.za(km+1) ) then
1913 exit intp
1914 else
1915 find_kb : do kk=kb,km
1916 if( zi(k).le.za(kk+1) ) then
1917 kb = kk
1918 exit find_kb
1919 else
1920 cycle find_kb
1921 endif
1922 enddo find_kb
1923 find_kt : do kk=kt,km
1924 if( zi(k+1).le.za(kk) ) then
1925 kt = kk
1926 exit find_kt
1927 else
1928 cycle find_kt
1929 endif
1930 enddo find_kt
1931 kt = kt - 1
1932 ! compute q with piecewise constant method
1933 if( kt.eq.kb ) then
1934 tl=(zi(k)-za(kb))/dza(kb)
1935 th=(zi(k+1)-za(kb))/dza(kb)
1936 tl2=tl*tl
1937 th2=th*th
1938 qqd=0.5*(qpi(kb)-qmi(kb))
1939 qqh=qqd*th2+qmi(kb)*th
1940 qql=qqd*tl2+qmi(kb)*tl
1941 qn(k) = (qqh-qql)/(th-tl)
1942 else if( kt.gt.kb ) then
1943 tl=(zi(k)-za(kb))/dza(kb)
1944 tl2=tl*tl
1945 qqd=0.5*(qpi(kb)-qmi(kb))
1946 qql=qqd*tl2+qmi(kb)*tl
1947 dql = qa(kb)-qql
1948 zsum = (1.-tl)*dza(kb)
1949 qsum = dql*dza(kb)
1950 if( kt-kb.gt.1 ) then
1951 do m=kb+1,kt-1
1952 zsum = zsum + dza(m)
1953 qsum = qsum + qa(m) * dza(m)
1954 enddo
1955 endif
1956 th=(zi(k+1)-za(kt))/dza(kt)
1957 th2=th*th
1958 qqd=0.5*(qpi(kt)-qmi(kt))
1959 dqh=qqd*th2+qmi(kt)*th
1960 zsum = zsum + th*dza(kt)
1961 qsum = qsum + dqh*dza(kt)
1962 qn(k) = qsum/zsum
1963 endif
1964 cycle intp
1965 endif
1966 !
1967 enddo intp
1968 !
1969 ! rain out
1970 sum_precip: do k=1,km
1971 if( za(k).lt.0.0 .and. za(k+1).lt.0.0 ) then
1972 precip(i) = precip(i) + qa(k)*dza(k)
1973 cycle sum_precip
1974 else if ( za(k).lt.0.0 .and. za(k+1).ge.0.0 ) then
1975 precip(i) = precip(i) + qa(k)*(0.0-za(k))
1976 exit sum_precip
1977 endif
1978 exit sum_precip
1979 enddo sum_precip
1980 !
1981 ! replace the new values
1982 rql(i,:) = qn(:)
1983 !
1984 ! ----------------------------------
1985 enddo i_loop
1986 !
1987 END SUBROUTINE nislfv_rain_plm
1988 END MODULE module_mp_wdm5