| blob | 08371a190a46e91cfc5ff980f5b0c991b42da7ed |
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_wsm3
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 :: avtr = 841.9 ! a constant for terminal velocity of rain
15 REAL, PARAMETER, PRIVATE :: bvtr = 0.8 ! a constant for terminal velocity of rain
16 REAL, PARAMETER, PRIVATE :: r0 = .8e-5 ! 8 microm in contrast to 10 micro m
17 REAL, PARAMETER, PRIVATE :: peaut = .55 ! collection efficiency
18 REAL, PARAMETER, PRIVATE :: xncr = 3.e8 ! maritime cloud in contrast to 3.e8 in tc80
19 REAL, PARAMETER, PRIVATE :: xmyu = 1.718e-5 ! the dynamic viscosity kgm-1s-1
20 REAL, PARAMETER, PRIVATE :: avts = 11.72 ! a constant for terminal velocity of snow
21 REAL, PARAMETER, PRIVATE :: bvts = .41 ! a constant for terminal velocity of snow
22 REAL, PARAMETER, PRIVATE :: n0smax = 1.e11 ! maximum n0s (t=-90C unlimited)
23 REAL, PARAMETER, PRIVATE :: lamdarmax = 8.e4 ! limited maximum value for slope parameter of rain
24 REAL, PARAMETER, PRIVATE :: lamdasmax = 1.e5 ! limited maximum value for slope parameter of snow
25 REAL, PARAMETER, PRIVATE :: lamdagmax = 6.e4 ! limited maximum value for slope parameter of graupel
26 REAL, PARAMETER, PRIVATE :: dicon = 11.9 ! constant for the cloud-ice diamter
27 REAL, PARAMETER, PRIVATE :: dimax = 500.e-6 ! limited maximum value for the cloud-ice diamter
28 REAL, PARAMETER, PRIVATE :: n0s = 2.e6 ! temperature dependent intercept parameter snow
29 REAL, PARAMETER, PRIVATE :: alpha = .12 ! .122 exponen factor for n0s
30 REAL, PARAMETER, PRIVATE :: qcrmin = 1.e-9 ! minimun values for qr, qs, and qg
31 REAL, SAVE :: &
32 qc0, qck1,bvtr1,bvtr2,bvtr3,bvtr4,g1pbr, &
33 g3pbr,g4pbr,g5pbro2,pvtr,eacrr,pacrr, &
34 precr1,precr2,xmmax,roqimax,bvts1, &
35 bvts2,bvts3,bvts4,g1pbs,g3pbs,g4pbs, &
36 g5pbso2,pvts,pacrs,precs1,precs2,pidn0r, &
37 pidn0s,xlv1,pi, &
38 rslopermax,rslopesmax,rslopegmax, &
39 rsloperbmax,rslopesbmax,rslopegbmax, &
40 rsloper2max,rslopes2max,rslopeg2max, &
41 rsloper3max,rslopes3max,rslopeg3max
42 !
43 ! Specifies code-inlining of fpvs function in WSM32D below. JM 20040507
44 !
45 CONTAINS
46 !===================================================================
47 !
48 SUBROUTINE wsm3(th, q, qci, qrs &
49 , w, den, pii, p, delz &
50 , delt,g, cpd, cpv, rd, rv, t0c &
51 , ep1, ep2, qmin &
52 , XLS, XLV0, XLF0, den0, denr &
53 , cliq,cice,psat &
54 , rain, rainncv &
55 , snow, snowncv &
56 , sr &
57 , ids,ide, jds,jde, kds,kde &
58 , ims,ime, jms,jme, kms,kme &
59 , its,ite, jts,jte, kts,kte &
60 )
61 !-------------------------------------------------------------------
62 IMPLICIT NONE
63 !-------------------------------------------------------------------
64 !
65 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
66 ims,ime, jms,jme, kms,kme , &
67 its,ite, jts,jte, kts,kte
68 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
69 INTENT(INOUT) :: &
70 th, &
71 q, &
72 qci, &
73 qrs
74 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
75 INTENT(IN ) :: w, &
76 den, &
77 pii, &
78 p, &
79 delz
80 REAL, INTENT(IN ) :: delt, &
81 g, &
82 rd, &
83 rv, &
84 t0c, &
85 den0, &
86 cpd, &
87 cpv, &
88 ep1, &
89 ep2, &
90 qmin, &
91 XLS, &
92 XLV0, &
93 XLF0, &
94 cliq, &
95 cice, &
96 psat, &
97 denr
98 REAL, DIMENSION( ims:ime , jms:jme ), &
99 INTENT(INOUT) :: rain, &
100 rainncv
101 REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
102 INTENT(INOUT) :: snow, &
103 snowncv, &
104 sr
105 ! LOCAL VAR
106 REAL, DIMENSION( its:ite , kts:kte ) :: t
107 INTEGER :: i,j,k
108 !-------------------------------------------------------------------
109 DO j=jts,jte
110 DO k=kts,kte
111 DO i=its,ite
112 t(i,k)=th(i,k,j)*pii(i,k,j)
113 ENDDO
114 ENDDO
115 CALL wsm32D(t, q(ims,kms,j), qci(ims,kms,j) &
116 ,qrs(ims,kms,j),w(ims,kms,j), den(ims,kms,j) &
117 ,p(ims,kms,j), delz(ims,kms,j) &
118 ,delt,g, cpd, cpv, rd, rv, t0c &
119 ,ep1, ep2, qmin &
120 ,XLS, XLV0, XLF0, den0, denr &
121 ,cliq,cice,psat &
122 ,j &
123 ,rain(ims,j), rainncv(ims,j) &
124 ,snow(ims,j),snowncv(ims,j) &
125 ,sr(ims,j) &
126 ,ids,ide, jds,jde, kds,kde &
127 ,ims,ime, jms,jme, kms,kme &
128 ,its,ite, jts,jte, kts,kte &
129 )
130 DO K=kts,kte
131 DO I=its,ite
132 th(i,k,j)=t(i,k)/pii(i,k,j)
133 ENDDO
134 ENDDO
135 ENDDO
136 END SUBROUTINE wsm3
137 !===================================================================
138 !
139 SUBROUTINE wsm32D(t, q &
140 ,qci, qrs,w, den, p, delz &
141 ,delt,g, cpd, cpv, rd, rv, t0c &
142 ,ep1, ep2, qmin &
143 ,XLS, XLV0, XLF0, den0, denr &
144 ,cliq,cice,psat &
145 ,lat &
146 ,rain, rainncv &
147 ,snow,snowncv &
148 ,sr &
149 ,ids,ide, jds,jde, kds,kde &
150 ,ims,ime, jms,jme, kms,kme &
151 ,its,ite, jts,jte, kts,kte &
152 )
153 !-------------------------------------------------------------------
154 IMPLICIT NONE
155 !-------------------------------------------------------------------
156 !
157 ! This code is a 3-class simple ice microphyiscs scheme (WSM3) of the
158 ! Single-Moment MicroPhyiscs (WSMMP). The WSMMP assumes that ice nuclei
159 ! number concentration is a function of temperature, and seperate assumption
160 ! is developed, in which ice crystal number concentration is a function
161 ! of ice amount. A theoretical background of the ice-microphysics and related
162 ! processes in the WSMMPs are described in Hong et al. (2004).
163 ! Production terms in the WSM6 scheme are described in Hong and Lim (2006).
164 ! All units are in m.k.s. and source/sink terms in kgkg-1s-1.
165 !
166 ! WSM3 cloud scheme
167 !
168 ! Developed by Song-You Hong (Yonsei Univ.), Jimy Dudhia (NCAR)
169 ! and Shu-Hua Chen (UC Davis)
170 ! Summer 2002
171 !
172 ! Implemented by Song-You Hong (Yonsei Univ.) and Jimy Dudhia (NCAR)
173 ! Summer 2003
174 !
175 ! History : semi-lagrangian scheme sedimentation(JH), and clean up
176 ! Hong, August 2009
177 !
178 ! Reference) Hong, Dudhia, Chen (HDC, 2004) Mon. Wea. Rev.
179 ! Dudhia (D89, 1989) J. Atmos. Sci.
180 ! Hong and Lim (HL, 2006) J. Korean Meteor. Soc.
181 ! Juang and Hong (JH, 2010) Mon. Wea. Rev.
182 !
183 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
184 ims,ime, jms,jme, kms,kme, &
185 its,ite, jts,jte, kts,kte, &
186 lat
187 REAL, DIMENSION( its:ite , kts:kte ), &
188 INTENT(INOUT) :: &
189 t
190 REAL, DIMENSION( ims:ime , kms:kme ), &
191 INTENT(INOUT) :: &
192 q, &
193 qci, &
194 qrs
195 REAL, DIMENSION( ims:ime , kms:kme ), &
196 INTENT(IN ) :: w, &
197 den, &
198 p, &
199 delz
200 REAL, INTENT(IN ) :: delt, &
201 g, &
202 cpd, &
203 cpv, &
204 t0c, &
205 den0, &
206 rd, &
207 rv, &
208 ep1, &
209 ep2, &
210 qmin, &
211 XLS, &
212 XLV0, &
213 XLF0, &
214 cliq, &
215 cice, &
216 psat, &
217 denr
218 REAL, DIMENSION( ims:ime ), &
219 INTENT(INOUT) :: rain, &
220 rainncv
221 REAL, DIMENSION( ims:ime ), OPTIONAL, &
222 INTENT(INOUT) :: snow, &
223 snowncv, &
224 sr
225 ! LOCAL VAR
226 REAL, DIMENSION( its:ite , kts:kte ) :: &
227 rh, &
228 qs, &
229 denfac, &
230 rslope, &
231 rslope2, &
232 rslope3, &
233 qrs_tmp, &
234 den_tmp, &
235 delz_tmp, &
236 rslopeb
237 REAL, DIMENSION( its:ite , kts:kte ) :: &
238 pgen, &
239 pisd, &
240 paut, &
241 pacr, &
242 pres, &
243 pcon
244 REAL, DIMENSION( its:ite , kts:kte ) :: &
245 fall, &
246 falk, &
247 xl, &
248 cpm, &
249 work1, &
250 work2, &
251 xni, &
252 qs0, &
253 denqci, &
254 denqrs, &
255 n0sfac, &
256 falkc, &
257 work1c, &
258 work2c, &
259 fallc
260 REAL, DIMENSION( its:ite ) :: delqrs,&
261 delqi
262 INTEGER, DIMENSION( its:ite ) :: kwork1,&
263 kwork2
264 INTEGER, DIMENSION( its:ite ) :: mstep, &
265 numdt
266 LOGICAL, DIMENSION( its:ite ) :: flgcld
267 REAL :: &
268 cpmcal, xlcal, diffus, &
269 viscos, xka, venfac, conden, diffac, &
270 x, y, z, a, b, c, d, e, &
271 fallsum, fallsum_qsi, vt2i,vt2s,acrfac, &
272 qdt, pvt, qik, delq, facq, qrsci, frzmlt, &
273 snomlt, hold, holdrs, facqci, supcol, coeres, &
274 supsat, dtcld, xmi, qciik, delqci, eacrs, satdt, &
275 qimax, diameter, xni0, roqi0, supice,holdc, holdci
276 INTEGER :: i, j, k, mstepmax, &
277 iprt, latd, lond, loop, loops, ifsat, kk, n, idim, kdim
278 ! Temporaries used for inlining fpvs function
279 REAL :: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
280 ! variables for optimization
281 REAL, DIMENSION( its:ite ) :: tvec1
282 !
283 !=================================================================
284 ! compute internal functions
285 !
286 cpmcal(x) = cpd*(1.-max(x,qmin))+max(x,qmin)*cpv
287 xlcal(x) = xlv0-xlv1*(x-t0c)
288 !----------------------------------------------------------------
289 ! diffus: diffusion coefficient of the water vapor
290 ! viscos: kinematic viscosity(m2s-1)
291 ! Optimizatin : A**B => exp(log(A)*(B))
292 !
293 diffus(x,y) = 8.794e-5 * exp(log(x)*(1.81)) / y ! 8.794e-5*x**1.81/y
294 viscos(x,y) = 1.496e-6 * (x*sqrt(x)) /(x+120.)/y ! 1.496e-6*x**1.5/(x+120.)/y
295 xka(x,y) = 1.414e3*viscos(x,y)*y
296 diffac(a,b,c,d,e) = d*a*a/(xka(c,d)*rv*c*c)+1./(e*diffus(c,b))
297 venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333))) &
298 /sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
299 conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+d*d/(rv*e)*c/(a*a))
300 !
301 idim = ite-its+1
302 kdim = kte-kts+1
303 !
304 !----------------------------------------------------------------
305 ! paddint 0 for negative values generated by dynamics
306 !
307 do k = kts, kte
308 do i = its, ite
309 qci(i,k) = max(qci(i,k),0.0)
310 qrs(i,k) = max(qrs(i,k),0.0)
311 enddo
312 enddo
313 !
314 !----------------------------------------------------------------
315 ! latent heat for phase changes and heat capacity. neglect the
316 ! changes during microphysical process calculation
317 ! emanuel(1994)
318 !
319 do k = kts, kte
320 do i = its, ite
321 cpm(i,k) = cpmcal(q(i,k))
322 xl(i,k) = xlcal(t(i,k))
323 enddo
324 enddo
325 do k = kts, kte
326 do i = its, ite
327 delz_tmp(i,k) = delz(i,k)
328 den_tmp(i,k) = den(i,k)
329 enddo
330 enddo
331 !
332 !----------------------------------------------------------------
333 ! compute the minor time steps.
334 !
335 loops = max(nint(delt/dtcldcr),1)
336 dtcld = delt/loops
337 if(delt.le.dtcldcr) dtcld = delt
338 !
339 do loop = 1,loops
340 !
341 !----------------------------------------------------------------
342 ! initialize the large scale variables
343 !
344 do i = its, ite
345 flgcld(i) = .true.
346 enddo
347 !
348 do k = kts, kte
349 CALL VREC( tvec1(its), den(its,k), ite-its+1)
350 do i = its, ite
351 tvec1(i) = tvec1(i)*den0
352 enddo
353 CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
354 enddo
355 !
356 ! Inline expansion for fpvs
357 ! qs(i,k) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
358 ! qs0(i,k) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
359 cvap = cpv
360 hvap=xlv0
361 hsub=xls
362 ttp=t0c+0.01
363 dldt=cvap-cliq
364 xa=-dldt/rv
365 xb=xa+hvap/(rv*ttp)
366 dldti=cvap-cice
367 xai=-dldti/rv
368 xbi=xai+hsub/(rv*ttp)
369 do k = kts, kte
370 do i = its, ite
371 tr=ttp/t(i,k)
372 if(t(i,k).lt.ttp) then
373 qs(i,k) =psat*(exp(log(tr)*(xai)))*exp(xbi*(1.-tr))
374 else
375 qs(i,k) =psat*(exp(log(tr)*(xa)))*exp(xb*(1.-tr))
376 endif
377 qs0(i,k) =psat*(exp(log(tr)*(xa)))*exp(xb*(1.-tr))
378 qs0(i,k) = (qs0(i,k)-qs(i,k))/qs(i,k)
379 qs(i,k) = ep2 * qs(i,k) / (p(i,k) - qs(i,k))
380 qs(i,k) = max(qs(i,k),qmin)
381 rh(i,k) = max(q(i,k) / qs(i,k),qmin)
382 enddo
383 enddo
384 !
385 !----------------------------------------------------------------
386 ! initialize the variables for microphysical physics
387 !
388 !
389 do k = kts, kte
390 do i = its, ite
391 pres(i,k) = 0.
392 paut(i,k) = 0.
393 pacr(i,k) = 0.
394 pgen(i,k) = 0.
395 pisd(i,k) = 0.
396 pcon(i,k) = 0.
397 fall(i,k) = 0.
398 falk(i,k) = 0.
399 fallc(i,k) = 0.
400 falkc(i,k) = 0.
401 xni(i,k) = 1.e3
402 enddo
403 enddo
404 !-------------------------------------------------------------
405 ! Ni: ice crystal number concentraiton [HDC 5c]
406 !-------------------------------------------------------------
407 do k = kts, kte
408 do i = its, ite
409 xni(i,k) = min(max(5.38e7 &
410 *exp(log((den(i,k)*max(qci(i,k),qmin)))*(0.75)),1.e3),1.e6)
411 enddo
412 enddo
413 !
414 !----------------------------------------------------------------
415 ! compute the fallout term:
416 ! first, vertical terminal velosity for minor loops
417 !---------------------------------------------------------------
418 do k = kts, kte
419 do i = its, ite
420 qrs_tmp(i,k) = qrs(i,k)
421 enddo
422 enddo
423 call slope_wsm3(qrs_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2,rslope3, &
424 work1,its,ite,kts,kte)
425 !
426 !
427 ! forward semi-laglangian scheme (JH), PCM (piecewise constant), (linear)
428 !
429 do k = kte, kts, -1
430 do i = its, ite
431 denqrs(i,k) = den(i,k)*qrs(i,k)
432 enddo
433 enddo
434 call nislfv_rain_plm(idim,kdim,den_tmp,denfac,t,delz_tmp,work1,denqrs, &
435 delqrs,dtcld,1,1)
436 do k = kts, kte
437 do i = its, ite
438 qrs(i,k) = max(denqrs(i,k)/den(i,k),0.)
439 fall(i,k) = denqrs(i,k)*work1(i,k)/delz(i,k)
440 enddo
441 enddo
442 do i = its, ite
443 fall(i,1) = delqrs(i)/delz(i,1)/dtcld
444 enddo
445 !---------------------------------------------------------------
446 ! Vice [ms-1] : fallout of ice crystal [HDC 5a]
447 !---------------------------------------------------------------
448 do k = kte, kts, -1
449 do i = its, ite
450 if(t(i,k).lt.t0c.and.qci(i,k).gt.0.) then
451 xmi = den(i,k)*qci(i,k)/xni(i,k)
452 diameter = max(dicon * sqrt(xmi), 1.e-25)
453 work1c(i,k) = 1.49e4*exp(log(diameter)*(1.31))
454 else
455 work1c(i,k) = 0.
456 endif
457 enddo
458 enddo
459 !
460 ! forward semi-laglangian scheme (JH), PCM (piecewise constant), (linear)
461 !
462 do k = kte, kts, -1
463 do i = its, ite
464 denqci(i,k) = den(i,k)*qci(i,k)
465 enddo
466 enddo
467 call nislfv_rain_plm(idim,kdim,den_tmp,denfac,t,delz_tmp,work1c,denqci, &
468 delqi,dtcld,1,0)
469 do k = kts, kte
470 do i = its, ite
471 qci(i,k) = max(denqci(i,k)/den(i,k),0.)
472 enddo
473 enddo
474 do i = its, ite
475 fallc(i,1) = delqi(i)/delz(i,1)/dtcld
476 enddo
477 !
478 !----------------------------------------------------------------
479 ! compute the freezing/melting term. [D89 B16-B17]
480 ! freezing occurs one layer above the melting level
481 !
482 do i = its, ite
483 mstep(i) = 0
484 enddo
485 do k = kts, kte
486 !
487 do i = its, ite
488 if(t(i,k).ge.t0c) then
489 mstep(i) = k
490 endif
491 enddo
492 enddo
493 !
494 do i = its, ite
495 kwork2(i) = mstep(i)
496 kwork1(i) = mstep(i)
497 if(mstep(i).ne.0) then
498 if (w(i,mstep(i)).gt.0.) then
499 kwork1(i) = mstep(i) + 1
500 endif
501 endif
502 enddo
503 !
504 do i = its, ite
505 k = kwork1(i)
506 kk = kwork2(i)
507 if(k*kk.ge.1) then
508 qrsci = qrs(i,k) + qci(i,k)
509 if(qrsci.gt.0..or.fall(i,kk).gt.0.) then
510 frzmlt = min(max(-w(i,k)*qrsci/delz(i,k),-qrsci/dtcld), &
511 qrsci/dtcld)
512 snomlt = min(max(fall(i,kk)/den(i,kk),-qrs(i,k)/dtcld), &
513 qrs(i,k)/dtcld)
514 if(k.eq.kk) then
515 t(i,k) = t(i,k) - xlf0/cpm(i,k)*(frzmlt+snomlt)*dtcld
516 else
517 t(i,k) = t(i,k) - xlf0/cpm(i,k)*frzmlt*dtcld
518 t(i,kk) = t(i,kk) - xlf0/cpm(i,kk)*snomlt*dtcld
519 endif
520 endif
521 endif
522 enddo
523 !
524 !----------------------------------------------------------------
525 ! rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
526 !
527 do i = its, ite
528 fallsum = fall(i,1)
529 fallsum_qsi = 0.
530 if((t0c-t(i,1)).gt.0) then
531 fallsum = fallsum+fallc(i,1)
532 fallsum_qsi = fall(i,1)+fallc(i,1)
533 endif
534 rainncv(i) = 0.
535 if(fallsum.gt.0.) then
536 rainncv(i) = fallsum*delz(i,1)/denr*dtcld*1000.
537 rain(i) = fallsum*delz(i,1)/denr*dtcld*1000. + rain(i)
538 endif
539 IF ( PRESENT (snowncv) .AND. PRESENT (snow)) THEN
540 snowncv(i) = 0.
541 if(fallsum_qsi.gt.0.) then
542 snowncv(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000.
543 snow(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. + snow(i)
544 endif
545 ENDIF
546 sr(i) = 0.
547 if(fallsum.gt.0.) sr(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. &
548 /(rainncv(i)+1.e-12)
549 enddo
550 !
551 !----------------------------------------------------------------
552 ! update the slope parameters for microphysics computation
553 !
554 do k = kts, kte
555 do i = its, ite
556 qrs_tmp(i,k) = qrs(i,k)
557 enddo
558 enddo
559 call slope_wsm3(qrs_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2,rslope3, &
560 work1,its,ite,kts,kte)
561 !
562 ! work1: the thermodynamic term in the denominator associated with
563 ! heat conduction and vapor diffusion
564 ! work2: parameter associated with the ventilation effects(y93)
565 !
566 do k = kts, kte
567 do i = its, ite
568 if(t(i,k).ge.t0c) then
569 work1(i,k) = diffac(xl(i,k),p(i,k),t(i,k),den(i,k),qs(i,k))
570 else
571 work1(i,k) = diffac(xls,p(i,k),t(i,k),den(i,k),qs(i,k))
572 endif
573 work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
574 enddo
575 enddo
576 !
577 do k = kts, kte
578 do i = its, ite
579 supsat = max(q(i,k),qmin)-qs(i,k)
580 satdt = supsat/dtcld
581 if(t(i,k).ge.t0c) then
582 !
583 !===============================================================
584 !
585 ! warm rain processes
586 !
587 ! - follows the processes in RH83 and LFO except for autoconcersion
588 !
589 !===============================================================
590 !---------------------------------------------------------------
591 ! praut: auto conversion rate from cloud to rain [HDC 16]
592 ! (C->R)
593 !---------------------------------------------------------------
594 if(qci(i,k).gt.qc0) then
595 ! paut(i,k) = qck1*qci(i,k)**(7./3.)
596 paut(i,k) = qck1*exp(log(qci(i,k))*((7./3.)))
597 paut(i,k) = min(paut(i,k),qci(i,k)/dtcld)
598 endif
599 !---------------------------------------------------------------
600 ! pracw: accretion of cloud water by rain [HL A40] [D89 B15]
601 ! (C->R)
602 !---------------------------------------------------------------
603 if(qrs(i,k).gt.qcrmin.and.qci(i,k).gt.qmin) then
604 pacr(i,k) = min(pacrr*rslope3(i,k)*rslopeb(i,k) &
605 *qci(i,k)*denfac(i,k),qci(i,k)/dtcld)
606 endif
607 !---------------------------------------------------------------
608 ! prevp: evaporation/condensation rate of rain [HDC 14]
609 ! (V->R or R->V)
610 !---------------------------------------------------------------
611 if(qrs(i,k).gt.0.) then
612 coeres = rslope2(i,k)*sqrt(rslope(i,k)*rslopeb(i,k))
613 pres(i,k) = (rh(i,k)-1.)*(precr1*rslope2(i,k) &
614 +precr2*work2(i,k)*coeres)/work1(i,k)
615 if(pres(i,k).lt.0.) then
616 pres(i,k) = max(pres(i,k),-qrs(i,k)/dtcld)
617 pres(i,k) = max(pres(i,k),satdt/2)
618 else
619 pres(i,k) = min(pres(i,k),satdt/2)
620 endif
621 endif
622 else
623 !
624 !===============================================================
625 !
626 ! cold rain processes
627 !
628 ! - follows the revised ice microphysics processes in HDC
629 ! - the processes same as in RH83 and LFO behave
630 ! following ice crystal hapits defined in HDC, inclduing
631 ! intercept parameter for snow (n0s), ice crystal number
632 ! concentration (ni), ice nuclei number concentration
633 ! (n0i), ice diameter (d)
634 !
635 !===============================================================
636 !
637 supcol = t0c-t(i,k)
638 n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
639 ifsat = 0
640 !-------------------------------------------------------------
641 ! Ni: ice crystal number concentraiton [HDC 5c]
642 !-------------------------------------------------------------
643 xni(i,k) = min(max(5.38e7 &
644 *exp(log((den(i,k)*max(qci(i,k),qmin)))*(0.75)),1.e3),1.e6)
645 eacrs = exp(0.07*(-supcol))
646 if(qrs(i,k).gt.qcrmin.and.qci(i,k).gt.qmin) then
647 xmi = den(i,k)*qci(i,k)/xni(i,k)
648 diameter = min(dicon * sqrt(xmi),dimax)
649 vt2i = 1.49e4*diameter**1.31
650 vt2s = pvts*rslopeb(i,k)*denfac(i,k)
651 !-------------------------------------------------------------
652 ! praci: Accretion of cloud ice by rain [HL A15] [LFO 25]
653 ! (T<T0: I->R)
654 !-------------------------------------------------------------
655 acrfac = 2.*rslope3(i,k)+2.*diameter*rslope2(i,k) &
656 +diameter**2*rslope(i,k)
657 pacr(i,k) = min(pi*qci(i,k)*eacrs*n0s*n0sfac(i,k) &
658 *abs(vt2s-vt2i)*acrfac/4.,qci(i,k)/dtcld)
659 endif
660 !-------------------------------------------------------------
661 ! pidep: Deposition/Sublimation rate of ice [HDC 9]
662 ! (T<T0: V->I or I->V)
663 !-------------------------------------------------------------
664 if(qci(i,k).gt.0.) then
665 xmi = den(i,k)*qci(i,k)/xni(i,k)
666 diameter = dicon * sqrt(xmi)
667 pisd(i,k) = 4.*diameter*xni(i,k)*(rh(i,k)-1.)/work1(i,k)
668 if(pisd(i,k).lt.0.) then
669 pisd(i,k) = max(pisd(i,k),satdt/2)
670 pisd(i,k) = max(pisd(i,k),-qci(i,k)/dtcld)
671 else
672 pisd(i,k) = min(pisd(i,k),satdt/2)
673 endif
674 if(abs(pisd(i,k)).ge.abs(satdt)) ifsat = 1
675 endif
676 !-------------------------------------------------------------
677 ! psdep: deposition/sublimation rate of snow [HDC 14]
678 ! (V->S or S->V)
679 !-------------------------------------------------------------
680 if(qrs(i,k).gt.0..and.ifsat.ne.1) then
681 coeres = rslope2(i,k)*sqrt(rslope(i,k)*rslopeb(i,k))
682 pres(i,k) = (rh(i,k)-1.)*n0sfac(i,k)*(precs1*rslope2(i,k) &
683 +precs2*work2(i,k)*coeres)/work1(i,k)
684 supice = satdt-pisd(i,k)
685 if(pres(i,k).lt.0.) then
686 pres(i,k) = max(pres(i,k),-qrs(i,k)/dtcld)
687 pres(i,k) = max(max(pres(i,k),satdt/2),supice)
688 else
689 pres(i,k) = min(min(pres(i,k),satdt/2),supice)
690 endif
691 if(abs(pisd(i,k)+pres(i,k)).ge.abs(satdt)) ifsat = 1
692 endif
693 !-------------------------------------------------------------
694 ! pigen: generation(nucleation) of ice from vapor [HDC 7-8]
695 ! (T<T0: V->I)
696 !-------------------------------------------------------------
697 if(supsat.gt.0.and.ifsat.ne.1) then
698 supice = satdt-pisd(i,k)-pres(i,k)
699 xni0 = 1.e3*exp(0.1*supcol)
700 roqi0 = 4.92e-11*exp(log(xni0)*(1.33))
701 pgen(i,k) = max(0.,(roqi0/den(i,k)-max(qci(i,k),0.))/dtcld)
702 pgen(i,k) = min(min(pgen(i,k),satdt),supice)
703 endif
704 !-------------------------------------------------------------
705 ! psaut: conversion(aggregation) of ice to snow [HDC 12]
706 ! (T<T0: I->S)
707 !-------------------------------------------------------------
708 if(qci(i,k).gt.0.) then
709 qimax = roqimax/den(i,k)
710 paut(i,k) = max(0.,(qci(i,k)-qimax)/dtcld)
711 endif
712 endif
713 enddo
714 enddo
715 !
716 !----------------------------------------------------------------
717 ! check mass conservation of generation terms and feedback to the
718 ! large scale
719 !
720 do k = kts, kte
721 do i = its, ite
722 qciik = max(qmin,qci(i,k))
723 delqci = (paut(i,k)+pacr(i,k)-pgen(i,k)-pisd(i,k))*dtcld
724 if(delqci.ge.qciik) then
725 facqci = qciik/delqci
726 paut(i,k) = paut(i,k)*facqci
727 pacr(i,k) = pacr(i,k)*facqci
728 pgen(i,k) = pgen(i,k)*facqci
729 pisd(i,k) = pisd(i,k)*facqci
730 endif
731 qik = max(qmin,q(i,k))
732 delq = (pres(i,k)+pgen(i,k)+pisd(i,k))*dtcld
733 if(delq.ge.qik) then
734 facq = qik/delq
735 pres(i,k) = pres(i,k)*facq
736 pgen(i,k) = pgen(i,k)*facq
737 pisd(i,k) = pisd(i,k)*facq
738 endif
739 work2(i,k) = -pres(i,k)-pgen(i,k)-pisd(i,k)
740 q(i,k) = q(i,k)+work2(i,k)*dtcld
741 qci(i,k) = max(qci(i,k)-(paut(i,k)+pacr(i,k)-pgen(i,k)-pisd(i,k)) &
742 *dtcld,0.)
743 qrs(i,k) = max(qrs(i,k)+(paut(i,k)+pacr(i,k)+pres(i,k))*dtcld,0.)
744 if(t(i,k).lt.t0c) then
745 t(i,k) = t(i,k)-xls*work2(i,k)/cpm(i,k)*dtcld
746 else
747 t(i,k) = t(i,k)-xl(i,k)*work2(i,k)/cpm(i,k)*dtcld
748 endif
749 enddo
750 enddo
751 !
752 cvap = cpv
753 hvap = xlv0
754 hsub = xls
755 ttp=t0c+0.01
756 dldt=cvap-cliq
757 xa=-dldt/rv
758 xb=xa+hvap/(rv*ttp)
759 dldti=cvap-cice
760 xai=-dldti/rv
761 xbi=xai+hsub/(rv*ttp)
762 do k = kts, kte
763 do i = its, ite
764 tr=ttp/t(i,k)
765 qs(i,k)=psat*(exp(log(tr)*(xa)))*exp(xb*(1.-tr))
766 qs(i,k) = ep2 * qs(i,k) / (p(i,k) - qs(i,k))
767 qs(i,k) = max(qs(i,k),qmin)
768 denfac(i,k) = sqrt(den0/den(i,k))
769 enddo
770 enddo
771 !
772 !----------------------------------------------------------------
773 ! pcond: condensational/evaporational rate of cloud water [HL A46] [RH83 A6]
774 ! if there exists additional water vapor condensated/if
775 ! evaporation of cloud water is not enough to remove subsaturation
776 !
777 do k = kts, kte
778 do i = its, ite
779 work1(i,k) = conden(t(i,k),q(i,k),qs(i,k),xl(i,k),cpm(i,k))
780 work2(i,k) = qci(i,k)+work1(i,k)
781 pcon(i,k) = min(max(work1(i,k),0.),max(q(i,k),0.))/dtcld
782 if(qci(i,k).gt.0..and.work1(i,k).lt.0.and.t(i,k).gt.t0c) &
783 pcon(i,k) = max(work1(i,k),-qci(i,k))/dtcld
784 q(i,k) = q(i,k)-pcon(i,k)*dtcld
785 qci(i,k) = max(qci(i,k)+pcon(i,k)*dtcld,0.)
786 t(i,k) = t(i,k)+pcon(i,k)*xl(i,k)/cpm(i,k)*dtcld
787 enddo
788 enddo
789 !
790 !----------------------------------------------------------------
791 ! padding for small values
792 !
793 do k = kts, kte
794 do i = its, ite
795 if(qci(i,k).le.qmin) qci(i,k) = 0.0
796 if(qrs(i,k).le.qcrmin) qrs(i,k) = 0.0
797 enddo
798 enddo
799 !
800 enddo ! big loops
801 END SUBROUTINE wsm32D
802 ! ...................................................................
803 REAL FUNCTION rgmma(x)
804 !-------------------------------------------------------------------
805 IMPLICIT NONE
806 !-------------------------------------------------------------------
807 ! rgmma function: use infinite product form
808 REAL :: euler
809 PARAMETER (euler=0.577215664901532)
810 REAL :: x, y
811 INTEGER :: i
812 if(x.eq.1.)then
813 rgmma=0.
814 else
815 rgmma=x*exp(euler*x)
816 do i=1,10000
817 y=float(i)
818 rgmma=rgmma*(1.000+x/y)*exp(-x/y)
819 enddo
820 rgmma=1./rgmma
821 endif
822 END FUNCTION rgmma
823 !
824 !--------------------------------------------------------------------------
825 REAL FUNCTION fpvs(t,ice,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c)
826 !--------------------------------------------------------------------------
827 IMPLICIT NONE
828 !--------------------------------------------------------------------------
829 REAL t,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c,dldt,xa,xb,dldti, &
830 xai,xbi,ttp,tr
831 INTEGER ice
832 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
833 ttp=t0c+0.01
834 dldt=cvap-cliq
835 xa=-dldt/rv
836 xb=xa+hvap/(rv*ttp)
837 dldti=cvap-cice
838 xai=-dldti/rv
839 xbi=xai+hsub/(rv*ttp)
840 tr=ttp/t
841 if(t.lt.ttp.and.ice.eq.1) then
842 fpvs=psat*(tr**xai)*exp(xbi*(1.-tr))
843 else
844 fpvs=psat*(tr**xa)*exp(xb*(1.-tr))
845 endif
846 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
847 END FUNCTION fpvs
848 !-------------------------------------------------------------------
849 SUBROUTINE wsm3init(den0,denr,dens,cl,cpv,allowed_to_read)
850 !-------------------------------------------------------------------
851 IMPLICIT NONE
852 !-------------------------------------------------------------------
853 !.... constants which may not be tunable
854 REAL, INTENT(IN) :: den0,denr,dens,cl,cpv
855 LOGICAL, INTENT(IN) :: allowed_to_read
856 !
857 pi = 4.*atan(1.)
858 xlv1 = cl-cpv
859 !
860 qc0 = 4./3.*pi*denr*r0**3*xncr/den0 ! 0.419e-3 -- .61e-3
861 qck1 = .104*9.8*peaut/(xncr*denr)**(1./3.)/xmyu*den0**(4./3.) ! 7.03
862 !
863 bvtr1 = 1.+bvtr
864 bvtr2 = 2.5+.5*bvtr
865 bvtr3 = 3.+bvtr
866 bvtr4 = 4.+bvtr
867 g1pbr = rgmma(bvtr1)
868 g3pbr = rgmma(bvtr3)
869 g4pbr = rgmma(bvtr4) ! 17.837825
870 g5pbro2 = rgmma(bvtr2) ! 1.8273
871 pvtr = avtr*g4pbr/6.
872 eacrr = 1.0
873 pacrr = pi*n0r*avtr*g3pbr*.25*eacrr
874 precr1 = 2.*pi*n0r*.78
875 precr2 = 2.*pi*n0r*.31*avtr**.5*g5pbro2
876 xmmax = (dimax/dicon)**2
877 roqimax = 2.08e22*dimax**8
878 !
879 bvts1 = 1.+bvts
880 bvts2 = 2.5+.5*bvts
881 bvts3 = 3.+bvts
882 bvts4 = 4.+bvts
883 g1pbs = rgmma(bvts1) !.8875
884 g3pbs = rgmma(bvts3)
885 g4pbs = rgmma(bvts4) ! 12.0786
886 g5pbso2 = rgmma(bvts2)
887 pvts = avts*g4pbs/6.
888 pacrs = pi*n0s*avts*g3pbs*.25
889 precs1 = 4.*n0s*.65
890 precs2 = 4.*n0s*.44*avts**.5*g5pbso2
891 pidn0r = pi*denr*n0r
892 pidn0s = pi*dens*n0s
893 !
894 rslopermax = 1./lamdarmax
895 rslopesmax = 1./lamdasmax
896 rsloperbmax = rslopermax ** bvtr
897 rslopesbmax = rslopesmax ** bvts
898 rsloper2max = rslopermax * rslopermax
899 rslopes2max = rslopesmax * rslopesmax
900 rsloper3max = rsloper2max * rslopermax
901 rslopes3max = rslopes2max * rslopesmax
902 !
903 END SUBROUTINE wsm3init
904 !
905 subroutine slope_wsm3(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3,vt,its,ite,kts,kte)
906 IMPLICIT NONE
907 INTEGER :: its,ite, jts,jte, kts,kte
908 REAL, DIMENSION( its:ite , kts:kte ) :: &
909 qrs, &
910 den, &
911 denfac, &
912 t, &
913 rslope, &
914 rslopeb, &
915 rslope2, &
916 rslope3, &
917 vt
918 REAL, PARAMETER :: t0c = 273.15
919 REAL, DIMENSION( its:ite , kts:kte ) :: &
920 n0sfac
921 REAL :: lamdar,lamdas,x, y, z, supcol, pvt
922 integer :: i, j, k
923 !----------------------------------------------------------------
924 ! size distributions: (x=mixing ratio, y=air density):
925 ! valid for mixing ratio > 1.e-9 kg/kg.
926 !
927 lamdar(x,y)= sqrt(sqrt(pidn0r/(x*y))) ! (pidn0r/(x*y))**.25
928 lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y))) ! (pidn0s*z/(x*y))**.25
929 !
930 do k = kts, kte
931 do i = its, ite
932 if(t(i,k).ge.t0c) then
933 pvt = pvtr
934 if(qrs(i,k).le.qcrmin)then
935 rslope(i,k) = rslopermax
936 rslopeb(i,k) = rsloperbmax
937 rslope2(i,k) = rsloper2max
938 rslope3(i,k) = rsloper3max
939 else
940 rslope(i,k) = 1./lamdar(qrs(i,k),den(i,k))
941 rslopeb(i,k) = exp(log(rslope(i,k))*(bvtr))
942 rslope2(i,k) = rslope(i,k)*rslope(i,k)
943 rslope3(i,k) = rslope2(i,k)*rslope(i,k)
944 endif
945 else
946 supcol = t0c-t(i,k)
947 n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
948 pvt = pvts
949 if(qrs(i,k).le.qcrmin)then
950 rslope(i,k) = rslopesmax
951 rslopeb(i,k) = rslopesbmax
952 rslope2(i,k) = rslopes2max
953 rslope3(i,k) = rslopes3max
954 else
955 rslope(i,k) = 1./lamdas(qrs(i,k),den(i,k),n0sfac(i,k))
956 rslopeb(i,k) = exp(log(rslope(i,k))*(bvts))
957 rslope2(i,k) = rslope(i,k)*rslope(i,k)
958 rslope3(i,k) = rslope2(i,k)*rslope(i,k)
959 endif
960 endif
961 vt(i,k) = pvt*rslopeb(i,k)*denfac(i,k)
962 if(qrs(i,k).le.0.0) vt(i,k) = 0.0
963 enddo
964 enddo
965 END subroutine slope_wsm3
966 !-------------------------------------------------------------------
967 SUBROUTINE nislfv_rain_pcm(im,km,denl,denfacl,tkl,dzl,wwl,rql,precip,dt,id,iter)
968 !-------------------------------------------------------------------
969 !
970 ! for non-iteration semi-Lagrangain forward advection for cloud
971 ! with mass conservation and positive definite advection
972 ! 2nd order interpolation with monotonic piecewise linear method
973 ! this routine is under assumption of decfl < 1 for semi_Lagrangian
974 !
975 ! dzl depth of model layer in meter
976 ! wwl terminal velocity at model layer m/s
977 ! rql cloud density*mixing ration
978 ! precip precipitation
979 ! dt time step
980 ! id kind of precip: 0 test case; 1 raindrop
981 ! iter how many time to guess mean terminal velocity: 0 pure forward.
982 ! 0 : use departure wind for advection
983 ! 1 : use mean wind for advection
984 ! > 1 : use mean wind after iter-1 iterations
985 !
986 ! author: hann-ming henry juang <henry.juang@noaa.gov>
987 ! implemented by song-you hong
988 !
989 implicit none
990 integer im,km,id
991 real dt
992 real dzl(im,km),wwl(im,km),rql(im,km),precip(im)
993 real denl(im,km),denfacl(im,km),tkl(im,km)
994 !
995 integer i,k,n,m,kk,kb,kt,iter
996 real tl,tl2,qql,dql,qqd
997 real th,th2,qqh,dqh
998 real zsum,qsum,dim,dip,c1,con1,fa1,fa2
999 real zsumt,qsumt,zsumb,qsumb
1000 real allold, allnew, zz, dzamin, cflmax, decfl
1001 real dz(km), ww(km), qq(km), wd(km), wa(km), was(km)
1002 real den(km), denfac(km), tk(km)
1003 real wi(km+1), zi(km+1), za(km+1)
1004 real qn(km), qr(km),tmp(km),tmp1(km),tmp2(km),tmp3(km)
1005 real dza(km+1), qa(km+1), qmi(km+1), qpi(km+1)
1006 !
1007 precip(:) = 0.0
1008 !
1009 i_loop : do i=1,im
1010 ! -----------------------------------
1011 dz(:) = dzl(i,:)
1012 qq(:) = rql(i,:)
1013 ww(:) = wwl(i,:)
1014 den(:) = denl(i,:)
1015 denfac(:) = denfacl(i,:)
1016 tk(:) = tkl(i,:)
1017 ! skip for no precipitation for all layers
1018 allold = 0.0
1019 do k=1,km
1020 allold = allold + qq(k)
1021 enddo
1022 if(allold.le.0.0) then
1023 cycle i_loop
1024 endif
1025 !
1026 ! compute interface values
1027 zi(1)=0.0
1028 do k=1,km
1029 zi(k+1) = zi(k)+dz(k)
1030 enddo
1031 !
1032 ! save departure wind
1033 wd(:) = ww(:)
1034 n=1
1035 100 continue
1036 ! pcm is 1st order, we should use 2nd order wi
1037 ! 2nd order interpolation to get wi
1038 wi(1) = ww(1)
1039 do k=2,km
1040 wi(k) = (ww(k)*dz(k-1)+ww(k-1)*dz(k))/(dz(k-1)+dz(k))
1041 enddo
1042 wi(km+1) = ww(km)
1043 !
1044 ! terminate of top of raingroup
1045 do k=2,km
1046 if( ww(k).eq.0.0 ) wi(k)=ww(k-1)
1047 enddo
1048 !
1049 ! diffusivity of wi
1050 con1 = 0.05
1051 do k=km,1,-1
1052 decfl = (wi(k+1)-wi(k))*dt/dz(k)
1053 if( decfl .gt. con1 ) then
1054 wi(k) = wi(k+1) - con1*dz(k)/dt
1055 endif
1056 enddo
1057 ! compute arrival point
1058 do k=1,km+1
1059 za(k) = zi(k) - wi(k)*dt
1060 enddo
1061 !
1062 do k=1,km
1063 dza(k) = za(k+1)-za(k)
1064 enddo
1065 dza(km+1) = zi(km+1) - za(km+1)
1066 !
1067 ! computer deformation at arrival point
1068 do k=1,km
1069 qa(k) = qq(k)*dz(k)/dza(k)
1070 qr(k) = qa(k)/den(k)
1071 enddo
1072 qa(km+1) = 0.0
1073 ! call maxmin(km,1,qa,' arrival points ')
1074 !
1075 ! compute arrival terminal velocity, and estimate mean terminal velocity
1076 ! then back to use mean terminal velocity
1077 if( n.le.iter ) then
1078 call slope_wsm3(qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,1,1,1,km)
1079 if( n.eq.2 ) wa(1:km) = 0.5*(wa(1:km)+was(1:km))
1080 do k=1,km
1081 !#ifdef DEBUG
1082 ! print*,' slope_wsm3 ',qr(k)*1000.,den(k),denfac(k),tk(k),tmp(k),tmp1(k),tmp2(k),ww(k),wa(k)
1083 !#endif
1084 ! mean wind is average of departure and new arrival winds
1085 ww(k) = 0.5* ( wd(k)+wa(k) )
1086 enddo
1087 was(:) = wa(:)
1088 n=n+1
1089 go to 100
1090 endif
1091 !
1092 !
1093 ! interpolation to regular point
1094 qn = 0.0
1095 kb=1
1096 kt=1
1097 intp : do k=1,km
1098 kb=max(kb-1,1)
1099 kt=max(kt-1,1)
1100 ! find kb and kt
1101 if( zi(k).ge.za(km+1) ) then
1102 exit intp
1103 else
1104 find_kb : do kk=kb,km
1105 if( zi(k).le.za(kk+1) ) then
1106 kb = kk
1107 exit find_kb
1108 else
1109 cycle find_kb
1110 endif
1111 enddo find_kb
1112 find_kt : do kk=kt,km
1113 if( zi(k+1).le.za(kk) ) then
1114 kt = kk
1115 exit find_kt
1116 else
1117 cycle find_kt
1118 endif
1119 enddo find_kt
1120 ! compute q with piecewise constant method
1121 if( kt-kb.eq.1 ) then
1122 qn(k) = qa(kb)
1123 else if( kt-kb.ge.2 ) then
1124 zsumb = za(kb+1)-zi(k)
1125 qsumb = qa(kb) * zsumb
1126 zsumt = zi(k+1)-za(kt-1)
1127 qsumt = qa(kt-1) * zsumt
1128 qsum = 0.0
1129 zsum = 0.0
1130 if( kt-kb.ge.3 ) then
1131 do m=kb+1,kt-2
1132 qsum = qsum + qa(m) * dza(m)
1133 zsum = zsum + dza(m)
1134 enddo
1135 endif
1136 qn(k) = (qsumb+qsum+qsumt)/(zsumb+zsum+zsumt)
1137 endif
1138 cycle intp
1139 endif
1140 !
1141 enddo intp
1142 !
1143 ! rain out
1144 sum_precip: do k=1,km
1145 if( za(k).lt.0.0 .and. za(k+1).lt.0.0 ) then
1146 precip(i) = precip(i) + qa(k)*dza(k)
1147 cycle sum_precip
1148 else if ( za(k).lt.0.0 .and. za(k+1).ge.0.0 ) then
1149 precip(i) = precip(i) + qa(k)*(0.0-za(k))
1150 exit sum_precip
1151 endif
1152 exit sum_precip
1153 enddo sum_precip
1154 !
1155 ! replace the new values
1156 rql(i,:) = qn(:)
1157 !
1158 ! ----------------------------------
1159 enddo i_loop
1160 !
1161 END SUBROUTINE nislfv_rain_pcm
1162 !-------------------------------------------------------------------
1163 SUBROUTINE nislfv_rain_plm(im,km,denl,denfacl,tkl,dzl,wwl,rql,precip,dt,id,iter)
1164 !-------------------------------------------------------------------
1165 !
1166 ! for non-iteration semi-Lagrangain forward advection for cloud
1167 ! with mass conservation and positive definite advection
1168 ! 2nd order interpolation with monotonic piecewise linear method
1169 ! this routine is under assumption of decfl < 1 for semi_Lagrangian
1170 !
1171 ! dzl depth of model layer in meter
1172 ! wwl terminal velocity at model layer m/s
1173 ! rql cloud density*mixing ration
1174 ! precip precipitation
1175 ! dt time step
1176 ! id kind of precip: 0 test case; 1 raindrop
1177 ! iter how many time to guess mean terminal velocity: 0 pure forward.
1178 ! 0 : use departure wind for advection
1179 ! 1 : use mean wind for advection
1180 ! > 1 : use mean wind after iter-1 iterations
1181 !
1182 ! author: hann-ming henry juang <henry.juang@noaa.gov>
1183 ! implemented by song-you hong
1184 !
1185 implicit none
1186 integer im,km,id
1187 real dt
1188 real dzl(im,km),wwl(im,km),rql(im,km),precip(im)
1189 real denl(im,km),denfacl(im,km),tkl(im,km)
1190 !
1191 integer i,k,n,m,kk,kb,kt,iter
1192 real tl,tl2,qql,dql,qqd
1193 real th,th2,qqh,dqh
1194 real zsum,qsum,dim,dip,c1,con1,fa1,fa2
1195 real allold, allnew, zz, dzamin, cflmax, decfl
1196 real dz(km), ww(km), qq(km), wd(km), wa(km), was(km)
1197 real den(km), denfac(km), tk(km)
1198 real wi(km+1), zi(km+1), za(km+1)
1199 real qn(km), qr(km),tmp(km),tmp1(km),tmp2(km),tmp3(km)
1200 real dza(km+1), qa(km+1), qmi(km+1), qpi(km+1)
1201 !
1202 precip(:) = 0.0
1203 !
1204 i_loop : do i=1,im
1205 ! -----------------------------------
1206 dz(:) = dzl(i,:)
1207 qq(:) = rql(i,:)
1208 ww(:) = wwl(i,:)
1209 den(:) = denl(i,:)
1210 denfac(:) = denfacl(i,:)
1211 tk(:) = tkl(i,:)
1212 ! skip for no precipitation for all layers
1213 allold = 0.0
1214 do k=1,km
1215 allold = allold + qq(k)
1216 enddo
1217 if(allold.le.0.0) then
1218 cycle i_loop
1219 endif
1220 !
1221 ! compute interface values
1222 zi(1)=0.0
1223 do k=1,km
1224 zi(k+1) = zi(k)+dz(k)
1225 enddo
1226 !
1227 ! save departure wind
1228 wd(:) = ww(:)
1229 n=1
1230 100 continue
1231 ! plm is 2nd order, we can use 2nd order wi or 3rd order wi
1232 ! 2nd order interpolation to get wi
1233 wi(1) = ww(1)
1234 wi(km+1) = ww(km)
1235 do k=2,km
1236 wi(k) = (ww(k)*dz(k-1)+ww(k-1)*dz(k))/(dz(k-1)+dz(k))
1237 enddo
1238 ! 3rd order interpolation to get wi
1239 fa1 = 9./16.
1240 fa2 = 1./16.
1241 wi(1) = ww(1)
1242 wi(2) = 0.5*(ww(2)+ww(1))
1243 do k=3,km-1
1244 wi(k) = fa1*(ww(k)+ww(k-1))-fa2*(ww(k+1)+ww(k-2))
1245 enddo
1246 wi(km) = 0.5*(ww(km)+ww(km-1))
1247 wi(km+1) = ww(km)
1248 !
1249 ! terminate of top of raingroup
1250 do k=2,km
1251 if( ww(k).eq.0.0 ) wi(k)=ww(k-1)
1252 enddo
1253 !
1254 ! diffusivity of wi
1255 con1 = 0.05
1256 do k=km,1,-1
1257 decfl = (wi(k+1)-wi(k))*dt/dz(k)
1258 if( decfl .gt. con1 ) then
1259 wi(k) = wi(k+1) - con1*dz(k)/dt
1260 endif
1261 enddo
1262 ! compute arrival point
1263 do k=1,km+1
1264 za(k) = zi(k) - wi(k)*dt
1265 enddo
1266 !
1267 do k=1,km
1268 dza(k) = za(k+1)-za(k)
1269 enddo
1270 dza(km+1) = zi(km+1) - za(km+1)
1271 !
1272 ! computer deformation at arrival point
1273 do k=1,km
1274 qa(k) = qq(k)*dz(k)/dza(k)
1275 qr(k) = qa(k)/den(k)
1276 enddo
1277 qa(km+1) = 0.0
1278 ! call maxmin(km,1,qa,' arrival points ')
1279 !
1280 ! compute arrival terminal velocity, and estimate mean terminal velocity
1281 ! then back to use mean terminal velocity
1282 if( n.le.iter ) then
1283 call slope_wsm3(qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,1,1,1,km)
1284 if( n.ge.2 ) wa(1:km)=0.5*(wa(1:km)+was(1:km))
1285 do k=1,km
1286 !#ifdef DEBUG
1287 ! print*,' slope_wsm3 ',qr(k)*1000.,den(k),denfac(k),tk(k),tmp(k),tmp1(k),tmp2(k),ww(k),wa(k)
1288 !#endif
1289 ! mean wind is average of departure and new arrival winds
1290 ww(k) = 0.5* ( wd(k)+wa(k) )
1291 enddo
1292 was(:) = wa(:)
1293 n=n+1
1294 go to 100
1295 endif
1296 !
1297 ! estimate values at arrival cell interface with monotone
1298 do k=2,km
1299 dip=(qa(k+1)-qa(k))/(dza(k+1)+dza(k))
1300 dim=(qa(k)-qa(k-1))/(dza(k-1)+dza(k))
1301 if( dip*dim.le.0.0 ) then
1302 qmi(k)=qa(k)
1303 qpi(k)=qa(k)
1304 else
1305 qpi(k)=qa(k)+0.5*(dip+dim)*dza(k)
1306 qmi(k)=2.0*qa(k)-qpi(k)
1307 if( qpi(k).lt.0.0 .or. qmi(k).lt.0.0 ) then
1308 qpi(k) = qa(k)
1309 qmi(k) = qa(k)
1310 endif
1311 endif
1312 enddo
1313 qpi(1)=qa(1)
1314 qmi(1)=qa(1)
1315 qmi(km+1)=qa(km+1)
1316 qpi(km+1)=qa(km+1)
1317 !
1318 ! interpolation to regular point
1319 qn = 0.0
1320 kb=1
1321 kt=1
1322 intp : do k=1,km
1323 kb=max(kb-1,1)
1324 kt=max(kt-1,1)
1325 ! find kb and kt
1326 if( zi(k).ge.za(km+1) ) then
1327 exit intp
1328 else
1329 find_kb : do kk=kb,km
1330 if( zi(k).le.za(kk+1) ) then
1331 kb = kk
1332 exit find_kb
1333 else
1334 cycle find_kb
1335 endif
1336 enddo find_kb
1337 find_kt : do kk=kt,km
1338 if( zi(k+1).le.za(kk) ) then
1339 kt = kk
1340 exit find_kt
1341 else
1342 cycle find_kt
1343 endif
1344 enddo find_kt
1345 kt = kt - 1
1346 ! compute q with piecewise constant method
1347 if( kt.eq.kb ) then
1348 tl=(zi(k)-za(kb))/dza(kb)
1349 th=(zi(k+1)-za(kb))/dza(kb)
1350 tl2=tl*tl
1351 th2=th*th
1352 qqd=0.5*(qpi(kb)-qmi(kb))
1353 qqh=qqd*th2+qmi(kb)*th
1354 qql=qqd*tl2+qmi(kb)*tl
1355 qn(k) = (qqh-qql)/(th-tl)
1356 else if( kt.gt.kb ) then
1357 tl=(zi(k)-za(kb))/dza(kb)
1358 tl2=tl*tl
1359 qqd=0.5*(qpi(kb)-qmi(kb))
1360 qql=qqd*tl2+qmi(kb)*tl
1361 dql = qa(kb)-qql
1362 zsum = (1.-tl)*dza(kb)
1363 qsum = dql*dza(kb)
1364 if( kt-kb.gt.1 ) then
1365 do m=kb+1,kt-1
1366 zsum = zsum + dza(m)
1367 qsum = qsum + qa(m) * dza(m)
1368 enddo
1369 endif
1370 th=(zi(k+1)-za(kt))/dza(kt)
1371 th2=th*th
1372 qqd=0.5*(qpi(kt)-qmi(kt))
1373 dqh=qqd*th2+qmi(kt)*th
1374 zsum = zsum + th*dza(kt)
1375 qsum = qsum + dqh*dza(kt)
1376 qn(k) = qsum/zsum
1377 endif
1378 cycle intp
1379 endif
1380 !
1381 enddo intp
1382 !
1383 ! rain out
1384 sum_precip: do k=1,km
1385 if( za(k).lt.0.0 .and. za(k+1).lt.0.0 ) then
1386 precip(i) = precip(i) + qa(k)*dza(k)
1387 cycle sum_precip
1388 else if ( za(k).lt.0.0 .and. za(k+1).ge.0.0 ) then
1389 precip(i) = precip(i) + qa(k)*(0.0-za(k))
1390 exit sum_precip
1391 endif
1392 exit sum_precip
1393 enddo sum_precip
1394 !
1395 ! replace the new values
1396 rql(i,:) = qn(:)
1397 !
1398 ! ----------------------------------
1399 enddo i_loop
1400 !
1401 END SUBROUTINE nislfv_rain_plm
1402 !
1403 END MODULE module_mp_wsm3