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