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