| blob | 9ee0e0dac51f7e0a194c1dc24fb22a37718a6e53 |
1 #if ( RWORDSIZE == 4 )
2 # define VREC vsrec
3 # define VSQRT vssqrt
4 #else
5 # define VREC vrec
6 # define VSQRT vsqrt
7 #endif
9 MODULE module_mp_wdm6
10 !
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 :: n0g = 4.e6 ! intercept parameter graupel
16 REAL, PARAMETER, PRIVATE :: avtr = 841.9 ! a constant for terminal velocity of rain
17 REAL, PARAMETER, PRIVATE :: bvtr = 0.8 ! a constant for terminal velocity of rain
18 REAL, PARAMETER, PRIVATE :: r0 = .8e-5 ! 8 microm in contrast to 10 micro m
19 REAL, PARAMETER, PRIVATE :: peaut = .55 ! collection efficiency
20 REAL, PARAMETER, PRIVATE :: xncr = 3.e8 ! maritime cloud in contrast to 3.e8 in tc80
21 REAL, PARAMETER, PRIVATE :: xmyu = 1.718e-5 ! the dynamic viscosity kgm-1s-1
22 REAL, PARAMETER, PRIVATE :: avts = 11.72 ! a constant for terminal velocity of snow
23 REAL, PARAMETER, PRIVATE :: bvts = .41 ! a constant for terminal velocity of snow
24 REAL, PARAMETER, PRIVATE :: avtg = 330. ! a constant for terminal velocity of graupel
25 REAL, PARAMETER, PRIVATE :: bvtg = 0.8 ! a constant for terminal velocity of graupel
26 REAL, PARAMETER, PRIVATE :: deng = 500. ! density of graupel
27 REAL, PARAMETER, PRIVATE :: n0smax = 1.e11 ! maximum n0s (t=-90C unlimited)
28 REAL, PARAMETER, PRIVATE :: lamdacmax = 1.e10 ! limited maximum value for slope parameter of cloud water
29 REAL, PARAMETER, PRIVATE :: lamdarmax = 1.e8 ! limited maximum value for slope parameter of rain
30 REAL, PARAMETER, PRIVATE :: lamdasmax = 1.e5 ! limited maximum value for slope parameter of snow
31 REAL, PARAMETER, PRIVATE :: lamdagmax = 6.e4 ! limited maximum value for slope parameter of graupel
32 REAL, PARAMETER, PRIVATE :: dicon = 11.9 ! constant for the cloud-ice diamter
33 REAL, PARAMETER, PRIVATE :: dimax = 500.e-6 ! limited maximum value for the cloud-ice diamter
34 REAL, PARAMETER, PRIVATE :: n0s = 2.e6 ! temperature dependent intercept parameter snow
35 REAL, PARAMETER, PRIVATE :: alpha = .12 ! .122 exponen factor for n0s
36 REAL, PARAMETER, PRIVATE :: pfrz1 = 100. ! constant in Biggs freezing
37 REAL, PARAMETER, PRIVATE :: pfrz2 = 0.66 ! constant in Biggs freezing
38 REAL, PARAMETER, PRIVATE :: qcrmin = 1.e-9 ! minimun values for qr, qs, and qg
39 REAL, PARAMETER, PRIVATE :: ncmin = 1.e1 ! minimum value for Nc
40 REAL, PARAMETER, PRIVATE :: nrmin = 1.e-2 ! minimum value for Nr
41 REAL, PARAMETER, PRIVATE :: eacrc = 1.0 ! Snow/cloud-water collection efficiency
42 REAL, PARAMETER, PRIVATE :: dens = 100.0 ! Density of snow
43 REAL, PARAMETER, PRIVATE :: qs0 = 6.e-4 ! threshold amount for aggretion to occur
44 !
45 REAL, PARAMETER, PRIVATE :: satmax = 1.0048 ! maximum saturation value for CCN activation
46 ! 1.008 for maritime /1.0048 for conti
47 REAL, PARAMETER, PRIVATE :: actk = 0.6 ! parameter for the CCN activation
48 REAL, PARAMETER, PRIVATE :: actr = 1.5 ! radius of activated CCN drops
49 REAL, PARAMETER, PRIVATE :: ncrk1 = 3.03e3 ! Long's collection kernel coefficient
50 REAL, PARAMETER, PRIVATE :: ncrk2 = 2.59e15 ! Long's collection kernel coefficient
51 REAL, PARAMETER, PRIVATE :: di100 = 1.e-4 ! parameter related with accretion and collection of cloud drops
52 REAL, PARAMETER, PRIVATE :: di600 = 6.e-4 ! parameter related with accretion and collection of cloud drops
53 REAL, PARAMETER, PRIVATE :: di2000 = 2000.e-6 ! parameter related with accretion and collection of cloud drops
54 REAL, PARAMETER, PRIVATE :: di82 = 82.e-6 ! dimater related with raindrops evaporation
55 REAL, PARAMETER, PRIVATE :: di15 = 15.e-6 ! auto conversion takes place beyond this diameter
56 !
57 REAL, SAVE :: &
58 qc0,qck1,pidnc,bvtr1,bvtr2,bvtr3,bvtr4,bvtr5, &
59 bvtr6,bvtr7, bvtr2o5,bvtr3o5, &
60 g1pbr,g2pbr,g3pbr,g4pbr,g5pbr,g6pbr,g7pbr, &
61 g5pbro2,g7pbro2,pi, &
62 pvtr,pvtrn,eacrr,pacrr,pidn0r,pidnr, &
63 precr1,precr2,xmmax,roqimax,bvts1,bvts2, &
64 bvts3,bvts4,g1pbs,g3pbs,g4pbs,g5pbso2, &
65 pvts,pacrs,precs1,precs2,pidn0s,xlv1,pacrc, &
66 bvtg1,bvtg2,bvtg3,bvtg4,g1pbg,g3pbg,g4pbg, &
67 g5pbgo2,pvtg,pacrg,precg1,precg2,pidn0g, &
68 rslopecmax,rslopec2max,rslopec3max, &
69 rslopermax,rslopesmax,rslopegmax, &
70 rsloperbmax,rslopesbmax,rslopegbmax, &
71 rsloper2max,rslopes2max,rslopeg2max, &
72 rsloper3max,rslopes3max,rslopeg3max
73 CONTAINS
74 !===================================================================
75 !
76 SUBROUTINE wdm6(th, q, qc, qr, qi, qs, qg, &
77 nn, nc, nr, &
78 den, pii, p, delz, &
79 delt,g, cpd, cpv, ccn0, rd, rv, t0c, &
80 ep1, ep2, qmin, &
81 XLS, XLV0, XLF0, den0, denr, &
82 cliq,cice,psat, &
83 rain, rainncv, &
84 snow, snowncv, &
85 sr, &
86 graupel, graupelncv, &
87 itimestep, &
88 ids,ide, jds,jde, kds,kde, &
89 ims,ime, jms,jme, kms,kme, &
90 its,ite, jts,jte, kts,kte &
91 )
92 !-------------------------------------------------------------------
93 IMPLICIT NONE
94 !-------------------------------------------------------------------
95 !
96 ! This code is a WRF double-moment 6-class GRAUPEL phase
97 ! microphyiscs scheme (WDM6). The WDM microphysics scheme predicts
98 ! number concentrations for warm rain species including clouds and
99 ! rain. cloud condensation nuclei (CCN) is also predicted.
100 ! The cold rain species including ice, snow, graupel follow the
101 ! WRF single-moment 6-class microphysics (WSM6, Hong and Lim 2006)
102 ! in which theoretical background for WSM ice phase microphysics is
103 ! based on Hong et al. (2004). A new mixed-phase terminal velocity
104 ! for precipitating ice is introduced in WSM6 (Dudhia et al. 2008).
105 ! The WDM scheme is described in Lim and Hong (2009).
106 ! All units are in m.k.s. and source/sink terms in kgkg-1s-1.
107 !
108 ! WDM6 cloud scheme
109 !
110 ! Coded by Kyo-Sun Lim and Song-You Hong (Yonsei Univ.) Fall 2008
111 !
112 ! Implemented by Kyo-Sun Lim and Jimy Dudhia (NCAR) Winter 2008
113 !
114 ! Reference) Lim and Hong (LH, 2010) Mon. Wea. Rev.
115 ! Hong, Dudhia, Chen (HDC, 2004) Mon. Wea. Rev.
116 ! Hong and Lim (HL, 2006) J. Korean Meteor. Soc.
117 ! Cohard and Pinty (CP, 2000) Quart. J. Roy. Meteor. Soc.
118 ! Khairoutdinov and Kogan (KK, 2000) Mon. Wea. Rev.
119 ! Dudhia, Hong and Lim (DHL, 2008) J. Meteor. Soc. Japan
120 !
121 ! Lin, Farley, Orville (LFO, 1983) J. Appl. Meteor.
122 ! Rutledge, Hobbs (RH83, 1983) J. Atmos. Sci.
123 ! Rutledge, Hobbs (RH84, 1984) J. Atmos. Sci.
124 !
125 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
126 ims,ime, jms,jme, kms,kme , &
127 its,ite, jts,jte, kts,kte
128 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
129 INTENT(INOUT) :: &
130 th, &
131 q, &
132 qc, &
133 qi, &
134 qr, &
135 qs, &
136 qg, &
137 nn, &
138 nc, &
139 nr
140 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
141 INTENT(IN ) :: &
142 den, &
143 pii, &
144 p, &
145 delz
146 REAL, INTENT(IN ) :: delt, &
147 g, &
148 rd, &
149 rv, &
150 t0c, &
151 den0, &
152 cpd, &
153 cpv, &
154 ccn0, &
155 ep1, &
156 ep2, &
157 qmin, &
158 XLS, &
159 XLV0, &
160 XLF0, &
161 cliq, &
162 cice, &
163 psat, &
164 denr
165 INTEGER, INTENT(IN ) :: itimestep
166 REAL, DIMENSION( ims:ime , jms:jme ), &
167 INTENT(INOUT) :: rain, &
168 rainncv, &
169 sr
170 REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
171 INTENT(INOUT) :: snow, &
172 snowncv
173 REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
174 INTENT(INOUT) :: graupel, &
175 graupelncv
176 ! LOCAL VAR
177 REAL, DIMENSION( its:ite , kts:kte ) :: t
178 REAL, DIMENSION( its:ite , kts:kte, 2 ) :: qci
179 REAL, DIMENSION( its:ite , kts:kte, 3 ) :: qrs, ncr
180 INTEGER :: i,j,k
181 !-------------------------------------------------------------------
182 IF (itimestep .eq. 1) THEN
183 DO j=jms,jme
184 DO k=kms,kme
185 DO i=ims,ime
186 nn(i,k,j) = ccn0
187 ENDDO
188 ENDDO
189 ENDDO
190 ENDIF
191 !
192 DO j=jts,jte
193 DO k=kts,kte
194 DO i=its,ite
195 t(i,k)=th(i,k,j)*pii(i,k,j)
196 qci(i,k,1) = qc(i,k,j)
197 qci(i,k,2) = qi(i,k,j)
198 qrs(i,k,1) = qr(i,k,j)
199 qrs(i,k,2) = qs(i,k,j)
200 qrs(i,k,3) = qg(i,k,j)
201 ncr(i,k,1) = nn(i,k,j)
202 ncr(i,k,2) = nc(i,k,j)
203 ncr(i,k,3) = nr(i,k,j)
204 ENDDO
205 ENDDO
206 ! Sending array starting locations of optional variables may cause
207 ! troubles, so we explicitly change the call.
208 CALL wdm62D(t, q(ims,kms,j), qci, qrs, ncr &
209 ,den(ims,kms,j) &
210 ,p(ims,kms,j), delz(ims,kms,j) &
211 ,delt,g, cpd, cpv, ccn0, rd, rv, t0c &
212 ,ep1, ep2, qmin &
213 ,XLS, XLV0, XLF0, den0, denr &
214 ,cliq,cice,psat &
215 ,j &
216 ,rain(ims,j),rainncv(ims,j) &
217 ,sr(ims,j) &
218 ,ids,ide, jds,jde, kds,kde &
219 ,ims,ime, jms,jme, kms,kme &
220 ,its,ite, jts,jte, kts,kte &
221 ,snow(ims,j),snowncv(ims,j) &
222 ,graupel(ims,j),graupelncv(ims,j) &
223 )
224 DO K=kts,kte
225 DO I=its,ite
226 th(i,k,j)=t(i,k)/pii(i,k,j)
227 qc(i,k,j) = qci(i,k,1)
228 qi(i,k,j) = qci(i,k,2)
229 qr(i,k,j) = qrs(i,k,1)
230 qs(i,k,j) = qrs(i,k,2)
231 qg(i,k,j) = qrs(i,k,3)
232 nn(i,k,j) = ncr(i,k,1)
233 nc(i,k,j) = ncr(i,k,2)
234 nr(i,k,j) = ncr(i,k,3)
235 ENDDO
236 ENDDO
237 ENDDO
238 END SUBROUTINE wdm6
239 !===================================================================
240 !
241 SUBROUTINE wdm62D(t, q, qci, qrs, ncr, den, p, delz &
242 ,delt,g, cpd, cpv, ccn0, rd, rv, t0c &
243 ,ep1, ep2, qmin &
244 ,XLS, XLV0, XLF0, den0, denr &
245 ,cliq,cice,psat &
246 ,lat &
247 ,rain,rainncv &
248 ,sr &
249 ,ids,ide, jds,jde, kds,kde &
250 ,ims,ime, jms,jme, kms,kme &
251 ,its,ite, jts,jte, kts,kte &
252 ,snow,snowncv &
253 ,graupel,graupelncv &
254 )
255 !-------------------------------------------------------------------
256 IMPLICIT NONE
257 !-------------------------------------------------------------------
258 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
259 ims,ime, jms,jme, kms,kme , &
260 its,ite, jts,jte, kts,kte , &
261 lat
262 REAL, DIMENSION( its:ite , kts:kte ), &
263 INTENT(INOUT) :: &
264 t
265 REAL, DIMENSION( its:ite , kts:kte, 2 ), &
266 INTENT(INOUT) :: &
267 qci
268 REAL, DIMENSION( its:ite , kts:kte, 3 ), &
269 INTENT(INOUT) :: &
270 qrs, &
271 ncr
272 REAL, DIMENSION( ims:ime , kms:kme ), &
273 INTENT(INOUT) :: &
274 q
275 REAL, DIMENSION( ims:ime , kms:kme ), &
276 INTENT(IN ) :: &
277 den, &
278 p, &
279 delz
280 REAL, INTENT(IN ) :: delt, &
281 g, &
282 cpd, &
283 cpv, &
284 ccn0, &
285 t0c, &
286 den0, &
287 rd, &
288 rv, &
289 ep1, &
290 ep2, &
291 qmin, &
292 XLS, &
293 XLV0, &
294 XLF0, &
295 cliq, &
296 cice, &
297 psat, &
298 denr
299 REAL, DIMENSION( ims:ime ), &
300 INTENT(INOUT) :: rain, &
301 rainncv, &
302 sr
303 REAL, DIMENSION( ims:ime ), OPTIONAL, &
304 INTENT(INOUT) :: snow, &
305 snowncv
306 REAL, DIMENSION( ims:ime ), OPTIONAL, &
307 INTENT(INOUT) :: graupel, &
308 graupelncv
309 ! LOCAL VAR
310 REAL, DIMENSION( its:ite , kts:kte , 3) :: &
311 rh, qs, rslope, rslope2, rslope3, rslopeb, &
312 falk, fall, work1, qrs_tmp
313 REAL, DIMENSION( its:ite , kts:kte ) :: &
314 rslopec, rslopec2,rslopec3
315 REAL, DIMENSION( its:ite , kts:kte, 2) :: &
316 avedia
317 REAL, DIMENSION( its:ite , kts:kte ) :: &
318 workn,falln,falkn
319 REAL, DIMENSION( its:ite , kts:kte ) :: &
320 worka,workr
321 REAL, DIMENSION( its:ite , kts:kte ) :: &
322 den_tmp, delz_tmp, ncr_tmp
323 REAL, DIMENSION( its:ite , kts:kte ) :: &
324 falkc, work1c, work2c, fallc
325 REAL, DIMENSION( its:ite , kts:kte ) :: &
326 pcact, prevp, psdep, pgdep, praut, psaut, pgaut, &
327 pracw, psacw, pgacw, pgacr, pgacs, psaci, pgmlt, praci, &
328 piacr, pracs, psacr, pgaci, pseml, pgeml, prevp_s
329 REAL, DIMENSION( its:ite , kts:kte ) :: paacw
330 REAL, DIMENSION( its:ite , kts:kte ) :: &
331 nraut, nracw, nrevp, ncevp, nccol, nrcol, &
332 nsacw, ngacw, niacr, nsacr, ngacr, naacw, &
333 nseml, ngeml, ncact
334 REAL, DIMENSION( its:ite , kts:kte ) :: &
335 pigen, pidep, pcond, xl, cpm, work2, psmlt, psevp, &
336 denfac, xni, pgevp,n0sfac, qsum, &
337 denqrs1, denqr1, denqrs2, denqrs3, denncr3, denqci
338 REAL, DIMENSION( its:ite ) :: &
339 delqrs1, delqrs2, delqrs3, delncr3, delqi
340 REAL :: gfac, sfac
341 ! variables for optimization
342 REAL, DIMENSION( its:ite ) :: tvec1
343 REAL :: temp
344 INTEGER, DIMENSION( its:ite ) :: mnstep, numndt
345 INTEGER, DIMENSION( its:ite ) :: mstep, numdt
346 LOGICAL, DIMENSION( its:ite ) :: flgcld
347 REAL :: &
348 cpmcal, xlcal, lamdac, &
349 diffus, &
350 viscos, xka, venfac, conden, diffac, &
351 x, y, z, a, b, c, d, e, &
352 ndt, qdt, holdrr, holdrs, holdrg, supcol, supcolt, &
353 pvt, coeres, supsat, dtcld, xmi, eacrs, satdt, &
354 qimax, diameter, xni0, roqi0, &
355 fallsum, fallsum_qsi, fallsum_qg, &
356 vt2i,vt2r,vt2s,vt2g,acrfac,egs,egi, &
357 xlwork2, factor, source, value, coecol, &
358 nfrzdtr, nfrzdtc, &
359 taucon, lencon, lenconcr, &
360 xlf, pfrzdtc, pfrzdtr, supice, alpha2, delta2, delta3
361 REAL :: vt2ave
362 REAL :: holdc, holdci
363 !
364 INTEGER :: i, j, k, mstepmax, &
365 iprt, latd, lond, loop, loops, ifsat, n, idim, kdim
366 ! Temporaries used for inlining fpvs function
367 REAL :: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
368 !
369 !=================================================================
370 ! compute internal functions
371 !
372 cpmcal(x) = cpd*(1.-max(x,qmin))+max(x,qmin)*cpv
373 xlcal(x) = xlv0-xlv1*(x-t0c)
374 !----------------------------------------------------------------
375 ! size distributions: (x=mixing ratio, y=air density):
376 ! valid for mixing ratio > 1.e-9 kg/kg.
377 !
378 ! Optimizatin : A**B => exp(log(A)*(B))
379 lamdac(x,y,z)= exp(log(((pidnc*z)/(x*y)))*((.33333333)))
380 !----------------------------------------------------------------
381 ! diffus: diffusion coefficient of the water vapor
382 ! viscos: kinematic viscosity(m2s-1)
383 !
384 diffus(x,y) = 8.794e-5 * exp(log(x)*(1.81)) / y ! 8.794e-5*x**1.81/y
385 viscos(x,y) = 1.496e-6 * (x*sqrt(x)) /(x+120.)/y ! 1.496e-6*x**1.5/(x+120.)/y
386 xka(x,y) = 1.414e3*viscos(x,y)*y
387 diffac(a,b,c,d,e) = d*a*a/(xka(c,d)*rv*c*c)+1./(e*diffus(c,b))
388 venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333))) &
389 /sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
390 conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+d*d/(rv*e)*c/(a*a))
391 !
392 idim = ite-its+1
393 kdim = kte-kts+1
394 !
395 !----------------------------------------------------------------
396 ! paddint 0 for negative values generated by dynamics
397 !
398 do k = kts, kte
399 do i = its, ite
400 qci(i,k,1) = max(qci(i,k,1),0.0)
401 qrs(i,k,1) = max(qrs(i,k,1),0.0)
402 qci(i,k,2) = max(qci(i,k,2),0.0)
403 qrs(i,k,2) = max(qrs(i,k,2),0.0)
404 qrs(i,k,3) = max(qrs(i,k,3),0.0)
405 ncr(i,k,1) = max(ncr(i,k,1),0.0)
406 ncr(i,k,2) = max(ncr(i,k,2),0.0)
407 ncr(i,k,3) = max(ncr(i,k,3),0.0)
408 enddo
409 enddo
410 !
411 !----------------------------------------------------------------
412 ! latent heat for phase changes and heat capacity. neglect the
413 ! changes during microphysical process calculation
414 ! emanuel(1994)
415 !
416 do k = kts, kte
417 do i = its, ite
418 cpm(i,k) = cpmcal(q(i,k))
419 xl(i,k) = xlcal(t(i,k))
420 enddo
421 enddo
422 do k = kts, kte
423 do i = its, ite
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 CALL VREC( tvec1(its), den(its,k), ite-its+1)
449 do i = its, ite
450 tvec1(i) = tvec1(i)*den0
451 enddo
452 CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
453 enddo
454 !
455 ! Inline expansion for fpvs
456 ! qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
457 ! qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
458 hsub = xls
459 hvap = xlv0
460 cvap = cpv
461 ttp=t0c+0.01
462 dldt=cvap-cliq
463 xa=-dldt/rv
464 xb=xa+hvap/(rv*ttp)
465 dldti=cvap-cice
466 xai=-dldti/rv
467 xbi=xai+hsub/(rv*ttp)
468 do k = kts, kte
469 do i = its, ite
470 tr=ttp/t(i,k)
471 qs(i,k,1)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
472 qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
473 qs(i,k,1) = max(qs(i,k,1),qmin)
474 rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
475 tr=ttp/t(i,k)
476 if(t(i,k).lt.ttp) then
477 qs(i,k,2)=psat*exp(log(tr)*(xai))*exp(xbi*(1.-tr))
478 else
479 qs(i,k,2)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
480 endif
481 qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
482 qs(i,k,2) = max(qs(i,k,2),qmin)
483 rh(i,k,2) = max(q(i,k) / qs(i,k,2),qmin)
484 enddo
485 enddo
486 !
487 !----------------------------------------------------------------
488 ! initialize the variables for microphysical physics
489 !
490 !
491 do k = kts, kte
492 do i = its, ite
493 prevp(i,k) = 0.
494 psdep(i,k) = 0.
495 pgdep(i,k) = 0.
496 praut(i,k) = 0.
497 psaut(i,k) = 0.
498 pgaut(i,k) = 0.
499 pracw(i,k) = 0.
500 praci(i,k) = 0.
501 piacr(i,k) = 0.
502 psaci(i,k) = 0.
503 psacw(i,k) = 0.
504 pracs(i,k) = 0.
505 psacr(i,k) = 0.
506 pgacw(i,k) = 0.
507 paacw(i,k) = 0.
508 pgaci(i,k) = 0.
509 pgacr(i,k) = 0.
510 pgacs(i,k) = 0.
511 pigen(i,k) = 0.
512 pidep(i,k) = 0.
513 pcond(i,k) = 0.
514 psmlt(i,k) = 0.
515 pgmlt(i,k) = 0.
516 pseml(i,k) = 0.
517 pgeml(i,k) = 0.
518 psevp(i,k) = 0.
519 pgevp(i,k) = 0.
520 pcact(i,k) = 0.
521 prevp_s(i,k) = 0.
522 falk(i,k,1) = 0.
523 falk(i,k,2) = 0.
524 falk(i,k,3) = 0.
525 fall(i,k,1) = 0.
526 fall(i,k,2) = 0.
527 fall(i,k,3) = 0.
528 fallc(i,k) = 0.
529 falkc(i,k) = 0.
530 falln(i,k) =0.
531 falkn(i,k) =0.
532 xni(i,k) = 1.e3
533 nsacw(i,k) = 0.
534 ngacw(i,k) = 0.
535 naacw(i,k) = 0.
536 niacr(i,k) = 0.
537 nsacr(i,k) = 0.
538 ngacr(i,k) = 0.
539 nseml(i,k) = 0.
540 ngeml(i,k) = 0.
541 nracw(i,k) = 0.
542 nccol(i,k) = 0.
543 nrcol(i,k) = 0.
544 ncact(i,k) = 0.
545 nraut(i,k) = 0.
546 nrevp(i,k) = 0.
547 ncevp(i,k) = 0.
548 enddo
549 enddo
550 do k = kts, kte
551 do i = its, ite
552 if(qci(i,k,1).le.qmin .or. ncr(i,k,2).le.ncmin ) then
553 rslopec(i,k) = rslopecmax
554 rslopec2(i,k) = rslopec2max
555 rslopec3(i,k) = rslopec3max
556 else
557 rslopec(i,k) = 1./lamdac(qci(i,k,1),den(i,k),ncr(i,k,2))
558 rslopec2(i,k) = rslopec(i,k)*rslopec(i,k)
559 rslopec3(i,k) = rslopec2(i,k)*rslopec(i,k)
560 endif
561 !-------------------------------------------------------------
562 ! Ni: ice crystal number concentraiton [HDC 5c]
563 !-------------------------------------------------------------
564 temp = (den(i,k)*max(qci(i,k,2),qmin))
565 temp = sqrt(sqrt(temp*temp*temp))
566 xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
567 enddo
568 enddo
569 !----------------------------------------------------------------
570 ! compute the fallout term:
571 ! first, vertical terminal velosity for minor loops
572 !----------------------------------------------------------------
573 do k = kts, kte
574 do i = its, ite
575 qrs_tmp(i,k,1) = qrs(i,k,1)
576 qrs_tmp(i,k,2) = qrs(i,k,2)
577 qrs_tmp(i,k,3) = qrs(i,k,3)
578 ncr_tmp(i,k) = ncr(i,k,3)
579 enddo
580 enddo
581 call slope_wdm6(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
582 rslope3,work1,workn,its,ite,kts,kte)
583 !
584 ! vt update for qr and nr
585 mstepmax = 1
586 numdt = 1
587 do k = kte, kts, -1
588 do i = its, ite
589 work1(i,k,1) = work1(i,k,1)/delz(i,k)
590 workn(i,k) = workn(i,k)/delz(i,k)
591 numdt(i) = max(nint(max(work1(i,k,1),workn(i,k))*dtcld+.5),1)
592 if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
593 enddo
594 enddo
595 do i = its, ite
596 if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
597 enddo
598 !
599 do n = 1, mstepmax
600 k = kte
601 do i = its, ite
602 if(n.le.mstep(i)) then
603 falk(i,k,1) = den(i,k)*qrs(i,k,1)*work1(i,k,1)/mstep(i)
604 falkn(i,k) = ncr(i,k,3)*workn(i,k)/mstep(i)
605 fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
606 falln(i,k) = falln(i,k)+falkn(i,k)
607 qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcld/den(i,k),0.)
608 ncr(i,k,3) = max(ncr(i,k,3)-falkn(i,k)*dtcld,0.)
609 endif
610 enddo
611 do k = kte-1, kts, -1
612 do i = its, ite
613 if(n.le.mstep(i)) then
614 falk(i,k,1) = den(i,k)*qrs(i,k,1)*work1(i,k,1)/mstep(i)
615 falkn(i,k) = ncr(i,k,3)*workn(i,k)/mstep(i)
616 fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
617 falln(i,k) = falln(i,k)+falkn(i,k)
618 qrs(i,k,1) = max(qrs(i,k,1)-(falk(i,k,1)-falk(i,k+1,1) &
619 *delz(i,k+1)/delz(i,k))*dtcld/den(i,k),0.)
620 ncr(i,k,3) = max(ncr(i,k,3)-(falkn(i,k)-falkn(i,k+1)*delz(i,k+1) &
621 /delz(i,k))*dtcld,0.)
622 endif
623 enddo
624 enddo
625 do k = kts, kte
626 do i = its, ite
627 qrs_tmp(i,k,1) = qrs(i,k,1)
628 ncr_tmp(i,k) = ncr(i,k,3)
629 enddo
630 enddo
631 call slope_rain(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
632 rslope3,work1,workn,its,ite,kts,kte)
633 do k = kte, kts, -1
634 do i = its, ite
635 work1(i,k,1) = work1(i,k,1)/delz(i,k)
636 workn(i,k) = workn(i,k)/delz(i,k)
637 enddo
638 enddo
639 enddo
640 ! for semi
641 do k = kte, kts, -1
642 do i = its, ite
643 qsum(i,k) = max( (qrs(i,k,2)+qrs(i,k,3)), 1.E-15)
644 if(qsum(i,k) .gt. 1.e-15 ) then
645 worka(i,k) = (work1(i,k,2)*qrs(i,k,2) + work1(i,k,3)*qrs(i,k,3)) &
646 /qsum(i,k)
647 else
648 worka(i,k) = 0.
649 endif
650 denqrs2(i,k) = den(i,k)*qrs(i,k,2)
651 denqrs3(i,k) = den(i,k)*qrs(i,k,3)
652 enddo
653 enddo
654 call nislfv_rain_plm6(idim,kdim,den_tmp,denfac,t,delz_tmp,worka, &
655 denqrs2,denqrs3,delqrs2,delqrs3,dtcld,1,1)
656 do k = kts, kte
657 do i = its, ite
658 qrs(i,k,2) = max(denqrs2(i,k)/den(i,k),0.)
659 qrs(i,k,3) = max(denqrs3(i,k)/den(i,k),0.)
660 fall(i,k,2) = denqrs2(i,k)*worka(i,k)/delz(i,k)
661 fall(i,k,3) = denqrs3(i,k)*worka(i,k)/delz(i,k)
662 enddo
663 enddo
664 do i = its, ite
665 fall(i,1,2) = delqrs2(i)/delz(i,1)/dtcld
666 fall(i,1,3) = delqrs3(i)/delz(i,1)/dtcld
667 enddo
668 do k = kts, kte
669 do i = its, ite
670 qrs_tmp(i,k,1) = qrs(i,k,1)
671 qrs_tmp(i,k,2) = qrs(i,k,2)
672 qrs_tmp(i,k,3) = qrs(i,k,3)
673 ncr_tmp(i,k) = ncr(i,k,3)
674 enddo
675 enddo
676 call slope_wdm6(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
677 rslope3,work1,workn,its,ite,kts,kte)
678 !
679 do k = kte, kts, -1
680 do i = its, ite
681 supcol = t0c-t(i,k)
682 n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
683 if(t(i,k).gt.t0c) then
684 !---------------------------------------------------------------
685 ! psmlt: melting of snow [HL A33] [RH83 A25]
686 ! (T>T0: QS->QR)
687 !---------------------------------------------------------------
688 xlf = xlf0
689 work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
690 if(qrs(i,k,2).gt.0.) then
691 coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
692 psmlt(i,k) = xka(t(i,k),den(i,k))/xlf*(t0c-t(i,k))*pi/2. &
693 *n0sfac(i,k)*(precs1*rslope2(i,k,2) &
694 +precs2*work2(i,k)*coeres)
695 psmlt(i,k) = min(max(psmlt(i,k)*dtcld/mstep(i),-qrs(i,k,2) &
696 /mstep(i)),0.)
697 qrs(i,k,2) = qrs(i,k,2) + psmlt(i,k)
698 qrs(i,k,1) = qrs(i,k,1) - psmlt(i,k)
699 !-------------------------------------------------------------------
700 ! nsmlt: melting of snow [LH A27]
701 ! (T>T0: ->NR)
702 !-------------------------------------------------------------------
703 if(qrs(i,k,2).gt.qcrmin) then
704 sfac = rslope(i,k,2)*n0s*n0sfac(i,k)/qrs(i,k,2)
705 ncr(i,k,3) = ncr(i,k,3) - sfac*psmlt(i,k)
706 endif
707 t(i,k) = t(i,k) + xlf/cpm(i,k)*psmlt(i,k)
708 endif
709 !---------------------------------------------------------------
710 ! pgmlt: melting of graupel [HL A23] [LFO 47]
711 ! (T>T0: QG->QR)
712 !---------------------------------------------------------------
713 if(qrs(i,k,3).gt.0.) then
714 coeres = rslope2(i,k,3)*sqrt(rslope(i,k,3)*rslopeb(i,k,3))
715 pgmlt(i,k) = xka(t(i,k),den(i,k))/xlf*(t0c-t(i,k))*(precg1 &
716 *rslope2(i,k,3) + precg2*work2(i,k)*coeres)
717 pgmlt(i,k) = min(max(pgmlt(i,k)*dtcld/mstep(i), &
718 -qrs(i,k,3)/mstep(i)),0.)
719 qrs(i,k,3) = qrs(i,k,3) + pgmlt(i,k)
720 qrs(i,k,1) = qrs(i,k,1) - pgmlt(i,k)
721 !-------------------------------------------------------------------
722 ! ngmlt: melting of graupel [LH A28]
723 ! (T>T0: ->NR)
724 !-------------------------------------------------------------------
725 if(qrs(i,k,3).gt.qcrmin) then
726 gfac = rslope(i,k,3)*n0g/qrs(i,k,3)
727 ncr(i,k,3) = ncr(i,k,3) - gfac*pgmlt(i,k)
728 endif
729 t(i,k) = t(i,k) + xlf/cpm(i,k)*pgmlt(i,k)
730 endif
731 endif
732 enddo
733 enddo
734 !---------------------------------------------------------------
735 ! Vice [ms-1] : fallout of ice crystal [HDC 5a]
736 !---------------------------------------------------------------
737 do k = kte, kts, -1
738 do i = its, ite
739 if(qci(i,k,2).le.0.) then
740 work1c(i,k) = 0.
741 else
742 xmi = den(i,k)*qci(i,k,2)/xni(i,k)
743 diameter = max(min(dicon * sqrt(xmi),dimax), 1.e-25)
744 work1c(i,k) = 1.49e4*exp(log(diameter)*(1.31))
745 endif
746 enddo
747 enddo
748 !
749 ! forward semi-laglangian scheme (JH), PCM (piecewise constant), (linear)
750 !
751 do k = kte, kts, -1
752 do i = its, ite
753 denqci(i,k) = den(i,k)*qci(i,k,2)
754 enddo
755 enddo
756 call nislfv_rain_plmr(idim,kdim,den_tmp,denfac,t,delz_tmp,work1c,denqci,denqci, &
757 delqi,dtcld,1,0,0)
758 do k = kts, kte
759 do i = its, ite
760 qci(i,k,2) = max(denqci(i,k)/den(i,k),0.)
761 enddo
762 enddo
763 do i = its, ite
764 fallc(i,1) = delqi(i)/delz(i,1)/dtcld
765 enddo
766 !
767 !----------------------------------------------------------------
768 ! rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
769 !
770 do i = its, ite
771 fallsum = fall(i,kts,1)+fall(i,kts,2)+fall(i,kts,3)+fallc(i,kts)
772 fallsum_qsi = fall(i,kts,2)+fallc(i,kts)
773 fallsum_qg = fall(i,kts,3)
774 rainncv(i) = 0.
775 if(fallsum.gt.0.) then
776 rainncv(i) = fallsum*delz(i,kts)/denr*dtcld*1000.
777 rain(i) = fallsum*delz(i,kts)/denr*dtcld*1000. + rain(i)
778 endif
779 IF ( PRESENT (snowncv) .AND. PRESENT (snow)) THEN
780 snowncv(i) = 0.
781 if(fallsum_qsi.gt.0.) then
782 snowncv(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000.
783 snow(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. + snow(i)
784 endif
785 ENDIF
786 IF ( PRESENT (graupelncv) .AND. PRESENT (graupel)) THEN
787 graupelncv(i) = 0.
788 if(fallsum_qg.gt.0.) then
789 graupelncv(i) = fallsum_qg*delz(i,kts)/denr*dtcld*1000.
790 graupel(i) = fallsum_qg*delz(i,kts)/denr*dtcld*1000. + graupel(i)
791 endif
792 ENDIF
793 sr(i) = 0.
794 if(fallsum.gt.0.)sr(i)=(fallsum_qsi*delz(i,kts)/denr*dtcld*1000. + &
795 fallsum_qg*delz(i,kts)/denr*dtcld*1000.) &
796 /(rainncv(i)+1.e-12)
797 enddo
798 !
799 !---------------------------------------------------------------
800 ! pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
801 ! (T>T0: QI->QC)
802 !---------------------------------------------------------------
803 do k = kts, kte
804 do i = its, ite
805 supcol = t0c-t(i,k)
806 xlf = xls-xl(i,k)
807 if(supcol.lt.0.) xlf = xlf0
808 if(supcol.lt.0 .and. qci(i,k,2).gt.0.) then
809 qci(i,k,1) = qci(i,k,1) + qci(i,k,2)
810 !---------------------------------------------------------------
811 ! nimlt: instantaneous melting of cloud ice [LH A18]
812 ! (T>T0: ->NC)
813 !--------------------------------------------------------------
814 ncr(i,k,2) = ncr(i,k,2) + xni(i,k)
815 t(i,k) = t(i,k) - xlf/cpm(i,k)*qci(i,k,2)
816 qci(i,k,2) = 0.
817 endif
818 !---------------------------------------------------------------
819 ! pihmf: homogeneous of cloud water below -40c [HL A45]
820 ! (T<-40C: QC->QI)
821 !---------------------------------------------------------------
822 if(supcol.gt.40. .and. qci(i,k,1).gt.0.) then
823 qci(i,k,2) = qci(i,k,2) + qci(i,k,1)
824 !---------------------------------------------------------------
825 ! nihmf: homogeneous of cloud water below -40c [LH A17]
826 ! (T<-40C: NC->)
827 !---------------------------------------------------------------
828 if(ncr(i,k,2).gt.0.) ncr(i,k,2) = 0.
829 t(i,k) = t(i,k) + xlf/cpm(i,k)*qci(i,k,1)
830 qci(i,k,1) = 0.
831 endif
832 !---------------------------------------------------------------
833 ! pihtf: heterogeneous of cloud water [HL A44]
834 ! (T0>T>-40C: QC->QI)
835 !---------------------------------------------------------------
836 if(supcol.gt.0. .and. qci(i,k,1).gt.qmin) then
837 supcolt=min(supcol,70.)
838 pfrzdtc = min(pi*pi*pfrz1*(exp(pfrz2*supcolt)-1.)*denr/den(i,k) &
839 *ncr(i,k,2)*rslopec3(i,k)*rslopec3(i,k)/18.*dtcld &
840 ,qci(i,k,1))
841 !---------------------------------------------------------------
842 ! nihtf: heterogeneous of cloud water [LH A16]
843 ! (T0>T>-40C: NC->)
844 !---------------------------------------------------------------
845 if(ncr(i,k,2).gt.ncmin) then
846 nfrzdtc = min(pi*pfrz1*(exp(pfrz2*supcolt)-1.)*ncr(i,k,2) &
847 *rslopec3(i,k)/6.*dtcld,ncr(i,k,2))
848 ncr(i,k,2) = ncr(i,k,2) - nfrzdtc
849 endif
850 qci(i,k,2) = qci(i,k,2) + pfrzdtc
851 t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtc
852 qci(i,k,1) = qci(i,k,1)-pfrzdtc
853 endif
854 !---------------------------------------------------------------
855 ! pgfrz: freezing of rain water [HL A20] [LFO 45]
856 ! (T<T0, QR->QG)
857 !---------------------------------------------------------------
858 if(supcol.gt.0. .and. qrs(i,k,1).gt.0.) then
859 supcolt=min(supcol,70.)
860 pfrzdtr = min(140.*(pi*pi)*pfrz1*ncr(i,k,3)*denr/den(i,k) &
861 *(exp(pfrz2*supcolt)-1.)*rslope3(i,k,1)*rslope3(i,k,1) &
862 *dtcld,qrs(i,k,1))
863 !---------------------------------------------------------------
864 ! ngfrz: freezing of rain water [LH A26]
865 ! (T<T0, NR-> )
866 !---------------------------------------------------------------
867 if(ncr(i,k,3).gt.nrmin) then
868 nfrzdtr = min(4.*pi*pfrz1*ncr(i,k,3)*(exp(pfrz2*supcolt)-1.) &
869 *rslope3(i,k,1)*dtcld, ncr(i,k,3))
870 ncr(i,k,3) = ncr(i,k,3) - nfrzdtr
871 endif
872 qrs(i,k,3) = qrs(i,k,3) + pfrzdtr
873 t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtr
874 qrs(i,k,1) = qrs(i,k,1) - pfrzdtr
875 endif
876 enddo
877 enddo
878 !
879 do k = kts, kte
880 do i = its, ite
881 ncr(i,k,2) = max(ncr(i,k,2),0.0)
882 ncr(i,k,3) = max(ncr(i,k,3),0.0)
883 enddo
884 enddo
885 !
886 !----------------------------------------------------------------
887 ! update the slope parameters for microphysics computation
888 !
889 do k = kts, kte
890 do i = its, ite
891 qrs_tmp(i,k,1) = qrs(i,k,1)
892 qrs_tmp(i,k,2) = qrs(i,k,2)
893 qrs_tmp(i,k,3) = qrs(i,k,3)
894 ncr_tmp(i,k) = ncr(i,k,3)
895 enddo
896 enddo
897 call slope_wdm6(qrs_tmp,ncr_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2, &
898 rslope3,work1,workn,its,ite,kts,kte)
899 do k = kts, kte
900 do i = its, ite
901 !-----------------------------------------------------------------
902 ! compute the mean-volume drop diameter [LH A10]
903 ! for raindrop distribution
904 !-----------------------------------------------------------------
905 avedia(i,k,2) = rslope(i,k,1)*((24.)**(.3333333))
906 !
907 if(qci(i,k,1).le.qmin .or. ncr(i,k,2).le.ncmin) then
908 rslopec(i,k) = rslopecmax
909 rslopec2(i,k) = rslopec2max
910 rslopec3(i,k) = rslopec3max
911 else
912 rslopec(i,k) = 1./lamdac(qci(i,k,1),den(i,k),ncr(i,k,2))
913 rslopec2(i,k) = rslopec(i,k)*rslopec(i,k)
914 rslopec3(i,k) = rslopec2(i,k)*rslopec(i,k)
915 endif
916 !--------------------------------------------------------------------
917 ! compute the mean-volume drop diameter [LH A7]
918 ! for cloud-droplet distribution
919 !--------------------------------------------------------------------
920 avedia(i,k,1) = rslopec(i,k)
921 enddo
922 enddo
923 !
924 do k = kts, kte
925 do i = its, ite
926 work1(i,k,1) = diffac(xl(i,k),p(i,k),t(i,k),den(i,k),qs(i,k,1))
927 work1(i,k,2) = diffac(xls,p(i,k),t(i,k),den(i,k),qs(i,k,2))
928 work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
929 enddo
930 enddo
931 !
932 !===============================================================
933 !
934 ! warm rain processes
935 !
936 ! - follows the double-moment processes in Lim and Hong
937 !
938 !===============================================================
939 !
940 do k = kts, kte
941 do i = its, ite
942 supsat = max(q(i,k),qmin)-qs(i,k,1)
943 satdt = supsat/dtcld
944 !---------------------------------------------------------------
945 ! praut: auto conversion rate from cloud to rain [LH 9] [CP 17]
946 ! (QC->QR)
947 !--------------------------------------------------------------
948 lencon = 2.7e-2*den(i,k)*qci(i,k,1)*(1.e20/16.*rslopec2(i,k) &
949 *rslopec2(i,k)-0.4)
950 lenconcr = max(1.2*lencon, qcrmin)
951 if(avedia(i,k,1).gt.di15) then
952 taucon = 3.7/den(i,k)/qci(i,k,1)/(0.5e6*rslopec(i,k)-7.5)
953 praut(i,k) = lencon/taucon
954 praut(i,k) = min(max(praut(i,k),0.),qci(i,k,1)/dtcld)
955 !---------------------------------------------------------------
956 ! nraut: auto conversion rate from cloud to rain [LH A6] [CP 18 & 19]
957 ! (NC->NR)
958 !---------------------------------------------------------------
959 nraut(i,k) = 3.5e9*den(i,k)*praut(i,k)
960 if(qrs(i,k,1).gt.lenconcr) &
961 nraut(i,k) = ncr(i,k,3)/qrs(i,k,1)*praut(i,k)
962 nraut(i,k) = min(nraut(i,k),ncr(i,k,2)/dtcld)
963 endif
964 !---------------------------------------------------------------
965 ! pracw: accretion of cloud water by rain [LH 10] [CP 22 & 23]
966 ! (QC->QR)
967 ! nracw: accretion of cloud water by rain [LH A9]
968 ! (NC->)
969 !---------------------------------------------------------------
970 if(qrs(i,k,1).ge.lenconcr) then
971 if(avedia(i,k,2).ge.di100) then
972 nracw(i,k) = min(ncrk1*ncr(i,k,2)*ncr(i,k,3)*(rslopec3(i,k) &
973 + 24.*rslope3(i,k,1)),ncr(i,k,2)/dtcld)
974 pracw(i,k) = min(pi/6.*(denr/den(i,k))*ncrk1*ncr(i,k,2) &
975 *ncr(i,k,3)*rslopec3(i,k)*(2.*rslopec3(i,k) &
976 + 24.*rslope3(i,k,1)),qci(i,k,1)/dtcld)
977 else
978 nracw(i,k) = min(ncrk2*ncr(i,k,2)*ncr(i,k,3)*(2.*rslopec3(i,k) &
979 *rslopec3(i,k)+5040.*rslope3(i,k,1) &
980 *rslope3(i,k,1)),ncr(i,k,2)/dtcld)
981 pracw(i,k) = min(pi/6.*(denr/den(i,k))*ncrk2*ncr(i,k,2) &
982 *ncr(i,k,3)*rslopec3(i,k)*(6.*rslopec3(i,k) &
983 *rslopec3(i,k)+5040.*rslope3(i,k,1)*rslope3(i,k,1)) &
984 ,qci(i,k,1)/dtcld)
985 endif
986 endif
987 !----------------------------------------------------------------
988 ! nccol: self collection of cloud water [LH A8] [CP 24 & 25]
989 ! (NC->)
990 !----------------------------------------------------------------
991 if(avedia(i,k,1).ge.di100) then
992 nccol(i,k) = ncrk1*ncr(i,k,2)*ncr(i,k,2)*rslopec3(i,k)
993 else
994 nccol(i,k) = 2.*ncrk2*ncr(i,k,2)*ncr(i,k,2)*rslopec3(i,k) &
995 *rslopec3(i,k)
996 endif
997 !----------------------------------------------------------------
998 ! nrcol: self collection of rain-drops and break-up [LH A21] [CP 24 & 25]
999 ! (NR->)
1000 !----------------------------------------------------------------
1001 if(qrs(i,k,1).ge.lenconcr) then
1002 if(avedia(i,k,2).lt.di100) then
1003 nrcol(i,k) = 5040.*ncrk2*ncr(i,k,3)*ncr(i,k,3)*rslope3(i,k,1) &
1004 *rslope3(i,k,1)
1005 elseif(avedia(i,k,2).ge.di100 .and. avedia(i,k,2).lt.di600) then
1006 nrcol(i,k) = 24.*ncrk1*ncr(i,k,3)*ncr(i,k,3)*rslope3(i,k,1)
1007 elseif(avedia(i,k,2).ge.di600 .and. avedia(i,k,2).lt.di2000) then
1008 coecol = -2.5e3*(avedia(i,k,2)-di600)
1009 nrcol(i,k) = 24.*exp(coecol)*ncrk1*ncr(i,k,3)*ncr(i,k,3) &
1010 *rslope3(i,k,1)
1011 else
1012 nrcol(i,k) = 0.
1013 endif
1014 endif
1015 !---------------------------------------------------------------
1016 ! prevp: evaporation/condensation rate of rain [HL A41]
1017 ! (QV->QR or QR->QV)
1018 !---------------------------------------------------------------
1019 if(qrs(i,k,1).gt.0.) then
1020 coeres = rslope(i,k,1)*sqrt(rslope(i,k,1)*rslopeb(i,k,1))
1021 prevp(i,k) = (rh(i,k,1)-1.)*ncr(i,k,3)*(precr1*rslope(i,k,1) &
1022 + precr2*work2(i,k)*coeres)/work1(i,k,1)
1023 if(prevp(i,k).lt.0.) then
1024 prevp(i,k) = max(prevp(i,k),-qrs(i,k,1)/dtcld)
1025 prevp(i,k) = max(prevp(i,k),satdt/2)
1026 !----------------------------------------------------------------
1027 ! Nrevp: evaporation/condensation rate of rain [LH A14]
1028 ! (NR->NC)
1029 !----------------------------------------------------------------
1030 if(avedia(i,k,2).le.di82) then
1031 nrevp(i,k) = ncr(i,k,3)/dtcld
1032 !----------------------------------------------------------------
1033 ! Prevp_s: evaporation/condensation rate of rain [LH A15] [KK 23]
1034 ! (QR->QC)
1035 !----------------------------------------------------------------
1036 prevp_s(i,k) = qrs(i,k,1)/dtcld
1037 endif
1038 else
1039 !
1040 prevp(i,k) = min(prevp(i,k),satdt/2)
1041 endif
1042 endif
1043 enddo
1044 enddo
1045 !
1046 !===============================================================
1047 !
1048 ! cold rain processes
1049 !
1050 ! - follows the revised ice microphysics processes in HDC
1051 ! - the processes same as in RH83 and RH84 and LFO behave
1052 ! following ice crystal hapits defined in HDC, inclduing
1053 ! intercept parameter for snow (n0s), ice crystal number
1054 ! concentration (ni), ice nuclei number concentration
1055 ! (n0i), ice diameter (d)
1056 !
1057 !===============================================================
1058 !
1059 do k = kts, kte
1060 do i = its, ite
1061 supcol = t0c-t(i,k)
1062 n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
1063 supsat = max(q(i,k),qmin)-qs(i,k,2)
1064 satdt = supsat/dtcld
1065 ifsat = 0
1066 !-------------------------------------------------------------
1067 ! Ni: ice crystal number concentraiton [HDC 5c]
1068 !-------------------------------------------------------------
1069 ! xni(i,k) = min(max(5.38e7*(den(i,k) &
1070 ! *max(qci(i,k,2),qmin))**0.75,1.e3),1.e6)
1071 temp = (den(i,k)*max(qci(i,k,2),qmin))
1072 temp = sqrt(sqrt(temp*temp*temp))
1073 xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
1074 eacrs = exp(0.07*(-supcol))
1075 !
1076 xmi = den(i,k)*qci(i,k,2)/xni(i,k)
1077 diameter = min(dicon * sqrt(xmi),dimax)
1078 vt2i = 1.49e4*diameter**1.31
1079 vt2r=pvtr*rslopeb(i,k,1)*denfac(i,k)
1080 vt2s=pvts*rslopeb(i,k,2)*denfac(i,k)
1081 vt2g=pvtg*rslopeb(i,k,3)*denfac(i,k)
1082 qsum(i,k) = max((qrs(i,k,2)+qrs(i,k,3)),1.e-15)
1083 if(qsum(i,k) .gt. 1.e-15) then
1084 vt2ave=(vt2s*qrs(i,k,2)+vt2g*qrs(i,k,3))/(qsum(i,k))
1085 else
1086 vt2ave=0.
1087 endif
1088 if(supcol.gt.0. .and. qci(i,k,2).gt.qmin) then
1089 if(qrs(i,k,1).gt.qcrmin) then
1090 !-------------------------------------------------------------
1091 ! praci: Accretion of cloud ice by rain [HL A15] [LFO 25]
1092 ! (T<T0: QI->QR)
1093 !-------------------------------------------------------------
1094 acrfac = 6.*rslope2(i,k,1)+4.*diameter*rslope(i,k,1) + diameter**2
1095 praci(i,k) = pi*qci(i,k,2)*ncr(i,k,3)*abs(vt2r-vt2i)*acrfac/4.
1096 praci(i,k) = min(praci(i,k),qci(i,k,2)/dtcld)
1097 !-------------------------------------------------------------
1098 ! piacr: Accretion of rain by cloud ice [HL A19] [LFO 26]
1099 ! (T<T0: QR->QS or QR->QG)
1100 !-------------------------------------------------------------
1101 piacr(i,k) = pi*pi*avtr*ncr(i,k,3)*denr*xni(i,k)*denfac(i,k) &
1102 *g7pbr*rslope3(i,k,1)*rslope2(i,k,1)*rslopeb(i,k,1) &
1103 /24./den(i,k)
1104 piacr(i,k) = min(piacr(i,k),qrs(i,k,1)/dtcld)
1105 endif
1106 !-------------------------------------------------------------
1107 ! niacr: Accretion of rain by cloud ice [LH A25]
1108 ! (T<T0: NR->)
1109 !-------------------------------------------------------------
1110 if(ncr(i,k,3).gt.nrmin) then
1111 niacr(i,k) = pi*avtr*ncr(i,k,3)*xni(i,k)*denfac(i,k)*g4pbr &
1112 *rslope2(i,k,1)*rslopeb(i,k,1)/4.
1113 niacr(i,k) = min(niacr(i,k),ncr(i,k,3)/dtcld)
1114 endif
1115 !-------------------------------------------------------------
1116 ! psaci: Accretion of cloud ice by snow [HDC 10]
1117 ! (T<T0: QI->QS)
1118 !-------------------------------------------------------------
1119 if(qrs(i,k,2).gt.qcrmin) then
1120 acrfac = 2.*rslope3(i,k,2)+2.*diameter*rslope2(i,k,2) &
1121 + diameter**2*rslope(i,k,2)
1122 psaci(i,k) = pi*qci(i,k,2)*eacrs*n0s*n0sfac(i,k) &
1123 *abs(vt2ave-vt2i)*acrfac/4.
1124 psaci(i,k) = min(psaci(i,k),qci(i,k,2)/dtcld)
1125 endif
1126 !-------------------------------------------------------------
1127 ! pgaci: Accretion of cloud ice by graupel [HL A17] [LFO 41]
1128 ! (T<T0: QI->QG)
1129 !-------------------------------------------------------------
1130 if(qrs(i,k,3).gt.qcrmin) then
1131 egi = exp(0.07*(-supcol))
1132 acrfac = 2.*rslope3(i,k,3)+2.*diameter*rslope2(i,k,3) &
1133 + diameter**2*rslope(i,k,3)
1134 pgaci(i,k) = pi*egi*qci(i,k,2)*n0g*abs(vt2ave-vt2i)*acrfac/4.
1135 pgaci(i,k) = min(pgaci(i,k),qci(i,k,2)/dtcld)
1136 endif
1137 endif
1138 !-------------------------------------------------------------
1139 ! psacw: Accretion of cloud water by snow [HL A7] [LFO 24]
1140 ! (T<T0: QC->QS, and T>=T0: QC->QR)
1141 !-------------------------------------------------------------
1142 if(qrs(i,k,2).gt.qcrmin .and. qci(i,k,1).gt.qmin) then
1143 psacw(i,k) = min(pacrc*n0sfac(i,k)*rslope3(i,k,2)*rslopeb(i,k,2) &
1144 *qci(i,k,1)*denfac(i,k),qci(i,k,1)/dtcld)
1145 endif
1146 !-------------------------------------------------------------
1147 ! nsacw: Accretion of cloud water by snow [LH A12]
1148 ! (NC ->)
1149 !-------------------------------------------------------------
1150 if(qrs(i,k,2).gt.qcrmin .and. ncr(i,k,2).gt.ncmin) then
1151 nsacw(i,k) = min(pacrc*n0sfac(i,k)*rslope3(i,k,2)*rslopeb(i,k,2) &
1152 *ncr(i,k,2)*denfac(i,k),ncr(i,k,2)/dtcld)
1153 endif
1154 !-------------------------------------------------------------
1155 ! pgacw: Accretion of cloud water by graupel [HL A6] [LFO 40]
1156 ! (T<T0: QC->QG, and T>=T0: QC->QR)
1157 !-------------------------------------------------------------
1158 if(qrs(i,k,3).gt.qcrmin .and. qci(i,k,1).gt.qmin) then
1159 pgacw(i,k) = min(pacrg*rslope3(i,k,3)*rslopeb(i,k,3)*qci(i,k,1) &
1160 *denfac(i,k),qci(i,k,1)/dtcld)
1161 endif
1162 !-------------------------------------------------------------
1163 ! ngacw: Accretion of cloud water by graupel [LH A13]
1164 ! (NC->
1165 !-------------------------------------------------------------
1166 if(qrs(i,k,3).gt.qcrmin .and. ncr(i,k,2).gt.ncmin) then
1167 ngacw(i,k) = min(pacrg*rslope3(i,k,3)*rslopeb(i,k,3)*ncr(i,k,2) &
1168 *denfac(i,k),ncr(i,k,2)/dtcld)
1169 endif
1170 !-------------------------------------------------------------
1171 ! paacw: Accretion of cloud water by averaged snow/graupel
1172 ! (T<T0: QC->QG or QS, and T>=T0: QC->QR)
1173 !-------------------------------------------------------------
1174 if(qrs(i,k,2).gt.qcrmin .and. qrs(i,k,3).gt.qcrmin) then
1175 paacw(i,k) = (qrs(i,k,2)*psacw(i,k)+qrs(i,k,3)*pgacw(i,k))/(qsum(i,k))
1176 !-------------------------------------------------------------
1177 ! naacw: Accretion of cloud water by averaged snow/graupel
1178 ! (Nc->)
1179 !-------------------------------------------------------------
1180 naacw(i,k) = (qrs(i,k,2)*nsacw(i,k)+qrs(i,k,3)*ngacw(i,k))/(qsum(i,k))
1181 endif
1182 !-------------------------------------------------------------
1183 ! pracs: Accretion of snow by rain [HL A11] [LFO 27]
1184 ! (T<T0: QS->QG)
1185 !-------------------------------------------------------------
1186 if(qrs(i,k,2).gt.qcrmin .and. qrs(i,k,1).gt.qcrmin) then
1187 if(supcol.gt.0) then
1188 acrfac = 5.*rslope3(i,k,2)*rslope3(i,k,2) &
1189 + 4.*rslope3(i,k,2)*rslope2(i,k,2)*rslope(i,k,1) &
1190 + 1.5*rslope2(i,k,2)*rslope2(i,k,2)*rslope2(i,k,1)
1191 pracs(i,k) = pi*pi*ncr(i,k,3)*n0s*n0sfac(i,k)*abs(vt2r-vt2ave) &
1192 *(dens/den(i,k))*acrfac
1193 pracs(i,k) = min(pracs(i,k),qrs(i,k,2)/dtcld)
1194 endif
1195 !-------------------------------------------------------------
1196 ! psacr: Accretion of rain by snow [HL A10] [LFO 28]
1197 ! (T<T0:QR->QS or QR->QG) (T>=T0: enhance melting of snow)
1198 !-------------------------------------------------------------
1199 acrfac = 30.*rslope3(i,k,1)*rslope2(i,k,1)*rslope(i,k,2) &
1200 + 5.*rslope2(i,k,1)*rslope2(i,k,1)*rslope2(i,k,2) &
1201 + 2.*rslope3(i,k,1)*rslope3(i,k,2)
1202 psacr(i,k) = pi*pi*ncr(i,k,3)*n0s*n0sfac(i,k)*abs(vt2ave-vt2r) &
1203 *(denr/den(i,k))*acrfac
1204 psacr(i,k) = min(psacr(i,k),qrs(i,k,1)/dtcld)
1205 endif
1206 if(qrs(i,k,2).gt.qcrmin .and. ncr(i,k,3).gt.nrmin) then
1207 !-------------------------------------------------------------
1208 ! nsacr: Accretion of rain by snow [LH A23]
1209 ! (T<T0:NR->)
1210 !-------------------------------------------------------------
1211 acrfac = 1.5*rslope2(i,k,1)*rslope(i,k,2) &
1212 + 1.0*rslope(i,k,1)*rslope2(i,k,2)+.5*rslope3(i,k,2)
1213 nsacr(i,k) = pi*ncr(i,k,3)*n0s*n0sfac(i,k)*abs(vt2ave-vt2r) &
1214 *acrfac
1215 nsacr(i,k) = min(nsacr(i,k),ncr(i,k,3)/dtcld)
1216 endif
1217 !-------------------------------------------------------------
1218 ! pgacr: Accretion of rain by graupel [HL A12] [LFO 42]
1219 ! (T<T0: QR->QG) (T>=T0: enhance melting of graupel)
1220 !-------------------------------------------------------------
1221 if(qrs(i,k,3).gt.qcrmin .and. qrs(i,k,1).gt.qcrmin) then
1222 acrfac = 30.*rslope3(i,k,1)*rslope2(i,k,1)*rslope(i,k,3) &
1223 + 5.*rslope2(i,k,1)*rslope2(i,k,1)*rslope2(i,k,3) &
1224 + 2.*rslope3(i,k,1)*rslope3(i,k,3)
1225 pgacr(i,k) = pi*pi*ncr(i,k,3)*n0g*abs(vt2ave-vt2r)*(denr/den(i,k)) &
1226 *acrfac
1227 pgacr(i,k) = min(pgacr(i,k),qrs(i,k,1)/dtcld)
1228 endif
1229 !-------------------------------------------------------------
1230 ! ngacr: Accretion of rain by graupel [LH A24]
1231 ! (T<T0: NR->)
1232 !-------------------------------------------------------------
1233 if(qrs(i,k,3).gt.qcrmin .and. ncr(i,k,3).gt.nrmin) then
1234 acrfac = 1.5*rslope2(i,k,1)*rslope(i,k,3) &
1235 + 1.0*rslope(i,k,1)*rslope2(i,k,3) + .5*rslope3(i,k,3)
1236 ngacr(i,k) = pi*ncr(i,k,3)*n0g*abs(vt2ave-vt2r)*acrfac
1237 ngacr(i,k) = min(ngacr(i,k),ncr(i,k,3)/dtcld)
1238 endif
1239 !
1240 !-------------------------------------------------------------
1241 ! pgacs: Accretion of snow by graupel [HL A13] [LFO 29]
1242 ! (QS->QG) : This process is eliminated in V3.0 with the
1243 ! new combined snow/graupel fall speeds
1244 !-------------------------------------------------------------
1245 if(qrs(i,k,3).gt.qcrmin .and. qrs(i,k,2).gt.qcrmin) then
1246 pgacs(i,k) = 0.
1247 endif
1248 if(supcol.le.0) then
1249 xlf = xlf0
1250 !-------------------------------------------------------------
1251 ! pseml: Enhanced melting of snow by accretion of water [HL A34]
1252 ! (T>=T0: QS->QR)
1253 !-------------------------------------------------------------
1254 if(qrs(i,k,2).gt.0.) &
1255 pseml(i,k) = min(max(cliq*supcol*(paacw(i,k)+psacr(i,k)) &
1256 /xlf,-qrs(i,k,2)/dtcld),0.)
1257 !--------------------------------------------------------------
1258 ! nseml: Enhanced melting of snow by accretion of water [LH A29]
1259 ! (T>=T0: ->NR)
1260 !--------------------------------------------------------------
1261 if (qrs(i,k,2).gt.qcrmin) then
1262 sfac = rslope(i,k,2)*n0s*n0sfac(i,k)/qrs(i,k,2)
1263 nseml(i,k) = -sfac*pseml(i,k)
1264 endif
1265 !-------------------------------------------------------------
1266 ! pgeml: Enhanced melting of graupel by accretion of water [HL A24] [RH84 A21-A22]
1267 ! (T>=T0: QG->QR)
1268 !-------------------------------------------------------------
1269 if(qrs(i,k,3).gt.0.) &
1270 pgeml(i,k) = min(max(cliq*supcol*(paacw(i,k)+pgacr(i,k))/xlf &
1271 ,-qrs(i,k,3)/dtcld),0.)
1272 !--------------------------------------------------------------
1273 ! ngeml: Enhanced melting of graupel by accretion of water [LH A30]
1274 ! (T>=T0: -> NR)
1275 !--------------------------------------------------------------
1276 if (qrs(i,k,3).gt.qcrmin) then
1277 gfac = rslope(i,k,3)*n0g/qrs(i,k,3)
1278 ngeml(i,k) = -gfac*pgeml(i,k)
1279 endif
1280 endif
1281 if(supcol.gt.0) then
1282 !-------------------------------------------------------------
1283 ! pidep: Deposition/Sublimation rate of ice [HDC 9]
1284 ! (T<T0: QV->QI or QI->QV)
1285 !-------------------------------------------------------------
1286 if(qci(i,k,2).gt.0. .and. ifsat.ne.1) then
1287 pidep(i,k) = 4.*diameter*xni(i,k)*(rh(i,k,2)-1.)/work1(i,k,2)
1288 supice = satdt-prevp(i,k)
1289 if(pidep(i,k).lt.0.) then
1290 pidep(i,k) = max(max(pidep(i,k),satdt/2),supice)
1291 pidep(i,k) = max(pidep(i,k),-qci(i,k,2)/dtcld)
1292 else
1293 pidep(i,k) = min(min(pidep(i,k),satdt/2),supice)
1294 endif
1295 if(abs(prevp(i,k)+pidep(i,k)).ge.abs(satdt)) ifsat = 1
1296 endif
1297 !-------------------------------------------------------------
1298 ! psdep: deposition/sublimation rate of snow [HDC 14]
1299 ! (T<T0: QV->QS or QS->QV)
1300 !-------------------------------------------------------------
1301 if(qrs(i,k,2).gt.0. .and. ifsat.ne.1) then
1302 coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
1303 psdep(i,k) = (rh(i,k,2)-1.)*n0sfac(i,k)*(precs1*rslope2(i,k,2) &
1304 + precs2*work2(i,k)*coeres)/work1(i,k,2)
1305 supice = satdt-prevp(i,k)-pidep(i,k)
1306 if(psdep(i,k).lt.0.) then
1307 psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)/dtcld)
1308 psdep(i,k) = max(max(psdep(i,k),satdt/2),supice)
1309 else
1310 psdep(i,k) = min(min(psdep(i,k),satdt/2),supice)
1311 endif
1312 if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)).ge.abs(satdt)) ifsat = 1
1313 endif
1314 !-------------------------------------------------------------
1315 ! pgdep: deposition/sublimation rate of graupel [HL A21] [LFO 46]
1316 ! (T<T0: QV->QG or QG->QV)
1317 !-------------------------------------------------------------
1318 if(qrs(i,k,3).gt.0. .and. ifsat.ne.1) then
1319 coeres = rslope2(i,k,3)*sqrt(rslope(i,k,3)*rslopeb(i,k,3))
1320 pgdep(i,k) = (rh(i,k,2)-1.)*(precg1*rslope2(i,k,3) &
1321 + precg2*work2(i,k)*coeres)/work1(i,k,2)
1322 supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)
1323 if(pgdep(i,k).lt.0.) then
1324 pgdep(i,k) = max(pgdep(i,k),-qrs(i,k,3)/dtcld)
1325 pgdep(i,k) = max(max(pgdep(i,k),satdt/2),supice)
1326 else
1327 pgdep(i,k) = min(min(pgdep(i,k),satdt/2),supice)
1328 endif
1329 if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)+pgdep(i,k)).ge. &
1330 abs(satdt)) ifsat = 1
1331 endif
1332 !-------------------------------------------------------------
1333 ! pigen: generation(nucleation) of ice from vapor [HL 50] [HDC 7-8]
1334 ! (T<T0: QV->QI)
1335 !-------------------------------------------------------------
1336 if(supsat.gt.0. .and. ifsat.ne.1) then
1337 supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)-pgdep(i,k)
1338 xni0 = 1.e3*exp(0.1*supcol)
1339 roqi0 = 4.92e-11*xni0**1.33
1340 pigen(i,k) = max(0.,(roqi0/den(i,k)-max(qci(i,k,2),0.))/dtcld)
1341 pigen(i,k) = min(min(pigen(i,k),satdt),supice)
1342 endif
1343 !
1344 !-------------------------------------------------------------
1345 ! psaut: conversion(aggregation) of ice to snow [HDC 12]
1346 ! (T<T0: QI->QS)
1347 !-------------------------------------------------------------
1348 if(qci(i,k,2).gt.0.) then
1349 qimax = roqimax/den(i,k)
1350 psaut(i,k) = max(0.,(qci(i,k,2)-qimax)/dtcld)
1351 endif
1352 !
1353 !-------------------------------------------------------------
1354 ! pgaut: conversion(aggregation) of snow to graupel [HL A4] [LFO 37]
1355 ! (T<T0: QS->QG)
1356 !-------------------------------------------------------------
1357 if(qrs(i,k,2).gt.0.) then
1358 alpha2 = 1.e-3*exp(0.09*(-supcol))
1359 pgaut(i,k) = min(max(0.,alpha2*(qrs(i,k,2)-qs0)),qrs(i,k,2)/dtcld)
1360 endif
1361 endif
1362 !
1363 !-------------------------------------------------------------
1364 ! psevp: Evaporation of melting snow [HL A35] [RH83 A27]
1365 ! (T>=T0: QS->QV)
1366 !-------------------------------------------------------------
1367 if(supcol.lt.0.) then
1368 if(qrs(i,k,2).gt.0. .and. rh(i,k,1).lt.1.) then
1369 coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
1370 psevp(i,k) = (rh(i,k,1)-1.)*n0sfac(i,k)*(precs1*rslope2(i,k,2) &
1371 +precs2*work2(i,k)*coeres)/work1(i,k,1)
1372 psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)/dtcld),0.)
1373 endif
1374 !-------------------------------------------------------------
1375 ! pgevp: Evaporation of melting graupel [HL A25] [RH84 A19]
1376 ! (T>=T0: QG->QV)
1377 !-------------------------------------------------------------
1378 if(qrs(i,k,3).gt.0. .and. rh(i,k,1).lt.1.) then
1379 coeres = rslope2(i,k,3)*sqrt(rslope(i,k,3)*rslopeb(i,k,3))
1380 pgevp(i,k) = (rh(i,k,1)-1.)*(precg1*rslope2(i,k,3) &
1381 + precg2*work2(i,k)*coeres)/work1(i,k,1)
1382 pgevp(i,k) = min(max(pgevp(i,k),-qrs(i,k,3)/dtcld),0.)
1383 endif
1384 endif
1385 enddo
1386 enddo
1387 !
1388 !
1389 !----------------------------------------------------------------
1390 ! check mass conservation of generation terms and feedback to the
1392 !
1393 do k = kts, kte
1394 do i = its, ite
1395 !
1396 delta2=0.
1397 delta3=0.
1398 if(qrs(i,k,1).lt.1.e-4 .and. qrs(i,k,2).lt.1.e-4) delta2=1.
1399 if(qrs(i,k,1).lt.1.e-4) delta3=1.
1400 if(t(i,k).le.t0c) then
1401 !
1402 ! cloud water
1403 !
1404 value = max(qmin,qci(i,k,1))
1405 source = (praut(i,k)+pracw(i,k)+paacw(i,k)+paacw(i,k)-prevp_s(i,k))&
1406 *dtcld
1407 if (source.gt.value) then
1408 factor = value/source
1409 praut(i,k) = praut(i,k)*factor
1410 pracw(i,k) = pracw(i,k)*factor
1411 paacw(i,k) = paacw(i,k)*factor
1412 prevp_s(i,k) = prevp_s(i,k)*factor
1413 endif
1414 !
1415 ! cloud ice
1416 !
1417 value = max(qmin,qci(i,k,2))
1418 source = (psaut(i,k)-pigen(i,k)-pidep(i,k)+praci(i,k)+psaci(i,k) &
1419 +pgaci(i,k))*dtcld
1420 if (source.gt.value) then
1421 factor = value/source
1422 psaut(i,k) = psaut(i,k)*factor
1423 pigen(i,k) = pigen(i,k)*factor
1424 pidep(i,k) = pidep(i,k)*factor
1425 praci(i,k) = praci(i,k)*factor
1426 psaci(i,k) = psaci(i,k)*factor
1427 pgaci(i,k) = pgaci(i,k)*factor
1428 endif
1429 !
1430 ! rain
1431 !
1432 value = max(qmin,qrs(i,k,1))
1433 source = (-praut(i,k)-prevp(i,k)-pracw(i,k)+piacr(i,k) &
1434 +prevp_s(i,k)+psacr(i,k)+pgacr(i,k))*dtcld
1435 if (source.gt.value) then
1436 factor = value/source
1437 praut(i,k) = praut(i,k)*factor
1438 prevp(i,k) = prevp(i,k)*factor
1439 pracw(i,k) = pracw(i,k)*factor
1440 piacr(i,k) = piacr(i,k)*factor
1441 psacr(i,k) = psacr(i,k)*factor
1442 pgacr(i,k) = pgacr(i,k)*factor
1443 prevp_s(i,k) = prevp_s(i,k)*factor
1444 endif
1445 !
1446 ! snow
1447 !
1448 value = max(qmin,qrs(i,k,2))
1449 source = -(psdep(i,k)+psaut(i,k)-pgaut(i,k)+paacw(i,k) &
1450 +piacr(i,k)*delta3+praci(i,k)*delta3 &
1451 -pracs(i,k)*(1.-delta2)+psacr(i,k)*delta2 &
1452 +psaci(i,k)-pgacs(i,k) )*dtcld
1453 if (source.gt.value) then
1454 factor = value/source
1455 psdep(i,k) = psdep(i,k)*factor
1456 psaut(i,k) = psaut(i,k)*factor
1457 pgaut(i,k) = pgaut(i,k)*factor
1458 paacw(i,k) = paacw(i,k)*factor
1459 piacr(i,k) = piacr(i,k)*factor
1460 praci(i,k) = praci(i,k)*factor
1461 psaci(i,k) = psaci(i,k)*factor
1462 pracs(i,k) = pracs(i,k)*factor
1463 psacr(i,k) = psacr(i,k)*factor
1464 pgacs(i,k) = pgacs(i,k)*factor
1465 endif
1466 !
1467 ! graupel
1468 !
1469 value = max(qmin,qrs(i,k,3))
1470 source = -(pgdep(i,k)+pgaut(i,k) &
1471 +piacr(i,k)*(1.-delta3)+praci(i,k)*(1.-delta3) &
1472 +psacr(i,k)*(1.-delta2)+pracs(i,k)*(1.-delta2) &
1473 +pgaci(i,k)+paacw(i,k)+pgacr(i,k)+pgacs(i,k))*dtcld
1474 if (source.gt.value) then
1475 factor = value/source
1476 pgdep(i,k) = pgdep(i,k)*factor
1477 pgaut(i,k) = pgaut(i,k)*factor
1478 piacr(i,k) = piacr(i,k)*factor
1479 praci(i,k) = praci(i,k)*factor
1480 psacr(i,k) = psacr(i,k)*factor
1481 pracs(i,k) = pracs(i,k)*factor
1482 paacw(i,k) = paacw(i,k)*factor
1483 pgaci(i,k) = pgaci(i,k)*factor
1484 pgacr(i,k) = pgacr(i,k)*factor
1485 pgacs(i,k) = pgacs(i,k)*factor
1486 endif
1487 !
1488 ! cloud
1489 !
1490 value = max(ncmin,ncr(i,k,2))
1491 source = (nraut(i,k)+nccol(i,k)+nracw(i,k) &
1492 +naacw(i,k)+naacw(i,k)-nrevp(i,k))*dtcld
1493 if (source.gt.value) then
1494 factor = value/source
1495 nraut(i,k) = nraut(i,k)*factor
1496 nccol(i,k) = nccol(i,k)*factor
1497 nracw(i,k) = nracw(i,k)*factor
1498 naacw(i,k) = naacw(i,k)*factor
1499 nrevp(i,k) = nrevp(i,k)*factor
1500 endif
1501 !
1502 ! rain
1503 !
1504 value = max(nrmin,ncr(i,k,3))
1505 source = (-nraut(i,k)+nrcol(i,k)+niacr(i,k)+nsacr(i,k)+ngacr(i,k) &
1506 +nrevp(i,k))*dtcld
1507 if (source.gt.value) then
1508 factor = value/source
1509 nraut(i,k) = nraut(i,k)*factor
1510 nrcol(i,k) = nrcol(i,k)*factor
1511 niacr(i,k) = niacr(i,k)*factor
1512 nsacr(i,k) = nsacr(i,k)*factor
1513 ngacr(i,k) = ngacr(i,k)*factor
1514 nrevp(i,k) = nrevp(i,k)*factor
1515 endif
1516 !
1517 work2(i,k)=-(prevp(i,k)+psdep(i,k)+pgdep(i,k)+pigen(i,k)+pidep(i,k))
1518 ! update
1519 q(i,k) = q(i,k)+work2(i,k)*dtcld
1520 qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
1521 +paacw(i,k)+paacw(i,k)+prevp_s(i,k))*dtcld,0.)
1522 qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
1523 +prevp(i,k)-piacr(i,k)-pgacr(i,k) &
1524 -psacr(i,k)-prevp_s(i,k))*dtcld,0.)
1525 qci(i,k,2) = max(qci(i,k,2)-(psaut(i,k)+praci(i,k) &
1526 +psaci(i,k)+pgaci(i,k)-pigen(i,k)-pidep(i,k)) &
1527 *dtcld,0.)
1528 qrs(i,k,2) = max(qrs(i,k,2)+(psdep(i,k)+psaut(i,k)+paacw(i,k) &
1529 -pgaut(i,k)+piacr(i,k)*delta3 &
1530 +praci(i,k)*delta3+psaci(i,k)-pgacs(i,k) &
1531 -pracs(i,k)*(1.-delta2)+psacr(i,k)*delta2) &
1532 *dtcld,0.)
1533 qrs(i,k,3) = max(qrs(i,k,3)+(pgdep(i,k)+pgaut(i,k) &
1534 +piacr(i,k)*(1.-delta3) &
1535 +praci(i,k)*(1.-delta3)+psacr(i,k)*(1.-delta2) &
1536 +pracs(i,k)*(1.-delta2)+pgaci(i,k)+paacw(i,k) &
1537 +pgacr(i,k)+pgacs(i,k))*dtcld,0.)
1538 ncr(i,k,2) = max(ncr(i,k,2)+(-nraut(i,k)-nccol(i,k)-nracw(i,k) &
1539 -naacw(i,k)-naacw(i,k)+nrevp(i,k))*dtcld,0.)
1540 ncr(i,k,3) = max(ncr(i,k,3)+(nraut(i,k)-nrcol(i,k)-niacr(i,k) &
1541 -nsacr(i,k)-ngacr(i,k)-nrevp(i,k))*dtcld,0.)
1542 xlf = xls-xl(i,k)
1543 xlwork2 = -xls*(psdep(i,k)+pgdep(i,k)+pidep(i,k)+pigen(i,k)) &
1544 -xl(i,k)*prevp(i,k)-xlf*(piacr(i,k)+paacw(i,k) &
1545 +paacw(i,k)+pgacr(i,k)+psacr(i,k))
1546 t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
1547 else
1548 !
1549 ! cloud water
1550 !
1551 value = max(qmin,qci(i,k,1))
1552 source= (praut(i,k)+pracw(i,k)+paacw(i,k)+paacw(i,k)-prevp_s(i,k)) &
1553 *dtcld
1554 if (source.gt.value) then
1555 factor = value/source
1556 praut(i,k) = praut(i,k)*factor
1557 pracw(i,k) = pracw(i,k)*factor
1558 paacw(i,k) = paacw(i,k)*factor
1559 prevp_s(i,k) = prevp_s(i,k)*factor
1560 endif
1561 !
1562 ! rain
1563 !
1564 value = max(qmin,qrs(i,k,1))
1565 source = (-paacw(i,k)-praut(i,k)+pseml(i,k)+pgeml(i,k) &
1566 +prevp_s(i,k)-pracw(i,k)-paacw(i,k)-prevp(i,k))*dtcld
1567 if (source.gt.value) then
1568 factor = value/source
1569 praut(i,k) = praut(i,k)*factor
1570 prevp(i,k) = prevp(i,k)*factor
1571 pracw(i,k) = pracw(i,k)*factor
1572 paacw(i,k) = paacw(i,k)*factor
1573 pseml(i,k) = pseml(i,k)*factor
1574 pgeml(i,k) = pgeml(i,k)*factor
1575 prevp_s(i,k) = prevp_s(i,k)*factor
1576 endif
1577 !
1578 ! snow
1579 !
1580 value = max(qcrmin,qrs(i,k,2))
1581 source=(pgacs(i,k)-pseml(i,k)-psevp(i,k))*dtcld
1582 if (source.gt.value) then
1583 factor = value/source
1584 pgacs(i,k) = pgacs(i,k)*factor
1585 psevp(i,k) = psevp(i,k)*factor
1586 pseml(i,k) = pseml(i,k)*factor
1587 endif
1588 !
1589 ! graupel
1590 !
1591 value = max(qcrmin,qrs(i,k,3))
1592 source=-(pgacs(i,k)+pgevp(i,k)+pgeml(i,k))*dtcld
1593 if (source.gt.value) then
1594 factor = value/source
1595 pgacs(i,k) = pgacs(i,k)*factor
1596 pgevp(i,k) = pgevp(i,k)*factor
1597 pgeml(i,k) = pgeml(i,k)*factor
1598 endif
1599 !
1600 ! cloud
1601 !
1602 value = max(ncmin,ncr(i,k,2))
1603 source = (+nraut(i,k)+nccol(i,k)+nracw(i,k)+naacw(i,k) &
1604 +naacw(i,k)-nrevp(i,k))*dtcld
1605 if (source.gt.value) then
1606 factor = value/source
1607 nraut(i,k) = nraut(i,k)*factor
1608 nccol(i,k) = nccol(i,k)*factor
1609 nracw(i,k) = nracw(i,k)*factor
1610 naacw(i,k) = naacw(i,k)*factor
1611 nrevp(i,k) = nrevp(i,k)*factor
1612 endif
1613 !
1614 ! rain
1615 !
1616 value = max(nrmin,ncr(i,k,3))
1617 source = (-nraut(i,k)+nrcol(i,k)-nseml(i,k)-ngeml(i,k) &
1618 +nrevp(i,k))*dtcld
1619 if (source.gt.value) then
1620 factor = value/source
1621 nraut(i,k) = nraut(i,k)*factor
1622 nrcol(i,k) = nrcol(i,k)*factor
1623 nrevp(i,k) = nrevp(i,k)*factor
1624 nseml(i,k) = nseml(i,k)*factor
1625 ngeml(i,k) = ngeml(i,k)*factor
1626 endif
1627 !
1628 work2(i,k)=-(prevp(i,k)+psevp(i,k)+pgevp(i,k))
1629 ! update
1630 q(i,k) = q(i,k)+work2(i,k)*dtcld
1631 qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
1632 +prevp_s(i,k)+paacw(i,k)+paacw(i,k))*dtcld,0.)
1633 qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
1634 +prevp(i,k)-prevp_s(i,k)+paacw(i,k)+paacw(i,k)-pseml(i,k) &
1635 -pgeml(i,k))*dtcld,0.)
1636 qrs(i,k,2) = max(qrs(i,k,2)+(psevp(i,k)-pgacs(i,k) &
1637 +pseml(i,k))*dtcld,0.)
1638 qrs(i,k,3) = max(qrs(i,k,3)+(pgacs(i,k)+pgevp(i,k) &
1639 +pgeml(i,k))*dtcld,0.)
1640 ncr(i,k,2) = max(ncr(i,k,2)+(-nraut(i,k)-nccol(i,k)-nracw(i,k) &
1641 -naacw(i,k)-naacw(i,k)+nrevp(i,k))*dtcld,0.)
1642 ncr(i,k,3) = max(ncr(i,k,3)+(nraut(i,k)-nrcol(i,k)+nseml(i,k) &
1643 +ngeml(i,k)-nrevp(i,k))*dtcld,0.)
1644 xlf = xls-xl(i,k)
1645 xlwork2 = -xl(i,k)*(prevp(i,k)+psevp(i,k)+pgevp(i,k)) &
1646 -xlf*(pseml(i,k)+pgeml(i,k))
1647 t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
1648 endif
1649 enddo
1650 enddo
1651 !
1652 ! Inline expansion for fpvs
1653 ! qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1654 ! qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1655 hsub = xls
1656 hvap = xlv0
1657 cvap = cpv
1658 ttp=t0c+0.01
1659 dldt=cvap-cliq
1660 xa=-dldt/rv
1661 xb=xa+hvap/(rv*ttp)
1662 dldti=cvap-cice
1663 xai=-dldti/rv
1664 xbi=xai+hsub/(rv*ttp)
1665 do k = kts, kte
1666 do i = its, ite
1667 tr=ttp/t(i,k)
1668 qs(i,k,1)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
1669 qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
1670 qs(i,k,1) = max(qs(i,k,1),qmin)
1671 tr=ttp/t(i,k)
1672 if(t(i,k).lt.ttp) then
1673 qs(i,k,2)=psat*exp(log(tr)*(xai))*exp(xbi*(1.-tr))
1674 else
1675 qs(i,k,2)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
1676 endif
1677 qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
1678 qs(i,k,2) = max(qs(i,k,2),qmin)
1679 rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
1680 enddo
1681 enddo
1682 !
1683 do k = kts, kte
1684 do i = its, ite
1685 !---------------------------------------------------------------
1686 ! rate of change of cloud drop concentration due to CCN activation
1687 ! pcact: QV -> QC [LH 8] [KK 14]
1688 ! ncact: NCCN -> NC [LH A2] [KK 12]
1689 !---------------------------------------------------------------
1690 if(rh(i,k,1).gt.1.) then
1691 ncact(i,k) = max(0.,((ncr(i,k,1)+ncr(i,k,2)) &
1692 *min(1.,(rh(i,k,1)/satmax)**actk) - ncr(i,k,2)))/dtcld
1693 ncact(i,k) =min(ncact(i,k),max(ncr(i,k,1),0.)/dtcld)
1694 pcact(i,k) = min(4.*pi*denr*(actr*1.E-6)**3*ncact(i,k)/ &
1695 (3.*den(i,k)),max(q(i,k),0.)/dtcld)
1696 q(i,k) = max(q(i,k)-pcact(i,k)*dtcld,0.)
1697 qci(i,k,1) = max(qci(i,k,1)+pcact(i,k)*dtcld,0.)
1698 ncr(i,k,1) = max(ncr(i,k,1)-ncact(i,k)*dtcld,0.)
1699 ncr(i,k,2) = max(ncr(i,k,2)+ncact(i,k)*dtcld,0.)
1700 t(i,k) = t(i,k)+pcact(i,k)*xl(i,k)/cpm(i,k)*dtcld
1701 endif
1702 !---------------------------------------------------------------
1703 ! pcond:condensational/evaporational rate of cloud water [HL A46] [RH83 A6]
1704 ! if there exists additional water vapor condensated/if
1705 ! evaporation of cloud water is not enough to remove subsaturation
1706 ! (QV->QC or QC->QV)
1707 !---------------------------------------------------------------
1708 tr=ttp/t(i,k)
1709 qs(i,k,1)=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
1710 qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
1711 qs(i,k,1) = max(qs(i,k,1),qmin)
1712 work1(i,k,1) = conden(t(i,k),q(i,k),qs(i,k,1),xl(i,k),cpm(i,k))
1713 work2(i,k) = qci(i,k,1)+work1(i,k,1)
1714 pcond(i,k) = min(max(work1(i,k,1)/dtcld,0.),max(q(i,k),0.)/dtcld)
1715 if(qci(i,k,1).gt.0. .and. work1(i,k,1).lt.0.) &
1716 pcond(i,k) = max(work1(i,k,1),-qci(i,k,1))/dtcld
1717 !----------------------------------------------------------------
1718 ! ncevp: evpration of Cloud number concentration [LH A3]
1719 ! (NC->NCCN)
1720 !----------------------------------------------------------------
1721 if(pcond(i,k).eq.-qci(i,k,1)/dtcld) then
1722 ncr(i,k,2) = 0.
1723 ncr(i,k,1) = ncr(i,k,1)+ncr(i,k,2)
1724 endif
1725 !
1726 q(i,k) = q(i,k)-pcond(i,k)*dtcld
1727 qci(i,k,1) = max(qci(i,k,1)+pcond(i,k)*dtcld,0.)
1728 t(i,k) = t(i,k)+pcond(i,k)*xl(i,k)/cpm(i,k)*dtcld
1729 enddo
1730 enddo
1731 !
1732 !----------------------------------------------------------------
1733 ! padding for small values
1734 !
1735 do k = kts, kte
1736 do i = its, ite
1737 if(qci(i,k,1).le.qmin) qci(i,k,1) = 0.0
1738 if(qci(i,k,2).le.qmin) qci(i,k,2) = 0.0
1739 enddo
1740 enddo
1741 enddo ! big loops
1742 END SUBROUTINE wdm62d
1743 ! ...................................................................
1744 REAL FUNCTION rgmma(x)
1745 !-------------------------------------------------------------------
1746 IMPLICIT NONE
1747 !-------------------------------------------------------------------
1748 ! rgmma function: use infinite product form
1749 REAL :: euler
1750 PARAMETER (euler=0.577215664901532)
1751 REAL :: x, y
1752 INTEGER :: i
1753 if(x.eq.1.)then
1754 rgmma=0.
1755 else
1756 rgmma=x*exp(euler*x)
1757 do i=1,10000
1758 y=float(i)
1759 rgmma=rgmma*(1.000+x/y)*exp(-x/y)
1760 enddo
1761 rgmma=1./rgmma
1762 endif
1763 END FUNCTION rgmma
1764 !
1765 !--------------------------------------------------------------------------
1766 REAL FUNCTION fpvs(t,ice,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c)
1767 !--------------------------------------------------------------------------
1768 IMPLICIT NONE
1769 !--------------------------------------------------------------------------
1770 REAL t,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c,dldt,xa,xb,dldti, &
1771 xai,xbi,ttp,tr
1772 INTEGER ice
1773 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1774 ttp=t0c+0.01
1775 dldt=cvap-cliq
1776 xa=-dldt/rv
1777 xb=xa+hvap/(rv*ttp)
1778 dldti=cvap-cice
1779 xai=-dldti/rv
1780 xbi=xai+hsub/(rv*ttp)
1781 tr=ttp/t
1782 if(t.lt.ttp .and. ice.eq.1) then
1783 fpvs=psat*(tr**xai)*exp(xbi*(1.-tr))
1784 else
1785 fpvs=psat*(tr**xa)*exp(xb*(1.-tr))
1786 endif
1787 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1788 END FUNCTION fpvs
1789 !-------------------------------------------------------------------
1790 SUBROUTINE wdm6init(den0,denr,dens,cl,cpv, ccn0, allowed_to_read)
1791 !-------------------------------------------------------------------
1792 IMPLICIT NONE
1793 !-------------------------------------------------------------------
1794 !.... constants which may not be tunable
1795 REAL, INTENT(IN) :: den0,denr,dens,cl,cpv,ccn0
1796 LOGICAL, INTENT(IN) :: allowed_to_read
1797 !
1798 pi = 4.*atan(1.)
1799 xlv1 = cl-cpv
1800 !
1801 qc0 = 4./3.*pi*denr*r0**3*xncr/den0 ! 0.419e-3 -- .61e-3
1802 qck1 = .104*9.8*peaut/(xncr*denr)**(1./3.)/xmyu*den0**(4./3.) ! 7.03
1803 pidnc = pi*denr/6.
1804 !
1805 bvtr1 = 1.+bvtr
1806 bvtr2 = 2.+bvtr
1807 bvtr3 = 3.+bvtr
1808 bvtr4 = 4.+bvtr
1809 bvtr5 = 5.+bvtr
1810 bvtr6 = 6.+bvtr
1811 bvtr7 = 7.+bvtr
1812 bvtr2o5 = 2.5+.5*bvtr
1813 bvtr3o5 = 3.5+.5*bvtr
1814 g1pbr = rgmma(bvtr1)
1815 g2pbr = rgmma(bvtr2)
1816 g3pbr = rgmma(bvtr3)
1817 g4pbr = rgmma(bvtr4) ! 17.837825
1818 g5pbr = rgmma(bvtr5)
1819 g6pbr = rgmma(bvtr6)
1820 g7pbr = rgmma(bvtr7)
1821 g5pbro2 = rgmma(bvtr2o5)
1822 g7pbro2 = rgmma(bvtr3o5)
1823 pvtr = avtr*g5pbr/24.
1824 pvtrn = avtr*g2pbr
1825 eacrr = 1.0
1826 pacrr = pi*n0r*avtr*g3pbr*.25*eacrr
1827 precr1 = 2.*pi*1.56
1828 precr2 = 2.*pi*.31*avtr**.5*g7pbro2
1829 pidn0r = pi*denr*n0r
1830 pidnr = 4.*pi*denr
1831 !
1832 xmmax = (dimax/dicon)**2
1833 roqimax = 2.08e22*dimax**8
1834 !
1835 bvts1 = 1.+bvts
1836 bvts2 = 2.5+.5*bvts
1837 bvts3 = 3.+bvts
1838 bvts4 = 4.+bvts
1839 g1pbs = rgmma(bvts1) !.8875
1840 g3pbs = rgmma(bvts3)
1841 g4pbs = rgmma(bvts4) ! 12.0786
1842 g5pbso2 = rgmma(bvts2)
1843 pvts = avts*g4pbs/6.
1844 pacrs = pi*n0s*avts*g3pbs*.25
1845 precs1 = 4.*n0s*.65
1846 precs2 = 4.*n0s*.44*avts**.5*g5pbso2
1847 pidn0s = pi*dens*n0s
1848 !
1849 pacrc = pi*n0s*avts*g3pbs*.25*eacrc
1850 !
1851 bvtg1 = 1.+bvtg
1852 bvtg2 = 2.5+.5*bvtg
1853 bvtg3 = 3.+bvtg
1854 bvtg4 = 4.+bvtg
1855 g1pbg = rgmma(bvtg1)
1856 g3pbg = rgmma(bvtg3)
1857 g4pbg = rgmma(bvtg4)
1858 g5pbgo2 = rgmma(bvtg2)
1859 pacrg = pi*n0g*avtg*g3pbg*.25
1860 pvtg = avtg*g4pbg/6.
1861 precg1 = 2.*pi*n0g*.78
1862 precg2 = 2.*pi*n0g*.31*avtg**.5*g5pbgo2
1863 pidn0g = pi*deng*n0g
1864 !
1865 rslopecmax = 1./lamdacmax
1866 rslopermax = 1./lamdarmax
1867 rslopesmax = 1./lamdasmax
1868 rslopegmax = 1./lamdagmax
1869 rsloperbmax = rslopermax ** bvtr
1870 rslopesbmax = rslopesmax ** bvts
1871 rslopegbmax = rslopegmax ** bvtg
1872 rslopec2max = rslopecmax * rslopecmax
1873 rsloper2max = rslopermax * rslopermax
1874 rslopes2max = rslopesmax * rslopesmax
1875 rslopeg2max = rslopegmax * rslopegmax
1876 rslopec3max = rslopec2max * rslopecmax
1877 rsloper3max = rsloper2max * rslopermax
1878 rslopes3max = rslopes2max * rslopesmax
1879 rslopeg3max = rslopeg2max * rslopegmax
1880 !
1881 END SUBROUTINE wdm6init
1882 !------------------------------------------------------------------------------
1883 subroutine slope_wdm6(qrs,ncr,den,denfac,t,rslope,rslopeb,rslope2,rslope3, &
1884 vt,vtn,its,ite,kts,kte)
1885 IMPLICIT NONE
1886 INTEGER :: its,ite, jts,jte, kts,kte
1887 REAL, DIMENSION( its:ite , kts:kte,3) :: &
1888 qrs, &
1889 rslope, &
1890 rslopeb, &
1891 rslope2, &
1892 rslope3, &
1893 vt
1894 REAL, DIMENSION( its:ite , kts:kte) :: &
1895 ncr, &
1896 vtn, &
1897 den, &
1898 denfac, &
1899 t
1900 REAL, PARAMETER :: t0c = 273.15
1901 REAL, DIMENSION( its:ite , kts:kte ) :: &
1902 n0sfac
1903 REAL :: lamdar, lamdas, lamdag, x, y, z, supcol
1904 integer :: i, j, k
1905 !----------------------------------------------------------------
1906 ! size distributions: (x=mixing ratio, y=air density):
1907 ! valid for mixing ratio > 1.e-9 kg/kg.
1908 !
1909 ! Optimizatin : A**B => exp(log(A)*(B))
1910 lamdar(x,y,z)= exp(log(((pidnr*z)/(x*y)))*((.33333333)))
1911 lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y))) ! (pidn0s*z/(x*y))**.25
1912 lamdag(x,y)= sqrt(sqrt(pidn0g/(x*y))) ! (pidn0g/(x*y))**.25
1913 !
1914 do k = kts, kte
1915 do i = its, ite
1916 supcol = t0c-t(i,k)
1917 !---------------------------------------------------------------
1918 ! n0s: Intercept parameter for snow [m-4] [HDC 6]
1919 !---------------------------------------------------------------
1920 n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
1921 if(qrs(i,k,1).le.qcrmin .or. ncr(i,k).le.nrmin ) then
1922 rslope(i,k,1) = rslopermax
1923 rslopeb(i,k,1) = rsloperbmax
1924 rslope2(i,k,1) = rsloper2max
1925 rslope3(i,k,1) = rsloper3max
1926 else
1927 rslope(i,k,1) = min(1./lamdar(qrs(i,k,1),den(i,k),ncr(i,k)),1.e-3)
1928 rslopeb(i,k,1) = rslope(i,k,1)**bvtr
1929 rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
1930 rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
1931 endif
1932 if(qrs(i,k,2).le.qcrmin) then
1933 rslope(i,k,2) = rslopesmax
1934 rslopeb(i,k,2) = rslopesbmax
1935 rslope2(i,k,2) = rslopes2max
1936 rslope3(i,k,2) = rslopes3max
1937 else
1938 rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k))
1939 rslopeb(i,k,2) = rslope(i,k,2)**bvts
1940 rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
1941 rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
1942 endif
1943 if(qrs(i,k,3).le.qcrmin) then
1944 rslope(i,k,3) = rslopegmax
1945 rslopeb(i,k,3) = rslopegbmax
1946 rslope2(i,k,3) = rslopeg2max
1947 rslope3(i,k,3) = rslopeg3max
1948 else
1949 rslope(i,k,3) = 1./lamdag(qrs(i,k,3),den(i,k))
1950 rslopeb(i,k,3) = rslope(i,k,3)**bvtg
1951 rslope2(i,k,3) = rslope(i,k,3)*rslope(i,k,3)
1952 rslope3(i,k,3) = rslope2(i,k,3)*rslope(i,k,3)
1953 endif
1954 vt(i,k,1) = pvtr*rslopeb(i,k,1)*denfac(i,k)
1955 vt(i,k,2) = pvts*rslopeb(i,k,2)*denfac(i,k)
1956 vt(i,k,3) = pvtg*rslopeb(i,k,3)*denfac(i,k)
1957 vtn(i,k) = pvtrn*rslopeb(i,k,1)*denfac(i,k)
1958 if(qrs(i,k,1).le.0.0) vt(i,k,1) = 0.0
1959 if(qrs(i,k,2).le.0.0) vt(i,k,2) = 0.0
1960 if(qrs(i,k,3).le.0.0) vt(i,k,3) = 0.0
1961 if(ncr(i,k).le.0.0) vtn(i,k) = 0.0
1962 enddo
1963 enddo
1964 END subroutine slope_wdm6
1965 !-----------------------------------------------------------------------------
1966 subroutine slope_rain(qrs,ncr,den,denfac,t,rslope,rslopeb,rslope2,rslope3, &
1967 vt,vtn,its,ite,kts,kte)
1968 IMPLICIT NONE
1969 INTEGER :: its,ite, jts,jte, kts,kte
1970 REAL, DIMENSION( its:ite , kts:kte) :: &
1971 qrs, &
1972 ncr, &
1973 rslope, &
1974 rslopeb, &
1975 rslope2, &
1976 rslope3, &
1977 vt, &
1978 vtn, &
1979 den, &
1980 denfac, &
1981 t
1982 REAL, PARAMETER :: t0c = 273.15
1983 REAL, DIMENSION( its:ite , kts:kte ) :: &
1984 n0sfac
1985 REAL :: lamdar, x, y, z, supcol
1986 integer :: i, j, k
1987 !----------------------------------------------------------------
1988 ! size distributions: (x=mixing ratio, y=air density):
1989 ! valid for mixing ratio > 1.e-9 kg/kg.
1990 lamdar(x,y,z)= exp(log(((pidnr*z)/(x*y)))*((.33333333)))
1991 !
1992 do k = kts, kte
1993 do i = its, ite
1994 if(qrs(i,k).le.qcrmin .or. ncr(i,k).le.nrmin) then
1995 rslope(i,k) = rslopermax
1996 rslopeb(i,k) = rsloperbmax
1997 rslope2(i,k) = rsloper2max
1998 rslope3(i,k) = rsloper3max
1999 else
2000 rslope(i,k) = min(1./lamdar(qrs(i,k),den(i,k),ncr(i,k)),1.e-3)
2001 rslopeb(i,k) = rslope(i,k)**bvtr
2002 rslope2(i,k) = rslope(i,k)*rslope(i,k)
2003 rslope3(i,k) = rslope2(i,k)*rslope(i,k)
2004 endif
2005 vt(i,k) = pvtr*rslopeb(i,k)*denfac(i,k)
2006 vtn(i,k) = pvtrn*rslopeb(i,k)*denfac(i,k)
2007 if(qrs(i,k).le.0.0) vt(i,k) = 0.0
2008 if(ncr(i,k).le.0.0) vtn(i,k) = 0.0
2009 enddo
2010 enddo
2011 END subroutine slope_rain
2012 !------------------------------------------------------------------------------
2013 subroutine slope_snow(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3, &
2014 vt,its,ite,kts,kte)
2015 IMPLICIT NONE
2016 INTEGER :: its,ite, jts,jte, kts,kte
2017 REAL, DIMENSION( its:ite , kts:kte) :: &
2018 qrs, &
2019 rslope, &
2020 rslopeb, &
2021 rslope2, &
2022 rslope3, &
2023 vt, &
2024 den, &
2025 denfac, &
2026 t
2027 REAL, PARAMETER :: t0c = 273.15
2028 REAL, DIMENSION( its:ite , kts:kte ) :: &
2029 n0sfac
2030 REAL :: lamdas, x, y, z, supcol
2031 integer :: i, j, k
2032 !----------------------------------------------------------------
2033 ! size distributions: (x=mixing ratio, y=air density):
2034 ! valid for mixing ratio > 1.e-9 kg/kg.
2035 lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y))) ! (pidn0s*z/(x*y))**.25
2036 !
2037 do k = kts, kte
2038 do i = its, ite
2039 supcol = t0c-t(i,k)
2040 !---------------------------------------------------------------
2041 ! n0s: Intercept parameter for snow [m-4] [HDC 6]
2042 !---------------------------------------------------------------
2043 n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
2044 if(qrs(i,k).le.qcrmin)then
2045 rslope(i,k) = rslopesmax
2046 rslopeb(i,k) = rslopesbmax
2047 rslope2(i,k) = rslopes2max
2048 rslope3(i,k) = rslopes3max
2049 else
2050 rslope(i,k) = 1./lamdas(qrs(i,k),den(i,k),n0sfac(i,k))
2051 rslopeb(i,k) = rslope(i,k)**bvts
2052 rslope2(i,k) = rslope(i,k)*rslope(i,k)
2053 rslope3(i,k) = rslope2(i,k)*rslope(i,k)
2054 endif
2055 vt(i,k) = pvts*rslopeb(i,k)*denfac(i,k)
2056 if(qrs(i,k).le.0.0) vt(i,k) = 0.0
2057 enddo
2058 enddo
2059 END subroutine slope_snow
2060 !----------------------------------------------------------------------------------
2061 subroutine slope_graup(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3, &
2062 vt,its,ite,kts,kte)
2063 IMPLICIT NONE
2064 INTEGER :: its,ite, jts,jte, kts,kte
2065 REAL, DIMENSION( its:ite , kts:kte) :: &
2066 qrs, &
2067 rslope, &
2068 rslopeb, &
2069 rslope2, &
2070 rslope3, &
2071 vt, &
2072 den, &
2073 denfac, &
2074 t
2075 REAL, PARAMETER :: t0c = 273.15
2076 REAL, DIMENSION( its:ite , kts:kte ) :: &
2077 n0sfac
2078 REAL :: lamdag, x, y, z, supcol
2079 integer :: i, j, k
2080 !----------------------------------------------------------------
2081 ! size distributions: (x=mixing ratio, y=air density):
2082 ! valid for mixing ratio > 1.e-9 kg/kg.
2083 lamdag(x,y)= sqrt(sqrt(pidn0g/(x*y))) ! (pidn0g/(x*y))**.25
2084 !
2085 do k = kts, kte
2086 do i = its, ite
2087 !---------------------------------------------------------------
2088 ! n0s: Intercept parameter for snow [m-4] [HDC 6]
2089 !---------------------------------------------------------------
2090 if(qrs(i,k).le.qcrmin)then
2091 rslope(i,k) = rslopegmax
2092 rslopeb(i,k) = rslopegbmax
2093 rslope2(i,k) = rslopeg2max
2094 rslope3(i,k) = rslopeg3max
2095 else
2096 rslope(i,k) = 1./lamdag(qrs(i,k),den(i,k))
2097 rslopeb(i,k) = rslope(i,k)**bvtg
2098 rslope2(i,k) = rslope(i,k)*rslope(i,k)
2099 rslope3(i,k) = rslope2(i,k)*rslope(i,k)
2100 endif
2101 vt(i,k) = pvtg*rslopeb(i,k)*denfac(i,k)
2102 if(qrs(i,k).le.0.0) vt(i,k) = 0.0
2103 enddo
2104 enddo
2105 END subroutine slope_graup
2106 !---------------------------------------------------------------------------------
2107 !-------------------------------------------------------------------
2108 SUBROUTINE nislfv_rain_plmr(im,km,denl,denfacl,tkl,dzl,wwl,rql,rnl,precip,dt,id,iter,rid)
2109 !-------------------------------------------------------------------
2110 !
2111 ! for non-iteration semi-Lagrangain forward advection for cloud
2112 ! with mass conservation and positive definite advection
2113 ! 2nd order interpolation with monotonic piecewise linear method
2114 ! this routine is under assumption of decfl < 1 for semi_Lagrangian
2115 !
2116 ! dzl depth of model layer in meter
2117 ! wwl terminal velocity at model layer m/s
2118 ! rql cloud density*mixing ration
2119 ! precip precipitation
2120 ! dt time step
2121 ! id kind of precip: 0 test case; 1 raindrop
2122 ! iter how many time to guess mean terminal velocity: 0 pure forward.
2123 ! 0 : use departure wind for advection
2124 ! 1 : use mean wind for advection
2125 ! > 1 : use mean wind after iter-1 iterations
2126 ! rid : 1 for number 0 for mixing ratio
2127 !
2128 ! author: hann-ming henry juang <henry.juang@noaa.gov>
2129 ! implemented by song-you hong
2130 !
2131 implicit none
2132 integer im,km,id
2133 real dt
2134 real dzl(im,km),wwl(im,km),rql(im,km),rnl(im,km),precip(im)
2135 real denl(im,km),denfacl(im,km),tkl(im,km)
2136 !
2137 integer i,k,n,m,kk,kb,kt,iter,rid
2138 real tl,tl2,qql,dql,qqd
2139 real th,th2,qqh,dqh
2140 real zsum,qsum,dim,dip,c1,con1,fa1,fa2
2141 real allold, allnew, zz, dzamin, cflmax, decfl
2142 real dz(km), ww(km), qq(km), nr(km), wd(km), wa(km), wa2(km), was(km)
2143 real den(km), denfac(km), tk(km)
2144 real wi(km+1), zi(km+1), za(km+1)
2145 real qn(km), qr(km),tmp(km),tmp1(km),tmp2(km),tmp3(km)
2146 real dza(km+1), qa(km+1), qmi(km+1), qpi(km+1)
2147 !
2148 precip(:) = 0.0
2149 !
2150 i_loop : do i=1,im
2151 ! -----------------------------------
2152 dz(:) = dzl(i,:)
2153 qq(:) = rql(i,:)
2154 nr(:) = rnl(i,:)
2155 if(rid .eq. 1) nr(:) = rnl(i,:)/denl(i,:)
2156 ww(:) = wwl(i,:)
2157 den(:) = denl(i,:)
2158 denfac(:) = denfacl(i,:)
2159 tk(:) = tkl(i,:)
2160 ! skip for no precipitation for all layers
2161 allold = 0.0
2162 do k=1,km
2163 allold = allold + qq(k)
2164 enddo
2165 if(allold.le.0.0) then
2166 cycle i_loop
2167 endif
2168 !
2169 ! compute interface values
2170 zi(1)=0.0
2171 do k=1,km
2172 zi(k+1) = zi(k)+dz(k)
2173 enddo
2174 !
2175 ! save departure wind
2176 wd(:) = ww(:)
2177 n=1
2178 100 continue
2179 ! plm is 2nd order, we can use 2nd order wi or 3rd order wi
2180 ! 2nd order interpolation to get wi
2181 wi(1) = ww(1)
2182 wi(km+1) = ww(km)
2183 do k=2,km
2184 wi(k) = (ww(k)*dz(k-1)+ww(k-1)*dz(k))/(dz(k-1)+dz(k))
2185 enddo
2186 ! 3rd order interpolation to get wi
2187 fa1 = 9./16.
2188 fa2 = 1./16.
2189 wi(1) = ww(1)
2190 wi(2) = 0.5*(ww(2)+ww(1))
2191 do k=3,km-1
2192 wi(k) = fa1*(ww(k)+ww(k-1))-fa2*(ww(k+1)+ww(k-2))
2193 enddo
2194 wi(km) = 0.5*(ww(km)+ww(km-1))
2195 wi(km+1) = ww(km)
2196 !
2197 ! terminate of top of raingroup
2198 do k=2,km
2199 if( ww(k).eq.0.0 ) wi(k)=ww(k-1)
2200 enddo
2201 !
2202 ! diffusivity of wi
2203 con1 = 0.05
2204 do k=km,1,-1
2205 decfl = (wi(k+1)-wi(k))*dt/dz(k)
2206 if( decfl .gt. con1 ) then
2207 wi(k) = wi(k+1) - con1*dz(k)/dt
2208 endif
2209 enddo
2210 ! compute arrival point
2211 do k=1,km+1
2212 za(k) = zi(k) - wi(k)*dt
2213 enddo
2214 !
2215 do k=1,km
2216 dza(k) = za(k+1)-za(k)
2217 enddo
2218 dza(km+1) = zi(km+1) - za(km+1)
2219 !
2220 ! computer deformation at arrival point
2221 do k=1,km
2222 qa(k) = qq(k)*dz(k)/dza(k)
2223 qr(k) = qa(k)/den(k)
2224 if(rid .eq. 1) qr(k) = qa(K)
2225 enddo
2226 qa(km+1) = 0.0
2227 ! call maxmin(km,1,qa,' arrival points ')
2228 !
2229 ! compute arrival terminal velocity, and estimate mean terminal velocity
2230 ! then back to use mean terminal velocity
2231 if( n.le.iter ) then
2232 if(rid.eq.1) then
2233 call slope_rain(nr,qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,wa2,1,1,1,km)
2234 else
2235 call slope_rain(qr,nr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,wa2,1,1,1,km)
2236 endif
2237 if(rid.eq.1) wa(:) = wa2(:)
2238 if( n.ge.2 ) wa(1:km)=0.5*(wa(1:km)+was(1:km))
2239 do k=1,km
2240 !#ifdef DEBUG
2241 ! print*,' slope_wsm3 ',qr(k)*1000.,den(k),denfac(k),tk(k),tmp(k),tmp1(k),tmp2(k),ww(k),wa(k)
2242 !#endif
2243 ! mean wind is average of departure and new arrival winds
2244 ww(k) = 0.5* ( wd(k)+wa(k) )
2245 enddo
2246 was(:) = wa(:)
2247 n=n+1
2248 go to 100
2249 endif
2250 !
2251 ! estimate values at arrival cell interface with monotone
2252 do k=2,km
2253 dip=(qa(k+1)-qa(k))/(dza(k+1)+dza(k))
2254 dim=(qa(k)-qa(k-1))/(dza(k-1)+dza(k))
2255 if( dip*dim.le.0.0 ) then
2256 qmi(k)=qa(k)
2257 qpi(k)=qa(k)
2258 else
2259 qpi(k)=qa(k)+0.5*(dip+dim)*dza(k)
2260 qmi(k)=2.0*qa(k)-qpi(k)
2261 if( qpi(k).lt.0.0 .or. qmi(k).lt.0.0 ) then
2262 qpi(k) = qa(k)
2263 qmi(k) = qa(k)
2264 endif
2265 endif
2266 enddo
2267 qpi(1)=qa(1)
2268 qmi(1)=qa(1)
2269 qmi(km+1)=qa(km+1)
2270 qpi(km+1)=qa(km+1)
2271 !
2272 ! interpolation to regular point
2273 qn = 0.0
2274 kb=1
2275 kt=1
2276 intp : do k=1,km
2277 kb=max(kb-1,1)
2278 kt=max(kt-1,1)
2279 ! find kb and kt
2280 if( zi(k).ge.za(km+1) ) then
2281 exit intp
2282 else
2283 find_kb : do kk=kb,km
2284 if( zi(k).le.za(kk+1) ) then
2285 kb = kk
2286 exit find_kb
2287 else
2288 cycle find_kb
2289 endif
2290 enddo find_kb
2291 find_kt : do kk=kt,km
2292 if( zi(k+1).le.za(kk) ) then
2293 kt = kk
2294 exit find_kt
2295 else
2296 cycle find_kt
2297 endif
2298 enddo find_kt
2299 kt = kt - 1
2300 ! compute q with piecewise constant method
2301 if( kt.eq.kb ) then
2302 tl=(zi(k)-za(kb))/dza(kb)
2303 th=(zi(k+1)-za(kb))/dza(kb)
2304 tl2=tl*tl
2305 th2=th*th
2306 qqd=0.5*(qpi(kb)-qmi(kb))
2307 qqh=qqd*th2+qmi(kb)*th
2308 qql=qqd*tl2+qmi(kb)*tl
2309 qn(k) = (qqh-qql)/(th-tl)
2310 else if( kt.gt.kb ) then
2311 tl=(zi(k)-za(kb))/dza(kb)
2312 tl2=tl*tl
2313 qqd=0.5*(qpi(kb)-qmi(kb))
2314 qql=qqd*tl2+qmi(kb)*tl
2315 dql = qa(kb)-qql
2316 zsum = (1.-tl)*dza(kb)
2317 qsum = dql*dza(kb)
2318 if( kt-kb.gt.1 ) then
2319 do m=kb+1,kt-1
2320 zsum = zsum + dza(m)
2321 qsum = qsum + qa(m) * dza(m)
2322 enddo
2323 endif
2324 th=(zi(k+1)-za(kt))/dza(kt)
2325 th2=th*th
2326 qqd=0.5*(qpi(kt)-qmi(kt))
2327 dqh=qqd*th2+qmi(kt)*th
2328 zsum = zsum + th*dza(kt)
2329 qsum = qsum + dqh*dza(kt)
2330 qn(k) = qsum/zsum
2331 endif
2332 cycle intp
2333 endif
2334 !
2335 enddo intp
2336 !
2337 ! rain out
2338 sum_precip: do k=1,km
2339 if( za(k).lt.0.0 .and. za(k+1).lt.0.0 ) then
2340 precip(i) = precip(i) + qa(k)*dza(k)
2341 cycle sum_precip
2342 else if ( za(k).lt.0.0 .and. za(k+1).ge.0.0 ) then
2343 precip(i) = precip(i) + qa(k)*(0.0-za(k))
2344 exit sum_precip
2345 endif
2346 exit sum_precip
2347 enddo sum_precip
2348 !
2349 ! replace the new values
2350 rql(i,:) = qn(:)
2351 !
2352 ! ----------------------------------
2353 enddo i_loop
2354 !
2355 END SUBROUTINE nislfv_rain_plmr
2356 !-------------------------------------------------------------------
2357 SUBROUTINE nislfv_rain_plm6(im,km,denl,denfacl,tkl,dzl,wwl,rql,rql2, precip1, precip2,dt,id,iter)
2358 !-------------------------------------------------------------------
2359 !
2360 ! for non-iteration semi-Lagrangain forward advection for cloud
2361 ! with mass conservation and positive definite advection
2362 ! 2nd order interpolation with monotonic piecewise linear method
2363 ! this routine is under assumption of decfl < 1 for semi_Lagrangian
2364 !
2365 ! dzl depth of model layer in meter
2366 ! wwl terminal velocity at model layer m/s
2367 ! rql cloud density*mixing ration
2368 ! precip precipitation
2369 ! dt time step
2370 ! id kind of precip: 0 test case; 1 raindrop
2371 ! iter how many time to guess mean terminal velocity: 0 pure forward.
2372 ! 0 : use departure wind for advection
2373 ! 1 : use mean wind for advection
2374 ! > 1 : use mean wind after iter-1 iterations
2375 !
2376 ! author: hann-ming henry juang <henry.juang@noaa.gov>
2377 ! implemented by song-you hong
2378 !
2379 implicit none
2380 integer im,km,id
2381 real dt
2382 real dzl(im,km),wwl(im,km),rql(im,km),rql2(im,km),precip(im),precip1(im),precip2(im)
2383 real denl(im,km),denfacl(im,km),tkl(im,km)
2384 !
2385 integer i,k,n,m,kk,kb,kt,iter,ist
2386 real tl,tl2,qql,dql,qqd
2387 real th,th2,qqh,dqh
2388 real zsum,qsum,dim,dip,c1,con1,fa1,fa2
2389 real allold, allnew, zz, dzamin, cflmax, decfl
2390 real dz(km), ww(km), qq(km), qq2(km), wd(km), wa(km), wa2(km), was(km)
2391 real den(km), denfac(km), tk(km)
2392 real wi(km+1), zi(km+1), za(km+1)
2393 real qn(km), qr(km),qr2(km),tmp(km),tmp1(km),tmp2(km),tmp3(km)
2394 real dza(km+1), qa(km+1), qa2(km+1),qmi(km+1), qpi(km+1)
2395 !
2396 precip(:) = 0.0
2397 precip1(:) = 0.0
2398 precip2(:) = 0.0
2399 !
2400 i_loop : do i=1,im
2401 ! -----------------------------------
2402 dz(:) = dzl(i,:)
2403 qq(:) = rql(i,:)
2404 qq2(:) = rql2(i,:)
2405 ww(:) = wwl(i,:)
2406 den(:) = denl(i,:)
2407 denfac(:) = denfacl(i,:)
2408 tk(:) = tkl(i,:)
2409 ! skip for no precipitation for all layers
2410 allold = 0.0
2411 do k=1,km
2412 allold = allold + qq(k)
2413 enddo
2414 if(allold.le.0.0) then
2415 cycle i_loop
2416 endif
2417 !
2418 ! compute interface values
2419 zi(1)=0.0
2420 do k=1,km
2421 zi(k+1) = zi(k)+dz(k)
2422 enddo
2423 !
2424 ! save departure wind
2425 wd(:) = ww(:)
2426 n=1
2427 100 continue
2428 ! plm is 2nd order, we can use 2nd order wi or 3rd order wi
2429 ! 2nd order interpolation to get wi
2430 wi(1) = ww(1)
2431 wi(km+1) = ww(km)
2432 do k=2,km
2433 wi(k) = (ww(k)*dz(k-1)+ww(k-1)*dz(k))/(dz(k-1)+dz(k))
2434 enddo
2435 ! 3rd order interpolation to get wi
2436 fa1 = 9./16.
2437 fa2 = 1./16.
2438 wi(1) = ww(1)
2439 wi(2) = 0.5*(ww(2)+ww(1))
2440 do k=3,km-1
2441 wi(k) = fa1*(ww(k)+ww(k-1))-fa2*(ww(k+1)+ww(k-2))
2442 enddo
2443 wi(km) = 0.5*(ww(km)+ww(km-1))
2444 wi(km+1) = ww(km)
2445 !
2446 ! terminate of top of raingroup
2447 do k=2,km
2448 if( ww(k).eq.0.0 ) wi(k)=ww(k-1)
2449 enddo
2450 !
2451 ! diffusivity of wi
2452 con1 = 0.05
2453 do k=km,1,-1
2454 decfl = (wi(k+1)-wi(k))*dt/dz(k)
2455 if( decfl .gt. con1 ) then
2456 wi(k) = wi(k+1) - con1*dz(k)/dt
2457 endif
2458 enddo
2459 ! compute arrival point
2460 do k=1,km+1
2461 za(k) = zi(k) - wi(k)*dt
2462 enddo
2463 !
2464 do k=1,km
2465 dza(k) = za(k+1)-za(k)
2466 enddo
2467 dza(km+1) = zi(km+1) - za(km+1)
2468 !
2469 ! computer deformation at arrival point
2470 do k=1,km
2471 qa(k) = qq(k)*dz(k)/dza(k)
2472 qa2(k) = qq2(k)*dz(k)/dza(k)
2473 qr(k) = qa(k)/den(k)
2474 qr2(k) = qa2(k)/den(k)
2475 enddo
2476 qa(km+1) = 0.0
2477 qa2(km+1) = 0.0
2478 ! call maxmin(km,1,qa,' arrival points ')
2479 !
2480 ! compute arrival terminal velocity, and estimate mean terminal velocity
2481 ! then back to use mean terminal velocity
2482 if( n.le.iter ) then
2483 call slope_snow(qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,1,1,1,km)
2484 call slope_graup(qr2,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa2,1,1,1,km)
2485 do k = 1, km
2486 tmp(k) = max((qr(k)+qr2(k)), 1.E-15)
2487 IF ( tmp(k) .gt. 1.e-15 ) THEN
2488 wa(k) = (wa(k)*qr(k) + wa2(k)*qr2(k))/tmp(k)
2489 ELSE
2490 wa(k) = 0.
2491 ENDIF
2492 enddo
2493 if( n.ge.2 ) wa(1:km)=0.5*(wa(1:km)+was(1:km))
2494 do k=1,km
2495 !#ifdef DEBUG
2496 ! print*,' slope_wsm3 ',qr(k)*1000.,den(k),denfac(k),tk(k),tmp(k),tmp1(k),tmp2(k), &
2497 ! ww(k),wa(k)
2498 !#endif
2499 ! mean wind is average of departure and new arrival winds
2500 ww(k) = 0.5* ( wd(k)+wa(k) )
2501 enddo
2502 was(:) = wa(:)
2503 n=n+1
2504 go to 100
2505 endif
2506 ist_loop : do ist = 1, 2
2507 if (ist.eq.2) then
2508 qa(:) = qa2(:)
2509 endif
2510 !
2511 precip(i) = 0.
2512 !
2513 ! estimate values at arrival cell interface with monotone
2514 do k=2,km
2515 dip=(qa(k+1)-qa(k))/(dza(k+1)+dza(k))
2516 dim=(qa(k)-qa(k-1))/(dza(k-1)+dza(k))
2517 if( dip*dim.le.0.0 ) then
2518 qmi(k)=qa(k)
2519 qpi(k)=qa(k)
2520 else
2521 qpi(k)=qa(k)+0.5*(dip+dim)*dza(k)
2522 qmi(k)=2.0*qa(k)-qpi(k)
2523 if( qpi(k).lt.0.0 .or. qmi(k).lt.0.0 ) then
2524 qpi(k) = qa(k)
2525 qmi(k) = qa(k)
2526 endif
2527 endif
2528 enddo
2529 qpi(1)=qa(1)
2530 qmi(1)=qa(1)
2531 qmi(km+1)=qa(km+1)
2532 qpi(km+1)=qa(km+1)
2533 !
2534 ! interpolation to regular point
2535 qn = 0.0
2536 kb=1
2537 kt=1
2538 intp : do k=1,km
2539 kb=max(kb-1,1)
2540 kt=max(kt-1,1)
2541 ! find kb and kt
2542 if( zi(k).ge.za(km+1) ) then
2543 exit intp
2544 else
2545 find_kb : do kk=kb,km
2546 if( zi(k).le.za(kk+1) ) then
2547 kb = kk
2548 exit find_kb
2549 else
2550 cycle find_kb
2551 endif
2552 enddo find_kb
2553 find_kt : do kk=kt,km
2554 if( zi(k+1).le.za(kk) ) then
2555 kt = kk
2556 exit find_kt
2557 else
2558 cycle find_kt
2559 endif
2560 enddo find_kt
2561 kt = kt - 1
2562 ! compute q with piecewise constant method
2563 if( kt.eq.kb ) then
2564 tl=(zi(k)-za(kb))/dza(kb)
2565 th=(zi(k+1)-za(kb))/dza(kb)
2566 tl2=tl*tl
2567 th2=th*th
2568 qqd=0.5*(qpi(kb)-qmi(kb))
2569 qqh=qqd*th2+qmi(kb)*th
2570 qql=qqd*tl2+qmi(kb)*tl
2571 qn(k) = (qqh-qql)/(th-tl)
2572 else if( kt.gt.kb ) then
2573 tl=(zi(k)-za(kb))/dza(kb)
2574 tl2=tl*tl
2575 qqd=0.5*(qpi(kb)-qmi(kb))
2576 qql=qqd*tl2+qmi(kb)*tl
2577 dql = qa(kb)-qql
2578 zsum = (1.-tl)*dza(kb)
2579 qsum = dql*dza(kb)
2580 if( kt-kb.gt.1 ) then
2581 do m=kb+1,kt-1
2582 zsum = zsum + dza(m)
2583 qsum = qsum + qa(m) * dza(m)
2584 enddo
2585 endif
2586 th=(zi(k+1)-za(kt))/dza(kt)
2587 th2=th*th
2588 qqd=0.5*(qpi(kt)-qmi(kt))
2589 dqh=qqd*th2+qmi(kt)*th
2590 zsum = zsum + th*dza(kt)
2591 qsum = qsum + dqh*dza(kt)
2592 qn(k) = qsum/zsum
2593 endif
2594 cycle intp
2595 endif
2596 !
2597 enddo intp
2598 !
2599 ! rain out
2600 sum_precip: do k=1,km
2601 if( za(k).lt.0.0 .and. za(k+1).lt.0.0 ) then
2602 precip(i) = precip(i) + qa(k)*dza(k)
2603 cycle sum_precip
2604 else if ( za(k).lt.0.0 .and. za(k+1).ge.0.0 ) then
2605 precip(i) = precip(i) + qa(k)*(0.0-za(k))
2606 exit sum_precip
2607 endif
2608 exit sum_precip
2609 enddo sum_precip
2610 !
2611 ! replace the new values
2612 if(ist.eq.1) then
2613 rql(i,:) = qn(:)
2614 precip1(i) = precip(i)
2615 else
2616 rql2(i,:) = qn(:)
2617 precip2(i) = precip(i)
2618 endif
2619 enddo ist_loop
2620 !
2621 ! ----------------------------------
2622 enddo i_loop
2623 !
2624 END SUBROUTINE nislfv_rain_plm6
2625 END MODULE module_mp_wdm6