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