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