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