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