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