| blob | e081a7b6e489a7bde28ebf636ea387c6d9c83ea2 |
1 #ifdef _ACCEL
2 # include "module_mp_wsm5_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 !Including inline expansion statistical function
13 MODULE module_mp_wsm5
14 !
15 USE module_mp_radar
16 USE module_model_constants, only : RE_QC_BG, RE_QI_BG, RE_QS_BG
17 !
18 REAL, PARAMETER, PRIVATE :: dtcldcr = 120. ! maximum time step for minor loops
19 REAL, PARAMETER, PRIVATE :: n0r = 8.e6 ! intercept parameter rain
20 REAL, PARAMETER, PRIVATE :: avtr = 841.9 ! a constant for terminal velocity of rain
21 REAL, PARAMETER, PRIVATE :: bvtr = 0.8 ! a constant for terminal velocity of rain
22 REAL, PARAMETER, PRIVATE :: r0 = .8e-5 ! 8 microm in contrast to 10 micro m
23 REAL, PARAMETER, PRIVATE :: peaut = .55 ! collection efficiency
24 REAL, PARAMETER, PRIVATE :: xncr = 3.e8 ! maritime cloud in contrast to 3.e8 in tc80
25 REAL, PARAMETER, PRIVATE :: xmyu = 1.718e-5 ! the dynamic viscosity kgm-1s-1
26 REAL, PARAMETER, PRIVATE :: avts = 11.72 ! a constant for terminal velocity of snow
27 REAL, PARAMETER, PRIVATE :: bvts = .41 ! a constant for terminal velocity of snow
28 REAL, PARAMETER, PRIVATE :: n0smax = 1.e11 ! maximum n0s (t=-90C unlimited)
29 REAL, PARAMETER, PRIVATE :: lamdarmax = 8.e4 ! limited maximum value for slope parameter of rain
30 REAL, PARAMETER, PRIVATE :: lamdasmax = 1.e5 ! limited maximum value for slope parameter of snow
31 REAL, PARAMETER, PRIVATE :: lamdagmax = 6.e4 ! limited maximum value for slope parameter of graupel
32 REAL, PARAMETER, PRIVATE :: dicon = 11.9 ! constant for the cloud-ice diamter
33 REAL, PARAMETER, PRIVATE :: dimax = 500.e-6 ! limited maximum value for the cloud-ice diamter
34 REAL, PARAMETER, PRIVATE :: n0s = 2.e6 ! temperature dependent intercept parameter snow
35 REAL, PARAMETER, PRIVATE :: alpha = .12 ! .122 exponen factor for n0s
36 REAL, PARAMETER, PRIVATE :: pfrz1 = 100. ! constant in Biggs freezing
37 REAL, PARAMETER, PRIVATE :: pfrz2 = 0.66 ! constant in Biggs freezing
38 REAL, PARAMETER, PRIVATE :: qcrmin = 1.e-9 ! minimun values for qr, qs, and qg
39 REAL, PARAMETER, PRIVATE :: eacrc = 1.0 ! Snow/cloud-water collection efficiency
40 REAL, SAVE :: &
41 qc0, qck1, pidnc, &
42 bvtr1,bvtr2,bvtr3,bvtr4,g1pbr, &
43 g3pbr,g4pbr,g5pbro2,pvtr,eacrr,pacrr, &
44 precr1,precr2,xmmax,roqimax,bvts1, &
45 bvts2,bvts3,bvts4,g1pbs,g3pbs,g4pbs, &
46 g5pbso2,pvts,pacrs,precs1,precs2,pidn0r, &
47 pidn0s,xlv1,pacrc,pi, &
48 rslopermax,rslopesmax,rslopegmax, &
49 rsloperbmax,rslopesbmax,rslopegbmax, &
50 rsloper2max,rslopes2max,rslopeg2max, &
51 rsloper3max,rslopes3max,rslopeg3max
52 !
53 ! Specifies code-inlining of fpvs function in WSM52D below. JM 20040507
54 !
55 CONTAINS
56 !===================================================================
57 !
58 SUBROUTINE wsm5(th, q, qc, qr, qi, qs &
59 ,den, pii, p, delz &
60 ,delt,g, cpd, cpv, rd, rv, t0c &
61 ,ep1, ep2, qmin &
62 ,XLS, XLV0, XLF0, den0, denr &
63 ,cliq,cice,psat &
64 ,rain, rainncv &
65 ,snow, snowncv &
66 ,sr &
67 ,refl_10cm, diagflag, do_radar_ref &
68 ,has_reqc, has_reqi, has_reqs & ! for radiation
69 ,re_cloud, re_ice, re_snow & ! for radiation
70 ,ids,ide, jds,jde, kds,kde &
71 ,ims,ime, jms,jme, kms,kme &
72 ,its,ite, jts,jte, kts,kte &
73 )
74 !-------------------------------------------------------------------
75 IMPLICIT NONE
76 !-------------------------------------------------------------------
77 !
78 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
79 ims,ime, jms,jme, kms,kme , &
80 its,ite, jts,jte, kts,kte
81 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
82 INTENT(INOUT) :: &
83 th, &
84 q, &
85 qc, &
86 qi, &
87 qr, &
88 qs
89 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
90 INTENT(IN ) :: &
91 den, &
92 pii, &
93 p, &
94 delz
95 REAL, INTENT(IN ) :: delt, &
96 g, &
97 rd, &
98 rv, &
99 t0c, &
100 den0, &
101 cpd, &
102 cpv, &
103 ep1, &
104 ep2, &
105 qmin, &
106 XLS, &
107 XLV0, &
108 XLF0, &
109 cliq, &
110 cice, &
111 psat, &
112 denr
113 REAL, DIMENSION( ims:ime , jms:jme ), &
114 INTENT(INOUT) :: rain, &
115 rainncv, &
116 sr
117 ! for radiation connecting
118 INTEGER, INTENT(IN):: &
119 has_reqc, &
120 has_reqi, &
121 has_reqs
122 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), &
123 INTENT(INOUT):: &
124 re_cloud, &
125 re_ice, &
126 re_snow
127 !+---+-----------------------------------------------------------------+
128 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: & ! GT
129 refl_10cm
130 !+---+-----------------------------------------------------------------+
132 REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
133 INTENT(INOUT) :: snow, &
134 snowncv
135 ! LOCAL VAR
136 #ifdef XEON_OPTIMIZED_WSM5
137 # include "mic-wsm5-3-5-locvar.h"
138 #else
139 REAL, DIMENSION( its:ite , kts:kte ) :: t
140 REAL, DIMENSION( its:ite , kts:kte, 2 ) :: qci, qrs
141 CHARACTER*256 :: emess
142 INTEGER :: i,j,k
143 #endif
145 !+---+-----------------------------------------------------------------+
146 REAL, DIMENSION(kts:kte):: qv1d, t1d, p1d, qr1d, qs1d, dBZ
147 LOGICAL, OPTIONAL, INTENT(IN) :: diagflag
148 INTEGER, OPTIONAL, INTENT(IN) :: do_radar_ref
149 !+---+-----------------------------------------------------------------+
151 ! to calculate effective radius for radiation
152 REAL, DIMENSION( kts:kte ) :: den1d
153 REAL, DIMENSION( kts:kte ) :: qc1d
154 REAL, DIMENSION( kts:kte ) :: qi1d
155 REAL, DIMENSION( kts:kte ) :: re_qc, re_qi, re_qs
157 !+---+-----------------------------------------------------------------+
159 #ifndef XEON_OPTIMIZED_WSM5
160 DO j=jts,jte
161 DO k=kts,kte
162 DO i=its,ite
163 t(i,k)=th(i,k,j)*pii(i,k,j)
164 qci(i,k,1) = qc(i,k,j)
165 qci(i,k,2) = qi(i,k,j)
166 qrs(i,k,1) = qr(i,k,j)
167 qrs(i,k,2) = qs(i,k,j)
168 ENDDO
169 ENDDO
170 ! Sending array starting locations of optional variables may cause
171 ! troubles, so we explicitly change the call.
172 CALL wsm52D(t, q(ims,kms,j), qci, qrs &
173 ,den(ims,kms,j) &
174 ,p(ims,kms,j), delz(ims,kms,j) &
175 ,delt,g, cpd, cpv, rd, rv, t0c &
176 ,ep1, ep2, qmin &
177 ,XLS, XLV0, XLF0, den0, denr &
178 ,cliq,cice,psat &
179 ,j &
180 ,rain(ims,j),rainncv(ims,j) &
181 ,sr(ims,j) &
182 ,ids,ide, jds,jde, kds,kde &
183 ,ims,ime, jms,jme, kms,kme &
184 ,its,ite, jts,jte, kts,kte &
185 ,snow,snowncv &
186 )
187 DO K=kts,kte
188 DO I=its,ite
189 th(i,k,j)=t(i,k)/pii(i,k,j)
190 qc(i,k,j) = qci(i,k,1)
191 qi(i,k,j) = qci(i,k,2)
192 qr(i,k,j) = qrs(i,k,1)
193 qs(i,k,j) = qrs(i,k,2)
194 ENDDO
195 ENDDO
197 !+---+-----------------------------------------------------------------+
198 IF ( PRESENT (diagflag) ) THEN
199 if (diagflag .and. do_radar_ref == 1) then
200 ! WRITE(emess,*)'calling refl10cm_wsm5 ',its, jts
201 ! CALL wrf_debug ( 0, emess )
202 DO I=its,ite
203 DO K=kts,kte
204 t1d(k)=th(i,k,j)*pii(i,k,j)
205 p1d(k)=p(i,k,j)
206 qv1d(k)=q(i,k,j)
207 qr1d(k)=qr(i,k,j)
208 qs1d(k)=qs(i,k,j)
209 ENDDO
210 call refl10cm_wsm5 (qv1d, qr1d, qs1d, &
211 t1d, p1d, dBZ, kts, kte, i, j)
212 do k = kts, kte
213 refl_10cm(i,k,j) = MAX(-35., dBZ(k))
214 enddo
215 ENDDO
216 endif
217 ENDIF
219 if (has_reqc.ne.0 .and. has_reqi.ne.0 .and. has_reqs.ne.0) then
220 do i=its,ite
221 do k=kts,kte
222 re_qc(k) = RE_QC_BG
223 re_qi(k) = RE_QI_BG
224 re_qs(k) = RE_QS_BG
226 t1d(k) = th(i,k,j)*pii(i,k,j)
227 den1d(k)= den(i,k,j)
228 qc1d(k) = qc(i,k,j)
229 qi1d(k) = qi(i,k,j)
230 qs1d(k) = qs(i,k,j)
231 enddo
232 call effectRad_wsm5(t1d, qc1d, qi1d, qs1d, den1d, &
233 qmin, t0c, re_qc, re_qi, re_qs, &
234 kts, kte, i, j)
235 do k=kts,kte
236 re_cloud(i,k,j) = MAX(RE_QC_BG, MIN(re_qc(k), 50.E-6))
237 re_ice(i,k,j) = MAX(RE_QI_BG, MIN(re_qi(k), 125.E-6))
238 re_snow(i,k,j) = MAX(RE_QS_BG, MIN(re_qs(k), 999.E-6))
239 enddo
240 enddo
241 endif ! has_reqc, etc...
242 !+---+-----------------------------------------------------------------+
244 ENDDO
245 #else
246 ! Code to call the XEON-optimized WSM5 is included here
247 ! this also contains OpenMP directives, so they are
248 ! removed from the driver in the case where WSM5SCHEME
249 ! is selected.
250 # include "mic-wsm5-3-5-callsite.h"
251 IF ( PRESENT (diagflag) ) THEN
252 IF (diagflag .and. do_radar_ref == 1) then
253 !$OMP PARALLEL DO &
254 !$OMP PRIVATE ( i,j,k,t1d,p1d,qv1d,qr1d,qs1d,dBZ )
255 DO j=jts,jte
256 DO I=its,ite
257 DO K=kts,kte
258 t1d(k)=th(i,k,j)*pii(i,k,j)
259 p1d(k)=p(i,k,j)
260 qv1d(k)=q(i,k,j)
261 qr1d(k)=qr(i,k,j)
262 qs1d(k)=qs(i,k,j)
263 ENDDO
264 call refl10cm_wsm5 (qv1d, qr1d, qs1d, &
265 t1d, p1d, dBZ, kts, kte, i, j)
266 DO k = kts, kte
267 refl_10cm(i,k,j) = MAX(-35., dBZ(k))
268 ENDDO
269 ENDDO
270 ENDDO
271 !$OMP END PARALLEL DO
272 ENDIF
273 ENDIF
274 #endif
276 END SUBROUTINE wsm5
278 #ifndef XEON_OPTIMIZED_WSM5
279 !===================================================================
280 !
281 SUBROUTINE wsm52D(t, q &
282 ,qci, qrs, den, p, delz &
283 ,delt,g, cpd, cpv, rd, rv, t0c &
284 ,ep1, ep2, qmin &
285 ,XLS, XLV0, XLF0, den0, denr &
286 ,cliq,cice,psat &
287 ,lat &
288 ,rain,rainncv &
289 ,sr &
290 ,ids,ide, jds,jde, kds,kde &
291 ,ims,ime, jms,jme, kms,kme &
292 ,its,ite, jts,jte, kts,kte &
293 ,snow,snowncv &
294 )
295 !-------------------------------------------------------------------
296 IMPLICIT NONE
297 !-------------------------------------------------------------------
298 !
299 ! This code is a 5-class mixed ice microphyiscs scheme (WSM5) of the
300 ! Single-Moment MicroPhyiscs (WSMMP). The WSMMP assumes that ice nuclei
301 ! number concentration is a function of temperature, and seperate assumption
302 ! is developed, in which ice crystal number concentration is a function
303 ! of ice amount. A theoretical background of the ice-microphysics and related
304 ! processes in the WSMMPs are described in Hong et al. (2004).
305 ! Production terms in the WSM6 scheme are described in Hong and Lim (2006).
306 ! All units are in m.k.s. and source/sink terms in kgkg-1s-1.
307 !
308 ! WSM5 cloud scheme
309 !
310 ! Coded by Song-You Hong (Yonsei Univ.)
311 ! Jimy Dudhia (NCAR) and Shu-Hua Chen (UC Davis)
312 ! Summer 2002
313 !
314 ! Implemented by Song-You Hong (Yonsei Univ.) and Jimy Dudhia (NCAR)
315 ! Summer 2003
316 !
317 ! further modifications :
318 ! semi-lagrangian sedimentation (JH,2010),hong, aug 2009
319 ! ==> higher accuracy and efficient at lower resolutions
320 ! reflectivity computation from greg thompson, lim, jun 2011
321 ! ==> only diagnostic, but with removal of too large drops
322 ! effective radius of hydrometeors, bae from kiaps, jan 2015
323 ! ==> consistency in solar insolation of rrtmg radiation
324 ! bug fix in melting terms, bae from kiaps, nov 2015
325 ! ==> density of air is divided, which has not been
326 !
327 ! Reference) Hong, Dudhia, Chen (HDC, 2004) Mon. Wea. Rev.
328 ! Rutledge, Hobbs (RH83, 1983) J. Atmos. Sci.
329 ! Hong and Lim (HL, 2006) J. Korean Meteor. Soc.
330 !
331 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
332 ims,ime, jms,jme, kms,kme , &
333 its,ite, jts,jte, kts,kte, &
334 lat
335 REAL, DIMENSION( its:ite , kts:kte ), &
336 INTENT(INOUT) :: &
337 t
338 REAL, DIMENSION( its:ite , kts:kte, 2 ), &
339 INTENT(INOUT) :: &
340 qci, &
341 qrs
342 REAL, DIMENSION( ims:ime , kms:kme ), &
343 INTENT(INOUT) :: &
344 q
345 REAL, DIMENSION( ims:ime , kms:kme ), &
346 INTENT(IN ) :: &
347 den, &
348 p, &
349 delz
350 REAL, INTENT(IN ) :: delt, &
351 g, &
352 cpd, &
353 cpv, &
354 t0c, &
355 den0, &
356 rd, &
357 rv, &
358 ep1, &
359 ep2, &
360 qmin, &
361 XLS, &
362 XLV0, &
363 XLF0, &
364 cliq, &
365 cice, &
366 psat, &
367 denr
368 REAL, DIMENSION( ims:ime ), &
369 INTENT(INOUT) :: rain, &
370 rainncv, &
371 sr
372 REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, &
373 INTENT(INOUT) :: snow, &
374 snowncv
375 ! LOCAL VAR
376 REAL, DIMENSION( its:ite , kts:kte , 2) :: &
377 rh, &
378 qs, &
379 rslope, &
380 rslope2, &
381 rslope3, &
382 rslopeb, &
383 qrs_tmp, &
384 falk, &
385 fall, &
386 work1
387 REAL, DIMENSION( its:ite , kts:kte ) :: &
388 falkc, &
389 fallc, &
390 xl, &
391 cpm, &
392 denfac, &
393 xni, &
394 denqrs1, &
395 denqrs2, &
396 denqci, &
397 n0sfac, &
398 work2, &
399 workr, &
400 works, &
401 work1c, &
402 work2c
403 REAL, DIMENSION( its:ite , kts:kte ) :: &
404 den_tmp, &
405 delz_tmp
406 REAL, DIMENSION( its:ite ) :: &
407 delqrs1, &
408 delqrs2, &
409 delqi
410 REAL, DIMENSION( its:ite , kts:kte ) :: &
411 pigen, &
412 pidep, &
413 psdep, &
414 praut, &
415 psaut, &
416 prevp, &
417 psevp, &
418 pracw, &
419 psacw, &
420 psaci, &
421 pcond, &
422 psmlt
423 INTEGER, DIMENSION( its:ite ) :: &
424 mstep, &
425 numdt
426 REAL, DIMENSION(its:ite) :: tstepsnow
427 REAL, DIMENSION(its:ite) :: rmstep
428 REAL dtcldden, rdelz, rdtcld
429 LOGICAL, DIMENSION( its:ite ) :: flgcld
430 #define WSM_NO_CONDITIONAL_IN_VECTOR
431 #ifdef WSM_NO_CONDITIONAL_IN_VECTOR
432 REAL, DIMENSION(its:ite) :: xal, xbl
433 #endif
434 REAL :: &
435 cpmcal, xlcal, diffus, &
436 viscos, xka, venfac, conden, diffac, &
437 x, y, z, a, b, c, d, e, &
438 qdt, holdrr, holdrs, supcol, supcolt, pvt, &
439 coeres, supsat, dtcld, xmi, eacrs, satdt, &
440 vt2i,vt2s,acrfac, &
441 qimax, diameter, xni0, roqi0, &
442 fallsum, fallsum_qsi, xlwork2, factor, source, &
443 value, xlf, pfrzdtc, pfrzdtr, supice, holdc, holdci
444 ! variables for optimization
445 REAL, DIMENSION( its:ite ) :: tvec1
446 REAL :: temp
447 INTEGER :: i, j, k, mstepmax, &
448 iprt, latd, lond, loop, loops, ifsat, n, idim, kdim
449 ! Temporaries used for inlining fpvs function
450 REAL :: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
451 REAL :: logtr
452 !
453 !=================================================================
454 ! compute internal functions
455 !
456 cpmcal(x) = cpd*(1.-max(x,qmin))+max(x,qmin)*cpv
457 xlcal(x) = xlv0-xlv1*(x-t0c)
458 !----------------------------------------------------------------
459 ! diffus: diffusion coefficient of the water vapor
460 ! viscos: kinematic viscosity(m2s-1)
461 ! diffus(x,y) = 8.794e-5 * exp(log(x)*(1.81)) / y ! 8.794e-5*x**1.81/y
462 ! viscos(x,y) = 1.496e-6 * (x*sqrt(x)) /(x+120.)/y ! 1.496e-6*x**1.5/(x+120.)/y
463 ! xka(x,y) = 1.414e3*viscos(x,y)*y
464 ! diffac(a,b,c,d,e) = d*a*a/(xka(c,d)*rv*c*c)+1./(e*diffus(c,b))
465 ! venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333))) &
466 ! /sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
467 ! conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+d*d/(rv*e)*c/(a*a))
468 !
469 !----------------------------------------------------------------
470 idim = ite-its+1
471 kdim = kte-kts+1
472 !
473 !----------------------------------------------------------------
474 ! paddint 0 for negative values generated by dynamics
475 !
476 do k = kts, kte
477 do i = its, ite
478 qci(i,k,1) = max(qci(i,k,1),0.0)
479 qrs(i,k,1) = max(qrs(i,k,1),0.0)
480 qci(i,k,2) = max(qci(i,k,2),0.0)
481 qrs(i,k,2) = max(qrs(i,k,2),0.0)
482 enddo
483 enddo
484 !
485 !----------------------------------------------------------------
486 ! latent heat for phase changes and heat capacity. neglect the
487 ! changes during microphysical process calculation
488 ! emanuel(1994)
489 !
490 do k = kts, kte
491 do i = its, ite
492 cpm(i,k) = cpmcal(q(i,k))
493 xl(i,k) = xlcal(t(i,k))
494 enddo
495 enddo
496 do k = kts, kte
497 do i = its, ite
498 delz_tmp(i,k) = delz(i,k)
499 den_tmp(i,k) = den(i,k)
500 enddo
501 enddo
502 !
503 !----------------------------------------------------------------
504 ! initialize the surface rain, snow
505 !
506 do i = its, ite
507 rainncv(i) = 0.
508 if(PRESENT (snowncv) .AND. PRESENT (snow)) snowncv(i,lat) = 0.
509 sr(i) = 0.
510 ! new local array to catch step snow
511 tstepsnow(i) = 0.
512 enddo
513 !
514 !----------------------------------------------------------------
515 ! compute the minor time steps.
516 !
517 loops = max(nint(delt/dtcldcr),1)
518 dtcld = delt/loops
519 if(delt.le.dtcldcr) dtcld = delt
520 !
521 do loop = 1,loops
522 !
523 !----------------------------------------------------------------
524 ! initialize the large scale variables
525 !
526 do i = its, ite
527 mstep(i) = 1
528 flgcld(i) = .true.
529 enddo
530 !
531 ! do k = kts, kte
532 ! do i = its, ite
533 ! denfac(i,k) = sqrt(den0/den(i,k))
534 ! enddo
535 ! enddo
536 do k = kts, kte
537 CALL VREC( tvec1(its), den(its,k), ite-its+1)
538 do i = its, ite
539 tvec1(i) = tvec1(i)*den0
540 enddo
541 CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
542 enddo
543 !
544 ! Inline expansion for fpvs
545 ! qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
546 ! qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
547 hsub = xls
548 hvap = xlv0
549 cvap = cpv
550 ttp=t0c+0.01
551 dldt=cvap-cliq
552 xa=-dldt/rv
553 xb=xa+hvap/(rv*ttp)
554 dldti=cvap-cice
555 xai=-dldti/rv
556 xbi=xai+hsub/(rv*ttp)
557 ! this is for compilers where the conditional inhibits vectorization
558 #ifdef WSM_NO_CONDITIONAL_IN_VECTOR
559 do k = kts, kte
560 do i = its, ite
561 if(t(i,k).lt.ttp) then
562 xal(i) = xai
563 xbl(i) = xbi
564 else
565 xal(i) = xa
566 xbl(i) = xb
567 endif
568 enddo
569 do i = its, ite
570 tr=ttp/t(i,k)
571 logtr=log(tr)
572 qs(i,k,1)=psat*exp(logtr*(xa)+xb*(1.-tr))
573 qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k))
574 qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
575 qs(i,k,1) = max(qs(i,k,1),qmin)
576 rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
577 qs(i,k,2)=psat*exp(logtr*(xal(i))+xbl(i)*(1.-tr))
578 qs(i,k,2) = min(qs(i,k,2),0.99*p(i,k))
579 qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
580 qs(i,k,2) = max(qs(i,k,2),qmin)
581 rh(i,k,2) = max(q(i,k) / qs(i,k,2),qmin)
582 #else
583 do k = kts, kte
584 do i = its, ite
585 tr=ttp/t(i,k)
586 logtr=log(tr)
587 qs(i,k,1)=psat*exp(logtr*(xa)+xb*(1.-tr))
588 qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k))
589 qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
590 qs(i,k,1) = max(qs(i,k,1),qmin)
591 rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
592 if(t(i,k).lt.ttp) then
593 qs(i,k,2)=psat*exp(logtr*(xai)+xbi*(1.-tr))
594 else
595 qs(i,k,2)=psat*exp(logtr*(xa)+xb*(1.-tr))
596 endif
597 qs(i,k,2) = min(qs(i,k,2),0.99*p(i,k))
598 qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
599 qs(i,k,2) = max(qs(i,k,2),qmin)
600 rh(i,k,2) = max(q(i,k) / qs(i,k,2),qmin)
601 enddo
602 enddo
603 #endif
604 enddo
605 enddo
606 !
607 !----------------------------------------------------------------
608 ! initialize the variables for microphysical physics
609 !
610 !
611 do k = kts, kte
612 do i = its, ite
613 prevp(i,k) = 0.
614 psdep(i,k) = 0.
615 praut(i,k) = 0.
616 psaut(i,k) = 0.
617 pracw(i,k) = 0.
618 psaci(i,k) = 0.
619 psacw(i,k) = 0.
620 pigen(i,k) = 0.
621 pidep(i,k) = 0.
622 pcond(i,k) = 0.
623 psmlt(i,k) = 0.
624 psevp(i,k) = 0.
625 falk(i,k,1) = 0.
626 falk(i,k,2) = 0.
627 fall(i,k,1) = 0.
628 fall(i,k,2) = 0.
629 fallc(i,k) = 0.
630 falkc(i,k) = 0.
631 xni(i,k) = 1.e3
632 enddo
633 enddo
634 !-------------------------------------------------------------
635 ! Ni: ice crystal number concentraiton [HDC 5c]
636 !-------------------------------------------------------------
637 do k = kts, kte
638 do i = its, ite
639 temp = (den(i,k)*max(qci(i,k,2),qmin))
640 temp = sqrt(sqrt(temp*temp*temp))
641 xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
642 enddo
643 enddo
644 !
645 !----------------------------------------------------------------
646 ! compute the fallout term:
647 ! first, vertical terminal velosity for minor loops
648 !----------------------------------------------------------------
649 do k = kts, kte
650 do i = its, ite
651 qrs_tmp(i,k,1) = qrs(i,k,1)
652 qrs_tmp(i,k,2) = qrs(i,k,2)
653 enddo
654 enddo
655 call slope_wsm5(qrs_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2,rslope3, &
656 work1,its,ite,kts,kte)
657 !
658 do k = kte, kts, -1
659 do i = its, ite
660 workr(i,k) = work1(i,k,1)
661 works(i,k) = work1(i,k,2)
662 denqrs1(i,k) = den(i,k)*qrs(i,k,1)
663 denqrs2(i,k) = den(i,k)*qrs(i,k,2)
664 if(qrs(i,k,1).le.0.0) workr(i,k) = 0.0
665 if(qrs(i,k,2).le.0.0) works(i,k) = 0.0
666 enddo
667 enddo
668 call nislfv_rain_plm(idim,kdim,den_tmp,denfac,t,delz_tmp,workr,denqrs1, &
669 delqrs1,dtcld,1,1)
670 call nislfv_rain_plm(idim,kdim,den_tmp,denfac,t,delz_tmp,works,denqrs2, &
671 delqrs2,dtcld,2,1)
672 do k = kts, kte
673 do i = its, ite
674 qrs(i,k,1) = max(denqrs1(i,k)/den(i,k),0.)
675 qrs(i,k,2) = max(denqrs2(i,k)/den(i,k),0.)
676 fall(i,k,1) = denqrs1(i,k)*workr(i,k)/delz(i,k)
677 fall(i,k,2) = denqrs2(i,k)*works(i,k)/delz(i,k)
678 enddo
679 enddo
680 do i = its, ite
681 fall(i,1,1) = delqrs1(i)/delz(i,1)/dtcld
682 fall(i,1,2) = delqrs2(i)/delz(i,1)/dtcld
683 enddo
684 do k = kts, kte
685 do i = its, ite
686 qrs_tmp(i,k,1) = qrs(i,k,1)
687 qrs_tmp(i,k,2) = qrs(i,k,2)
688 enddo
689 enddo
690 call slope_wsm5(qrs_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2,rslope3, &
691 work1,its,ite,kts,kte)
692 do k = kte, kts, -1
693 do i = its, ite
694 supcol = t0c-t(i,k)
695 n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
696 if(t(i,k).gt.t0c.and.qrs(i,k,2).gt.0.) then
697 !----------------------------------------------------------------
698 ! psmlt: melting of snow [HL A33] [RH83 A25]
699 ! (T>T0: S->R)
700 !----------------------------------------------------------------
701 xlf = xlf0
702 ! work2(i,k)= venfac(p(i,k),t(i,k),den(i,k))
703 work2(i,k)= (exp(log(((1.496e-6*((t(i,k))*sqrt(t(i,k))) &
704 /((t(i,k))+120.)/(den(i,k)))/(8.794e-5 &
705 *exp(log(t(i,k))*(1.81))/p(i,k)))) &
706 *((.3333333)))/sqrt((1.496e-6*((t(i,k)) &
707 *sqrt(t(i,k)))/((t(i,k))+120.)/(den(i,k)))) &
708 *sqrt(sqrt(den0/(den(i,k)))))
709 coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
710 psmlt(i,k) = (1.414e3*(1.496e-6*((t(i,k))*sqrt(t(i,k))) &
711 /((t(i,k))+120.)/(den(i,k)) )*(den(i,k))) &
712 /xlf*(t0c-t(i,k))*pi/2. &
713 *n0sfac(i,k)*(precs1*rslope2(i,k,2)+precs2 &
714 *work2(i,k)*coeres)/den(i,k)
715 psmlt(i,k) = min(max(psmlt(i,k)*dtcld/mstep(i), &
716 -qrs(i,k,2)/mstep(i)),0.)
717 qrs(i,k,2) = qrs(i,k,2) + psmlt(i,k)
718 qrs(i,k,1) = qrs(i,k,1) - psmlt(i,k)
719 t(i,k) = t(i,k) + xlf/cpm(i,k)*psmlt(i,k)
720 endif
721 enddo
722 enddo
723 !---------------------------------------------------------------
724 ! Vice [ms-1] : fallout of ice crystal [HDC 5a]
725 !---------------------------------------------------------------
726 do k = kte, kts, -1
727 do i = its, ite
728 if(qci(i,k,2).le.0.) then
729 work1c(i,k) = 0.
730 else
731 xmi = den(i,k)*qci(i,k,2)/xni(i,k)
732 diameter = max(min(dicon * sqrt(xmi),dimax), 1.e-25)
733 work1c(i,k) = 1.49e4*exp(log(diameter)*(1.31))
734 endif
735 enddo
736 enddo
737 !
738 ! forward semi-laglangian scheme (JH), PCM (piecewise constant), (linear)
739 !
740 do k = kte, kts, -1
741 do i = its, ite
742 denqci(i,k) = den(i,k)*qci(i,k,2)
743 enddo
744 enddo
745 call nislfv_rain_plm(idim,kdim,den_tmp,denfac,t,delz_tmp,work1c,denqci, &
746 delqi,dtcld,1,0)
747 do k = kts, kte
748 do i = its, ite
749 qci(i,k,2) = max(denqci(i,k)/den(i,k),0.)
750 enddo
751 enddo
752 do i = its, ite
753 fallc(i,1) = delqi(i)/delz(i,1)/dtcld
754 enddo
755 !
756 !----------------------------------------------------------------
757 ! rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
758 !
759 do i = its, ite
760 fallsum = fall(i,1,1)+fall(i,1,2)+fallc(i,1)
761 fallsum_qsi = fall(i,1,2)+fallc(i,1)
762 if(fallsum.gt.0.) then
763 rainncv(i) = fallsum*delz(i,1)/denr*dtcld*1000. + rainncv(i)
764 rain(i) = fallsum*delz(i,1)/denr*dtcld*1000. + rain(i)
765 endif
766 if(fallsum_qsi.gt.0.) then
767 tstepsnow(i) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. &
768 +tstepsnow(i)
769 IF ( PRESENT (snowncv) .AND. PRESENT (snow)) THEN
770 snowncv(i,lat) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. &
771 +snowncv(i,lat)
772 snow(i,lat) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. + snow(i,lat)
773 ENDIF
774 endif
775 IF ( PRESENT (snowncv) ) THEN
776 if(fallsum.gt.0.)sr(i)=snowncv(i,lat)/(rainncv(i)+1.e-12)
777 ELSE
778 if(fallsum.gt.0.)sr(i)=tstepsnow(i)/(rainncv(i)+1.e-12)
779 ENDIF
780 enddo
781 !
782 !---------------------------------------------------------------
783 ! pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
784 ! (T>T0: I->C)
785 !---------------------------------------------------------------
786 do k = kts, kte
787 do i = its, ite
788 supcol = t0c-t(i,k)
789 xlf = xls-xl(i,k)
790 if(supcol.lt.0.) xlf = xlf0
791 if(supcol.lt.0.and.qci(i,k,2).gt.0.) then
792 qci(i,k,1) = qci(i,k,1) + qci(i,k,2)
793 t(i,k) = t(i,k) - xlf/cpm(i,k)*qci(i,k,2)
794 qci(i,k,2) = 0.
795 endif
796 !---------------------------------------------------------------
797 ! pihmf: homogeneous freezing of cloud water below -40c [HL A45]
798 ! (T<-40C: C->I)
799 !---------------------------------------------------------------
800 if(supcol.gt.40..and.qci(i,k,1).gt.0.) then
801 qci(i,k,2) = qci(i,k,2) + qci(i,k,1)
802 t(i,k) = t(i,k) + xlf/cpm(i,k)*qci(i,k,1)
803 qci(i,k,1) = 0.
804 endif
805 !---------------------------------------------------------------
806 ! pihtf: heterogeneous freezing of cloud water [HL A44]
807 ! (T0>T>-40C: C->I)
808 !---------------------------------------------------------------
809 if(supcol.gt.0..and.qci(i,k,1).gt.0.) then
810 supcolt=min(supcol,50.)
811 ! pfrzdtc = min(pfrz1*(exp(pfrz2*supcol)-1.) &
812 ! *den(i,k)/denr/xncr*qci(i,k,1)**2*dtcld,qci(i,k,1))
813 pfrzdtc = min(pfrz1*(exp(pfrz2*supcolt)-1.) &
814 *den(i,k)/denr/xncr*qci(i,k,1)*qci(i,k,1)*dtcld,qci(i,k,1))
815 qci(i,k,2) = qci(i,k,2) + pfrzdtc
816 t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtc
817 qci(i,k,1) = qci(i,k,1)-pfrzdtc
818 endif
819 !---------------------------------------------------------------
820 ! psfrz: freezing of rain water [HL A20] [LFO 45]
821 ! (T<T0, R->S)
822 !---------------------------------------------------------------
823 if(supcol.gt.0..and.qrs(i,k,1).gt.0.) then
824 supcolt=min(supcol,50.)
825 ! pfrzdtr = min(20.*pi**2*pfrz1*n0r*denr/den(i,k) &
826 ! *(exp(pfrz2*supcol)-1.)*rslope(i,k,1)**7*dtcld, &
827 ! qrs(i,k,1))
828 temp = rslope(i,k,1)
829 temp = temp*temp*temp*temp*temp*temp*temp
830 pfrzdtr = min(20.*(pi*pi)*pfrz1*n0r*denr/den(i,k) &
831 *(exp(pfrz2*supcolt)-1.)*temp*dtcld, &
832 qrs(i,k,1))
833 qrs(i,k,2) = qrs(i,k,2) + pfrzdtr
834 t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtr
835 qrs(i,k,1) = qrs(i,k,1)-pfrzdtr
836 endif
837 enddo
838 enddo
839 !
840 !----------------------------------------------------------------
841 ! update the slope parameters for microphysics computation
842 !
843 do k = kts, kte
844 do i = its, ite
845 qrs_tmp(i,k,1) = qrs(i,k,1)
846 qrs_tmp(i,k,2) = qrs(i,k,2)
847 enddo
848 enddo
849 call slope_wsm5(qrs_tmp,den_tmp,denfac,t,rslope,rslopeb,rslope2,rslope3, &
850 work1,its,ite,kts,kte)
851 !------------------------------------------------------------------
852 ! work1: the thermodynamic term in the denominator associated with
853 ! heat conduction and vapor diffusion
854 ! (ry88, y93, h85)
855 ! work2: parameter associated with the ventilation effects(y93)
856 !
857 do k = kts, kte
858 do i = its, ite
859 ! work1(i,k,1) = diffac(xl(i,k),p(i,k),t(i,k),den(i,k),qs(i,k,1))
860 work1(i,k,1) = ((((den(i,k))*(xl(i,k))*(xl(i,k)))*((t(i,k))+120.) &
861 *(den(i,k)))/(1.414e3*(1.496e-6*((t(i,k))*sqrt(t(i,k))))&
862 *(den(i,k))*(rv*(t(i,k))*(t(i,k))))) &
863 + p(i,k)/((qs(i,k,1))*(8.794e-5*exp(log(t(i,k))*(1.81))))
864 ! work1(i,k,2) = diffac(xls,p(i,k),t(i,k),den(i,k),qs(i,k,2))
865 work1(i,k,2) = ((((den(i,k))*(xls)*(xls))*((t(i,k))+120.)*(den(i,k)))&
866 /(1.414e3*(1.496e-6*((t(i,k))*sqrt(t(i,k))))*(den(i,k)) &
867 *(rv*(t(i,k))*(t(i,k)))) &
868 + p(i,k)/(qs(i,k,2)*(8.794e-5*exp(log(t(i,k))*(1.81)))))
869 ! work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
870 work2(i,k) = (exp(.3333333*log(((1.496e-6 * ((t(i,k))*sqrt(t(i,k)))) &
871 *p(i,k))/(((t(i,k))+120.)*den(i,k)*(8.794e-5 &
872 *exp(log(t(i,k))*(1.81))))))*sqrt(sqrt(den0/(den(i,k))))) &
873 /sqrt((1.496e-6*((t(i,k))*sqrt(t(i,k)))) &
874 /(((t(i,k))+120.)*den(i,k)))
875 enddo
876 enddo
877 !
878 !===============================================================
879 !
880 ! warm rain processes
881 !
882 ! - follows the processes in RH83 and LFO except for autoconcersion
883 !
884 !===============================================================
885 !
886 do k = kts, kte
887 do i = its, ite
888 supsat = max(q(i,k),qmin)-qs(i,k,1)
889 satdt = supsat/dtcld
890 !---------------------------------------------------------------
891 ! praut: auto conversion rate from cloud to rain [HDC 16]
892 ! (C->R)
893 !---------------------------------------------------------------
894 if(qci(i,k,1).gt.qc0) then
895 praut(i,k) = qck1*exp(log(qci(i,k,1))*((7./3.)))
896 praut(i,k) = min(praut(i,k),qci(i,k,1)/dtcld)
897 endif
898 !---------------------------------------------------------------
899 ! pracw: accretion of cloud water by rain [HL A40] [LFO 51]
900 ! (C->R)
901 !---------------------------------------------------------------
902 if(qrs(i,k,1).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
903 pracw(i,k) = min(pacrr*rslope3(i,k,1)*rslopeb(i,k,1) &
904 *qci(i,k,1)*denfac(i,k),qci(i,k,1)/dtcld)
905 endif
906 !---------------------------------------------------------------
907 ! prevp: evaporation/condensation rate of rain [HDC 14]
908 ! (V->R or R->V)
909 !---------------------------------------------------------------
910 if(qrs(i,k,1).gt.0.) then
911 coeres = rslope2(i,k,1)*sqrt(rslope(i,k,1)*rslopeb(i,k,1))
912 prevp(i,k) = (rh(i,k,1)-1.)*(precr1*rslope2(i,k,1) &
913 +precr2*work2(i,k)*coeres)/work1(i,k,1)
914 if(prevp(i,k).lt.0.) then
915 prevp(i,k) = max(prevp(i,k),-qrs(i,k,1)/dtcld)
916 prevp(i,k) = max(prevp(i,k),satdt/2)
917 else
918 prevp(i,k) = min(prevp(i,k),satdt/2)
919 endif
920 endif
921 enddo
922 enddo
923 !
924 !===============================================================
925 !
926 ! cold rain processes
927 !
928 ! - follows the revised ice microphysics processes in HDC
929 ! - the processes same as in RH83 and RH84 and LFO behave
930 ! following ice crystal hapits defined in HDC, inclduing
931 ! intercept parameter for snow (n0s), ice crystal number
932 ! concentration (ni), ice nuclei number concentration
933 ! (n0i), ice diameter (d)
934 !
935 !===============================================================
936 !
937 rdtcld = 1./dtcld
938 do k = kts, kte
939 do i = its, ite
940 supcol = t0c-t(i,k)
941 n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
942 supsat = max(q(i,k),qmin)-qs(i,k,2)
943 satdt = supsat/dtcld
944 ifsat = 0
945 !-------------------------------------------------------------
946 ! Ni: ice crystal number concentraiton [HDC 5c]
947 !-------------------------------------------------------------
948 ! xni(i,k) = min(max(5.38e7*(den(i,k) &
949 ! *max(qci(i,k,2),qmin))**0.75,1.e3),1.e6)
950 temp = (den(i,k)*max(qci(i,k,2),qmin))
951 temp = sqrt(sqrt(temp*temp*temp))
952 xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
953 eacrs = exp(0.07*(-supcol))
954 !
955 if(supcol.gt.0) then
956 if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,2).gt.qmin) then
957 xmi = den(i,k)*qci(i,k,2)/xni(i,k)
958 diameter = min(dicon * sqrt(xmi),dimax)
959 vt2i = 1.49e4*diameter**1.31
960 vt2s = pvts*rslopeb(i,k,2)*denfac(i,k)
961 !-------------------------------------------------------------
962 ! psaci: Accretion of cloud ice by rain [HDC 10]
963 ! (T<T0: I->S)
964 !-------------------------------------------------------------
965 acrfac = 2.*rslope3(i,k,2)+2.*diameter*rslope2(i,k,2) &
966 +diameter**2*rslope(i,k,2)
967 psaci(i,k) = pi*qci(i,k,2)*eacrs*n0s*n0sfac(i,k) &
968 *abs(vt2s-vt2i)*acrfac/4.
969 endif
970 endif
971 !-------------------------------------------------------------
972 ! psacw: Accretion of cloud water by snow [HL A7] [LFO 24]
973 ! (T<T0: C->S, and T>=T0: C->R)
974 !-------------------------------------------------------------
975 if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
976 psacw(i,k) = min(pacrc*n0sfac(i,k)*rslope3(i,k,2) &
977 *rslopeb(i,k,2)*qci(i,k,1)*denfac(i,k) &
978 ! ,qci(i,k,1)/dtcld)
979 ,qci(i,k,1)*rdtcld)
980 endif
981 if(supcol .gt. 0) then
982 !-------------------------------------------------------------
983 ! pidep: Deposition/Sublimation rate of ice [HDC 9]
984 ! (T<T0: V->I or I->V)
985 !-------------------------------------------------------------
986 if(qci(i,k,2).gt.0.and.ifsat.ne.1) then
987 xmi = den(i,k)*qci(i,k,2)/xni(i,k)
988 diameter = dicon * sqrt(xmi)
989 pidep(i,k) = 4.*diameter*xni(i,k)*(rh(i,k,2)-1.)/work1(i,k,2)
990 supice = satdt-prevp(i,k)
991 if(pidep(i,k).lt.0.) then
992 ! pidep(i,k) = max(max(pidep(i,k),satdt/2),supice)
993 ! pidep(i,k) = max(pidep(i,k),-qci(i,k,2)/dtcld)
994 pidep(i,k) = max(max(pidep(i,k),satdt*.5),supice)
995 pidep(i,k) = max(pidep(i,k),-qci(i,k,2)*rdtcld)
996 else
997 ! pidep(i,k) = min(min(pidep(i,k),satdt/2),supice)
998 pidep(i,k) = min(min(pidep(i,k),satdt*.5),supice)
999 endif
1000 if(abs(prevp(i,k)+pidep(i,k)).ge.abs(satdt)) ifsat = 1
1001 endif
1002 !-------------------------------------------------------------
1003 ! psdep: deposition/sublimation rate of snow [HDC 14]
1004 ! (V->S or S->V)
1005 !-------------------------------------------------------------
1006 if(qrs(i,k,2).gt.0..and.ifsat.ne.1) then
1007 coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
1008 psdep(i,k) = (rh(i,k,2)-1.)*n0sfac(i,k) &
1009 *(precs1*rslope2(i,k,2)+precs2 &
1010 *work2(i,k)*coeres)/work1(i,k,2)
1011 supice = satdt-prevp(i,k)-pidep(i,k)
1012 if(psdep(i,k).lt.0.) then
1013 ! psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)/dtcld)
1014 ! psdep(i,k) = max(max(psdep(i,k),satdt/2),supice)
1015 psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)*rdtcld)
1016 psdep(i,k) = max(max(psdep(i,k),satdt*.5),supice)
1017 else
1018 ! psdep(i,k) = min(min(psdep(i,k),satdt/2),supice)
1019 psdep(i,k) = min(min(psdep(i,k),satdt*.5),supice)
1020 endif
1021 if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)).ge.abs(satdt)) &
1022 ifsat = 1
1023 endif
1024 !-------------------------------------------------------------
1025 ! pigen: generation(nucleation) of ice from vapor [HL A50] [HDC 7-8]
1026 ! (T<T0: V->I)
1027 !-------------------------------------------------------------
1028 if(supsat.gt.0.and.ifsat.ne.1) then
1029 supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)
1030 xni0 = 1.e3*exp(0.1*supcol)
1031 roqi0 = 4.92e-11*exp(log(xni0)*(1.33))
1032 pigen(i,k) = max(0.,(roqi0/den(i,k)-max(qci(i,k,2),0.)) &
1033 ! /dtcld)
1034 *rdtcld)
1035 pigen(i,k) = min(min(pigen(i,k),satdt),supice)
1036 endif
1037 !
1038 !-------------------------------------------------------------
1039 ! psaut: conversion(aggregation) of ice to snow [HDC 12]
1040 ! (T<T0: I->S)
1041 !-------------------------------------------------------------
1042 if(qci(i,k,2).gt.0.) then
1043 qimax = roqimax/den(i,k)
1044 ! psaut(i,k) = max(0.,(qci(i,k,2)-qimax)/dtcld)
1045 psaut(i,k) = max(0.,(qci(i,k,2)-qimax)*rdtcld)
1046 endif
1047 endif
1048 !-------------------------------------------------------------
1049 ! psevp: Evaporation of melting snow [HL A35] [RH83 A27]
1050 ! (T>T0: S->V)
1051 !-------------------------------------------------------------
1052 if(supcol.lt.0.) then
1053 if(qrs(i,k,2).gt.0..and.rh(i,k,1).lt.1.) &
1054 psevp(i,k) = psdep(i,k)*work1(i,k,2)/work1(i,k,1)
1055 ! psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)/dtcld),0.)
1056 psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)*rdtcld),0.)
1057 endif
1058 enddo
1059 enddo
1060 !
1061 !
1062 !----------------------------------------------------------------
1063 ! check mass conservation of generation terms and feedback to the
1064 ! large scale
1065 !
1066 do k = kts, kte
1067 do i = its, ite
1068 if(t(i,k).le.t0c) then
1069 !
1070 ! cloud water
1071 !
1072 value = max(qmin,qci(i,k,1))
1073 source = (praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
1074 if (source.gt.value) then
1075 factor = value/source
1076 praut(i,k) = praut(i,k)*factor
1077 pracw(i,k) = pracw(i,k)*factor
1078 psacw(i,k) = psacw(i,k)*factor
1079 endif
1080 !
1081 ! cloud ice
1082 !
1083 value = max(qmin,qci(i,k,2))
1084 source = (psaut(i,k)+psaci(i,k)-pigen(i,k)-pidep(i,k))*dtcld
1085 if (source.gt.value) then
1086 factor = value/source
1087 psaut(i,k) = psaut(i,k)*factor
1088 psaci(i,k) = psaci(i,k)*factor
1089 pigen(i,k) = pigen(i,k)*factor
1090 pidep(i,k) = pidep(i,k)*factor
1091 endif
1092 !
1093 ! rain
1094 !
1095 !
1096 value = max(qmin,qrs(i,k,1))
1097 source = (-praut(i,k)-pracw(i,k)-prevp(i,k))*dtcld
1098 if (source.gt.value) then
1099 factor = value/source
1100 praut(i,k) = praut(i,k)*factor
1101 pracw(i,k) = pracw(i,k)*factor
1102 prevp(i,k) = prevp(i,k)*factor
1103 endif
1104 !
1105 ! snow
1106 !
1107 value = max(qmin,qrs(i,k,2))
1108 source = (-psdep(i,k)-psaut(i,k)-psaci(i,k)-psacw(i,k))*dtcld
1109 if (source.gt.value) then
1110 factor = value/source
1111 psdep(i,k) = psdep(i,k)*factor
1112 psaut(i,k) = psaut(i,k)*factor
1113 psaci(i,k) = psaci(i,k)*factor
1114 psacw(i,k) = psacw(i,k)*factor
1115 endif
1116 !
1117 work2(i,k)=-(prevp(i,k)+psdep(i,k)+pigen(i,k)+pidep(i,k))
1118 ! update
1119 q(i,k) = q(i,k)+work2(i,k)*dtcld
1120 qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
1121 +psacw(i,k))*dtcld,0.)
1122 qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
1123 +prevp(i,k))*dtcld,0.)
1124 qci(i,k,2) = max(qci(i,k,2)-(psaut(i,k)+psaci(i,k) &
1125 -pigen(i,k)-pidep(i,k))*dtcld,0.)
1126 qrs(i,k,2) = max(qrs(i,k,2)+(psdep(i,k)+psaut(i,k) &
1127 +psaci(i,k)+psacw(i,k))*dtcld,0.)
1128 xlf = xls-xl(i,k)
1129 xlwork2 = -xls*(psdep(i,k)+pidep(i,k)+pigen(i,k)) &
1130 -xl(i,k)*prevp(i,k)-xlf*psacw(i,k)
1131 t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
1132 else
1133 !
1134 ! cloud water
1135 !
1136 value = max(qmin,qci(i,k,1))
1137 source=(praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
1138 if (source.gt.value) then
1139 factor = value/source
1140 praut(i,k) = praut(i,k)*factor
1141 pracw(i,k) = pracw(i,k)*factor
1142 psacw(i,k) = psacw(i,k)*factor
1143 endif
1144 !
1145 ! rain
1146 !
1147 value = max(qmin,qrs(i,k,1))
1148 source = (-praut(i,k)-pracw(i,k)-prevp(i,k)-psacw(i,k))*dtcld
1149 if (source.gt.value) then
1150 factor = value/source
1151 praut(i,k) = praut(i,k)*factor
1152 pracw(i,k) = pracw(i,k)*factor
1153 prevp(i,k) = prevp(i,k)*factor
1154 psacw(i,k) = psacw(i,k)*factor
1155 endif
1156 !
1157 ! snow
1158 !
1159 value = max(qcrmin,qrs(i,k,2))
1160 source=(-psevp(i,k))*dtcld
1161 if (source.gt.value) then
1162 factor = value/source
1163 psevp(i,k) = psevp(i,k)*factor
1164 endif
1165 work2(i,k)=-(prevp(i,k)+psevp(i,k))
1166 ! update
1167 q(i,k) = q(i,k)+work2(i,k)*dtcld
1168 qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
1169 +psacw(i,k))*dtcld,0.)
1170 qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
1171 +prevp(i,k) +psacw(i,k))*dtcld,0.)
1172 qrs(i,k,2) = max(qrs(i,k,2)+psevp(i,k)*dtcld,0.)
1173 xlf = xls-xl(i,k)
1174 xlwork2 = -xl(i,k)*(prevp(i,k)+psevp(i,k))
1175 t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
1176 endif
1177 enddo
1178 enddo
1179 !
1180 ! Inline expansion for fpvs
1181 ! qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1182 ! qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1183 hsub = xls
1184 hvap = xlv0
1185 cvap = cpv
1186 ttp=t0c+0.01
1187 dldt=cvap-cliq
1188 xa=-dldt/rv
1189 xb=xa+hvap/(rv*ttp)
1190 dldti=cvap-cice
1191 xai=-dldti/rv
1192 xbi=xai+hsub/(rv*ttp)
1193 do k = kts, kte
1194 do i = its, ite
1195 tr=ttp/t(i,k)
1196 logtr = log(tr)
1197 qs(i,k,1)=psat*exp(logtr*(xa)+xb*(1.-tr))
1198 qs(i,k,1) = min(qs(i,k,1),0.99*p(i,k))
1199 qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
1200 qs(i,k,1) = max(qs(i,k,1),qmin)
1201 enddo
1202 enddo
1203 !
1204 !----------------------------------------------------------------
1205 ! pcond: condensational/evaporational rate of cloud water [HL A46] [RH83 A6]
1206 ! if there exists additional water vapor condensated/if
1207 ! evaporation of cloud water is not enough to remove subsaturation
1208 !
1209 do k = kts, kte
1210 do i = its, ite
1211 ! work1(i,k,1) = conden(t(i,k),q(i,k),qs(i,k,1),xl(i,k),cpm(i,k))
1212 work1(i,k,1) = ((max(q(i,k),qmin)-(qs(i,k,1)))/(1.+(xl(i,k)) &
1213 *(xl(i,k))/(rv*(cpm(i,k)))*(qs(i,k,1)) &
1214 /((t(i,k))*(t(i,k)))))
1215 work2(i,k) = qci(i,k,1)+work1(i,k,1)
1216 pcond(i,k) = min(max(work1(i,k,1)/dtcld,0.),max(q(i,k),0.)/dtcld)
1217 if(qci(i,k,1).gt.0..and.work1(i,k,1).lt.0.) &
1218 pcond(i,k) = max(work1(i,k,1),-qci(i,k,1))/dtcld
1219 q(i,k) = q(i,k)-pcond(i,k)*dtcld
1220 qci(i,k,1) = max(qci(i,k,1)+pcond(i,k)*dtcld,0.)
1221 t(i,k) = t(i,k)+pcond(i,k)*xl(i,k)/cpm(i,k)*dtcld
1222 enddo
1223 enddo
1224 !
1225 !
1226 !----------------------------------------------------------------
1227 ! padding for small values
1228 !
1229 do k = kts, kte
1230 do i = its, ite
1231 if(qci(i,k,1).le.qmin) qci(i,k,1) = 0.0
1232 if(qci(i,k,2).le.qmin) qci(i,k,2) = 0.0
1233 enddo
1234 enddo
1235 enddo ! big loops
1236 END SUBROUTINE wsm52d
1237 ! ...................................................................
1238 !--------------------------------------------------------------------------
1239 REAL FUNCTION fpvs(t,ice,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c)
1240 !--------------------------------------------------------------------------
1241 IMPLICIT NONE
1242 !--------------------------------------------------------------------------
1243 REAL t,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c,dldt,xa,xb,dldti, &
1244 xai,xbi,ttp,tr
1245 INTEGER ice
1246 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1247 ttp=t0c+0.01
1248 dldt=cvap-cliq
1249 xa=-dldt/rv
1250 xb=xa+hvap/(rv*ttp)
1251 dldti=cvap-cice
1252 xai=-dldti/rv
1253 xbi=xai+hsub/(rv*ttp)
1254 tr=ttp/t
1255 if(t.lt.ttp.and.ice.eq.1) then
1256 fpvs=psat*exp(log(tr)*(xai))*exp(xbi*(1.-tr))
1257 else
1258 fpvs=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
1259 endif
1260 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1261 END FUNCTION fpvs
1262 !------------------------------------------------------------------------------
1263 subroutine slope_wsm5(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3, &
1264 vt,its,ite,kts,kte)
1265 IMPLICIT NONE
1266 INTEGER :: its,ite, jts,jte, kts,kte
1267 REAL, DIMENSION( its:ite , kts:kte,2) :: &
1268 qrs, &
1269 rslope, &
1270 rslopeb, &
1271 rslope2, &
1272 rslope3, &
1273 vt
1274 REAL, DIMENSION( its:ite , kts:kte) :: &
1275 den, &
1276 denfac, &
1277 t
1278 REAL, PARAMETER :: t0c = 273.15
1279 REAL, DIMENSION( its:ite , kts:kte ) :: &
1280 n0sfac
1281 REAL :: lamdar, lamdas, x, y, z, supcol
1282 integer :: i, j, k
1283 !----------------------------------------------------------------
1284 ! size distributions: (x=mixing ratio, y=air density):
1285 ! valid for mixing ratio > 1.e-9 kg/kg.
1286 lamdar(x,y)= sqrt(sqrt(pidn0r/(x*y))) ! (pidn0r/(x*y))**.25
1287 lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y))) ! (pidn0s*z/(x*y))**.25
1288 !
1289 do k = kts, kte
1290 do i = its, ite
1291 supcol = t0c-t(i,k)
1292 !---------------------------------------------------------------
1293 ! n0s: Intercept parameter for snow [m-4] [HDC 6]
1294 !---------------------------------------------------------------
1295 n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
1296 if(qrs(i,k,1).le.qcrmin)then
1297 rslope(i,k,1) = rslopermax
1298 rslopeb(i,k,1) = rsloperbmax
1299 rslope2(i,k,1) = rsloper2max
1300 rslope3(i,k,1) = rsloper3max
1301 else
1302 rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k))
1303 rslopeb(i,k,1) = exp(log(rslope(i,k,1))*(bvtr))
1304 rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
1305 rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
1306 endif
1307 if(qrs(i,k,2).le.qcrmin)then
1308 rslope(i,k,2) = rslopesmax
1309 rslopeb(i,k,2) = rslopesbmax
1310 rslope2(i,k,2) = rslopes2max
1311 rslope3(i,k,2) = rslopes3max
1312 else
1313 rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k))
1314 rslopeb(i,k,2) = exp(log(rslope(i,k,2))*(bvts))
1315 rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
1316 rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
1317 endif
1318 vt(i,k,1) = pvtr*rslopeb(i,k,1)*denfac(i,k)
1319 vt(i,k,2) = pvts*rslopeb(i,k,2)*denfac(i,k)
1320 if(qrs(i,k,1).le.0.0) vt(i,k,1) = 0.0
1321 if(qrs(i,k,2).le.0.0) vt(i,k,2) = 0.0
1322 enddo
1323 enddo
1324 END subroutine slope_wsm5
1325 !-----------------------------------------------------------------------------
1326 subroutine slope_rain(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3, &
1327 vt,its,ite,kts,kte)
1328 IMPLICIT NONE
1329 INTEGER :: its,ite, jts,jte, kts,kte
1330 REAL, DIMENSION( its:ite , kts:kte) :: &
1331 qrs, &
1332 rslope, &
1333 rslopeb, &
1334 rslope2, &
1335 rslope3, &
1336 vt, &
1337 den, &
1338 denfac, &
1339 t
1340 REAL, PARAMETER :: t0c = 273.15
1341 REAL, DIMENSION( its:ite , kts:kte ) :: &
1342 n0sfac
1343 REAL :: lamdar, x, y, z, supcol
1344 integer :: i, j, k
1345 !----------------------------------------------------------------
1346 ! size distributions: (x=mixing ratio, y=air density):
1347 ! valid for mixing ratio > 1.e-9 kg/kg.
1348 lamdar(x,y)= sqrt(sqrt(pidn0r/(x*y))) ! (pidn0r/(x*y))**.25
1349 !
1350 do k = kts, kte
1351 do i = its, ite
1352 if(qrs(i,k).le.qcrmin)then
1353 rslope(i,k) = rslopermax
1354 rslopeb(i,k) = rsloperbmax
1355 rslope2(i,k) = rsloper2max
1356 rslope3(i,k) = rsloper3max
1357 else
1358 rslope(i,k) = 1./lamdar(qrs(i,k),den(i,k))
1359 rslopeb(i,k) = rslope(i,k)**bvtr
1360 rslope2(i,k) = rslope(i,k)*rslope(i,k)
1361 rslope3(i,k) = rslope2(i,k)*rslope(i,k)
1362 endif
1363 vt(i,k) = pvtr*rslopeb(i,k)*denfac(i,k)
1364 if(qrs(i,k).le.0.0) vt(i,k) = 0.0
1365 enddo
1366 enddo
1367 END subroutine slope_rain
1368 !------------------------------------------------------------------------------
1369 subroutine slope_snow(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3, &
1370 vt,its,ite,kts,kte)
1371 IMPLICIT NONE
1372 INTEGER :: its,ite, jts,jte, kts,kte
1373 REAL, DIMENSION( its:ite , kts:kte) :: &
1374 qrs, &
1375 rslope, &
1376 rslopeb, &
1377 rslope2, &
1378 rslope3, &
1379 vt, &
1380 den, &
1381 denfac, &
1382 t
1383 REAL, PARAMETER :: t0c = 273.15
1384 REAL, DIMENSION( its:ite , kts:kte ) :: &
1385 n0sfac
1386 REAL :: lamdas, x, y, z, supcol
1387 integer :: i, j, k
1388 !----------------------------------------------------------------
1389 ! size distributions: (x=mixing ratio, y=air density):
1390 ! valid for mixing ratio > 1.e-9 kg/kg.
1391 lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y))) ! (pidn0s*z/(x*y))**.25
1392 !
1393 do k = kts, kte
1394 do i = its, ite
1395 supcol = t0c-t(i,k)
1396 !---------------------------------------------------------------
1397 ! n0s: Intercept parameter for snow [m-4] [HDC 6]
1398 !---------------------------------------------------------------
1399 n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
1400 if(qrs(i,k).le.qcrmin)then
1401 rslope(i,k) = rslopesmax
1402 rslopeb(i,k) = rslopesbmax
1403 rslope2(i,k) = rslopes2max
1404 rslope3(i,k) = rslopes3max
1405 else
1406 rslope(i,k) = 1./lamdas(qrs(i,k),den(i,k),n0sfac(i,k))
1407 rslopeb(i,k) = rslope(i,k)**bvts
1408 rslope2(i,k) = rslope(i,k)*rslope(i,k)
1409 rslope3(i,k) = rslope2(i,k)*rslope(i,k)
1410 endif
1411 vt(i,k) = pvts*rslopeb(i,k)*denfac(i,k)
1412 if(qrs(i,k).le.0.0) vt(i,k) = 0.0
1413 enddo
1414 enddo
1415 END subroutine slope_snow
1416 !-------------------------------------------------------------------
1417 SUBROUTINE nislfv_rain_plm(im,km,denl,denfacl,tkl,dzl,wwl,rql,precip,dt,id,iter)
1418 !-------------------------------------------------------------------
1419 !
1420 ! for non-iteration semi-Lagrangain forward advection for cloud
1421 ! with mass conservation and positive definite advection
1422 ! 2nd order interpolation with monotonic piecewise linear method
1423 ! this routine is under assumption of decfl < 1 for semi_Lagrangian
1424 !
1425 ! dzl depth of model layer in meter
1426 ! wwl terminal velocity at model layer m/s
1427 ! rql cloud density*mixing ration
1428 ! precip precipitation
1429 ! dt time step
1430 ! id kind of precip: 0 test case; 1 raindrop 2: snow
1431 ! iter how many time to guess mean terminal velocity: 0 pure forward.
1432 ! 0 : use departure wind for advection
1433 ! 1 : use mean wind for advection
1434 ! > 1 : use mean wind after iter-1 iterations
1435 !
1436 ! author: hann-ming henry juang <henry.juang@noaa.gov>
1437 ! implemented by song-you hong
1438 !
1439 implicit none
1440 integer im,km,id
1441 real dt
1442 real dzl(im,km),wwl(im,km),rql(im,km),precip(im)
1443 real denl(im,km),denfacl(im,km),tkl(im,km)
1444 !
1445 integer i,k,n,m,kk,kb,kt,iter
1446 real tl,tl2,qql,dql,qqd
1447 real th,th2,qqh,dqh
1448 real zsum,qsum,dim,dip,c1,con1,fa1,fa2
1449 real allold, allnew, zz, dzamin, cflmax, decfl
1450 real dz(km), ww(km), qq(km), wd(km), wa(km), was(km)
1451 real den(km), denfac(km), tk(km)
1452 real wi(km+1), zi(km+1), za(km+1)
1453 real qn(km), qr(km),tmp(km),tmp1(km),tmp2(km),tmp3(km)
1454 real dza(km+1), qa(km+1), qmi(km+1), qpi(km+1)
1455 !
1456 precip(:) = 0.0
1457 !
1458 i_loop : do i=1,im
1459 ! -----------------------------------
1460 dz(:) = dzl(i,:)
1461 qq(:) = rql(i,:)
1462 ww(:) = wwl(i,:)
1463 den(:) = denl(i,:)
1464 denfac(:) = denfacl(i,:)
1465 tk(:) = tkl(i,:)
1466 ! skip for no precipitation for all layers
1467 allold = 0.0
1468 do k=1,km
1469 allold = allold + qq(k)
1470 enddo
1471 if(allold.le.0.0) then
1472 cycle i_loop
1473 endif
1474 !
1475 ! compute interface values
1476 zi(1)=0.0
1477 do k=1,km
1478 zi(k+1) = zi(k)+dz(k)
1479 enddo
1480 !
1481 ! save departure wind
1482 wd(:) = ww(:)
1483 n=1
1484 100 continue
1485 ! plm is 2nd order, we can use 2nd order wi or 3rd order wi
1486 ! 2nd order interpolation to get wi
1487 wi(1) = ww(1)
1488 wi(km+1) = ww(km)
1489 do k=2,km
1490 wi(k) = (ww(k)*dz(k-1)+ww(k-1)*dz(k))/(dz(k-1)+dz(k))
1491 enddo
1492 ! 3rd order interpolation to get wi
1493 fa1 = 9./16.
1494 fa2 = 1./16.
1495 wi(1) = ww(1)
1496 wi(2) = 0.5*(ww(2)+ww(1))
1497 do k=3,km-1
1498 wi(k) = fa1*(ww(k)+ww(k-1))-fa2*(ww(k+1)+ww(k-2))
1499 enddo
1500 wi(km) = 0.5*(ww(km)+ww(km-1))
1501 wi(km+1) = ww(km)
1502 !
1503 ! terminate of top of raingroup
1504 do k=2,km
1505 if( ww(k).eq.0.0 ) wi(k)=ww(k-1)
1506 enddo
1507 !
1508 ! diffusivity of wi
1509 con1 = 0.05
1510 do k=km,1,-1
1511 decfl = (wi(k+1)-wi(k))*dt/dz(k)
1512 if( decfl .gt. con1 ) then
1513 wi(k) = wi(k+1) - con1*dz(k)/dt
1514 endif
1515 enddo
1516 ! compute arrival point
1517 do k=1,km+1
1518 za(k) = zi(k) - wi(k)*dt
1519 enddo
1520 !
1521 do k=1,km
1522 dza(k) = za(k+1)-za(k)
1523 enddo
1524 dza(km+1) = zi(km+1) - za(km+1)
1525 !
1526 ! computer deformation at arrival point
1527 do k=1,km
1528 qa(k) = qq(k)*dz(k)/dza(k)
1529 qr(k) = qa(k)/den(k)
1530 enddo
1531 qa(km+1) = 0.0
1532 ! call maxmin(km,1,qa,' arrival points ')
1533 !
1534 ! compute arrival terminal velocity, and estimate mean terminal velocity
1535 ! then back to use mean terminal velocity
1536 if( n.le.iter ) then
1537 if (id.eq.1) then
1538 call slope_rain(qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,1,1,1,km)
1539 else
1540 call slope_snow(qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,1,1,1,km)
1541 endif
1542 if( n.ge.2 ) wa(1:km)=0.5*(wa(1:km)+was(1:km))
1543 do k=1,km
1544 !#ifdef DEBUG
1545 ! print*,' slope_wsm3 ',qr(k)*1000.,den(k),denfac(k),tk(k),tmp(k),tmp1(k),tmp2(k),ww(k),wa(k)
1546 !#endif
1547 ! mean wind is average of departure and new arrival winds
1548 ww(k) = 0.5* ( wd(k)+wa(k) )
1549 enddo
1550 was(:) = wa(:)
1551 n=n+1
1552 go to 100
1553 endif
1554 !
1555 ! estimate values at arrival cell interface with monotone
1556 do k=2,km
1557 dip=(qa(k+1)-qa(k))/(dza(k+1)+dza(k))
1558 dim=(qa(k)-qa(k-1))/(dza(k-1)+dza(k))
1559 if( dip*dim.le.0.0 ) then
1560 qmi(k)=qa(k)
1561 qpi(k)=qa(k)
1562 else
1563 qpi(k)=qa(k)+0.5*(dip+dim)*dza(k)
1564 qmi(k)=2.0*qa(k)-qpi(k)
1565 if( qpi(k).lt.0.0 .or. qmi(k).lt.0.0 ) then
1566 qpi(k) = qa(k)
1567 qmi(k) = qa(k)
1568 endif
1569 endif
1570 enddo
1571 qpi(1)=qa(1)
1572 qmi(1)=qa(1)
1573 qmi(km+1)=qa(km+1)
1574 qpi(km+1)=qa(km+1)
1575 !
1576 ! interpolation to regular point
1577 qn = 0.0
1578 kb=1
1579 kt=1
1580 intp : do k=1,km
1581 kb=max(kb-1,1)
1582 kt=max(kt-1,1)
1583 ! find kb and kt
1584 if( zi(k).ge.za(km+1) ) then
1585 exit intp
1586 else
1587 find_kb : do kk=kb,km
1588 if( zi(k).le.za(kk+1) ) then
1589 kb = kk
1590 exit find_kb
1591 else
1592 cycle find_kb
1593 endif
1594 enddo find_kb
1595 find_kt : do kk=kt,km
1596 if( zi(k+1).le.za(kk) ) then
1597 kt = kk
1598 exit find_kt
1599 else
1600 cycle find_kt
1601 endif
1602 enddo find_kt
1603 kt = kt - 1
1604 ! compute q with piecewise constant method
1605 if( kt.eq.kb ) then
1606 tl=(zi(k)-za(kb))/dza(kb)
1607 th=(zi(k+1)-za(kb))/dza(kb)
1608 tl2=tl*tl
1609 th2=th*th
1610 qqd=0.5*(qpi(kb)-qmi(kb))
1611 qqh=qqd*th2+qmi(kb)*th
1612 qql=qqd*tl2+qmi(kb)*tl
1613 qn(k) = (qqh-qql)/(th-tl)
1614 else if( kt.gt.kb ) then
1615 tl=(zi(k)-za(kb))/dza(kb)
1616 tl2=tl*tl
1617 qqd=0.5*(qpi(kb)-qmi(kb))
1618 qql=qqd*tl2+qmi(kb)*tl
1619 dql = qa(kb)-qql
1620 zsum = (1.-tl)*dza(kb)
1621 qsum = dql*dza(kb)
1622 if( kt-kb.gt.1 ) then
1623 do m=kb+1,kt-1
1624 zsum = zsum + dza(m)
1625 qsum = qsum + qa(m) * dza(m)
1626 enddo
1627 endif
1628 th=(zi(k+1)-za(kt))/dza(kt)
1629 th2=th*th
1630 qqd=0.5*(qpi(kt)-qmi(kt))
1631 dqh=qqd*th2+qmi(kt)*th
1632 zsum = zsum + th*dza(kt)
1633 qsum = qsum + dqh*dza(kt)
1634 qn(k) = qsum/zsum
1635 endif
1636 cycle intp
1637 endif
1638 !
1639 enddo intp
1640 !
1641 ! rain out
1642 sum_precip: do k=1,km
1643 if( za(k).lt.0.0 .and. za(k+1).lt.0.0 ) then
1644 precip(i) = precip(i) + qa(k)*dza(k)
1645 cycle sum_precip
1646 else if ( za(k).lt.0.0 .and. za(k+1).ge.0.0 ) then
1647 precip(i) = precip(i) + qa(k)*(0.0-za(k))
1648 exit sum_precip
1649 endif
1650 exit sum_precip
1651 enddo sum_precip
1652 !
1653 ! replace the new values
1654 rql(i,:) = qn(:)
1655 !
1656 ! ----------------------------------
1657 enddo i_loop
1658 !
1659 END SUBROUTINE nislfv_rain_plm
1661 # else
1662 # include "mic-wsm5-3-5-code.h"
1663 # endif
1664 ! end of #ifndef XEON_OPTIMIZED_WSM5
1666 ! remainder of routines are common to original and MIC version
1668 !+---+-----------------------------------------------------------------+
1669 subroutine refl10cm_wsm5 (qv1d, qr1d, qs1d, &
1670 t1d, p1d, dBZ, kts, kte, ii, jj)
1672 IMPLICIT NONE
1674 !..Sub arguments
1675 INTEGER, INTENT(IN):: kts, kte, ii, jj
1676 REAL, DIMENSION(kts:kte), INTENT(IN):: &
1677 qv1d, qr1d, qs1d, t1d, p1d
1678 REAL, DIMENSION(kts:kte), INTENT(INOUT):: dBZ
1680 !..Local variables
1681 REAL, DIMENSION(kts:kte):: temp, pres, qv, rho
1682 REAL, DIMENSION(kts:kte):: rr, rs
1683 REAL:: temp_C
1685 DOUBLE PRECISION, DIMENSION(kts:kte):: ilamr, ilams
1686 DOUBLE PRECISION, DIMENSION(kts:kte):: N0_r, N0_s
1687 DOUBLE PRECISION:: lamr, lams
1688 LOGICAL, DIMENSION(kts:kte):: L_qr, L_qs
1690 REAL, DIMENSION(kts:kte):: ze_rain, ze_snow
1691 DOUBLE PRECISION:: fmelt_s
1693 INTEGER:: i, k, k_0, kbot, n
1694 LOGICAL:: melti
1696 DOUBLE PRECISION:: cback, x, eta, f_d
1697 REAL, PARAMETER:: R=287.
1699 !+---+
1701 do k = kts, kte
1702 dBZ(k) = -35.0
1703 enddo
1705 !+---+-----------------------------------------------------------------+
1706 !..Put column of data into local arrays.
1707 !+---+-----------------------------------------------------------------+
1708 do k = kts, kte
1709 temp(k) = t1d(k)
1710 temp_C = min(-0.001, temp(K)-273.15)
1711 qv(k) = MAX(1.E-10, qv1d(k))
1712 pres(k) = p1d(k)
1713 rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
1715 if (qr1d(k) .gt. 1.E-9) then
1716 rr(k) = qr1d(k)*rho(k)
1717 N0_r(k) = n0r
1718 lamr = (xam_r*xcrg(3)*N0_r(k)/rr(k))**(1./xcre(1))
1719 ilamr(k) = 1./lamr
1720 L_qr(k) = .true.
1721 else
1722 rr(k) = 1.E-12
1723 L_qr(k) = .false.
1724 endif
1726 if (qs1d(k) .gt. 1.E-9) then
1727 rs(k) = qs1d(k)*rho(k)
1728 N0_s(k) = min(n0smax, n0s*exp(-alpha*temp_C))
1729 lams = (xam_s*xcsg(3)*N0_s(k)/rs(k))**(1./xcse(1))
1730 ilams(k) = 1./lams
1731 L_qs(k) = .true.
1732 else
1733 rs(k) = 1.E-12
1734 L_qs(k) = .false.
1735 endif
1736 enddo
1738 !+---+-----------------------------------------------------------------+
1739 !..Locate K-level of start of melting (k_0 is level above).
1740 !+---+-----------------------------------------------------------------+
1741 melti = .false.
1742 k_0 = kts
1743 do k = kte-1, kts, -1
1744 if ( (temp(k).gt.273.15) .and. L_qr(k) .and. L_qs(k+1) ) then
1745 k_0 = MAX(k+1, k_0)
1746 melti=.true.
1747 goto 195
1748 endif
1749 enddo
1750 195 continue
1752 !+---+-----------------------------------------------------------------+
1753 !..Assume Rayleigh approximation at 10 cm wavelength. Rain (all temps)
1754 !.. and non-water-coated snow and graupel when below freezing are
1755 !.. simple. Integrations of m(D)*m(D)*N(D)*dD.
1756 !+---+-----------------------------------------------------------------+
1758 do k = kts, kte
1759 ze_rain(k) = 1.e-22
1760 ze_snow(k) = 1.e-22
1761 if (L_qr(k)) ze_rain(k) = N0_r(k)*xcrg(4)*ilamr(k)**xcre(4)
1762 if (L_qs(k)) ze_snow(k) = (0.176/0.93) * (6.0/PI)*(6.0/PI) &
1763 * (xam_s/900.0)*(xam_s/900.0) &
1764 * N0_s(k)*xcsg(4)*ilams(k)**xcse(4)
1765 enddo
1768 !+---+-----------------------------------------------------------------+
1769 !..Special case of melting ice (snow/graupel) particles. Assume the
1770 !.. ice is surrounded by the liquid water. Fraction of meltwater is
1771 !.. extremely simple based on amount found above the melting level.
1772 !.. Uses code from Uli Blahak (rayleigh_soak_wetgraupel and supporting
1773 !.. routines).
1774 !+---+-----------------------------------------------------------------+
1776 if (melti .and. k_0.ge.kts+1) then
1777 do k = k_0-1, kts, -1
1779 !..Reflectivity contributed by melting snow
1780 if (L_qs(k) .and. L_qs(k_0) ) then
1781 fmelt_s = MAX(0.005d0, MIN(1.0d0-rs(k)/rs(k_0), 0.99d0))
1782 eta = 0.d0
1783 lams = 1./ilams(k)
1784 do n = 1, nrbins
1785 x = xam_s * xxDs(n)**xbm_s
1786 call rayleigh_soak_wetgraupel (x,DBLE(xocms),DBLE(xobms), &
1787 fmelt_s, melt_outside_s, m_w_0, m_i_0, lamda_radar, &
1788 CBACK, mixingrulestring_s, matrixstring_s, &
1789 inclusionstring_s, hoststring_s, &
1790 hostmatrixstring_s, hostinclusionstring_s)
1791 f_d = N0_s(k)*xxDs(n)**xmu_s * DEXP(-lams*xxDs(n))
1792 eta = eta + f_d * CBACK * simpson(n) * xdts(n)
1793 enddo
1794 ze_snow(k) = SNGL(lamda4 / (pi5 * K_w) * eta)
1795 endif
1796 enddo
1797 endif
1799 do k = kte, kts, -1
1800 dBZ(k) = 10.*log10((ze_rain(k)+ze_snow(k))*1.d18)
1801 enddo
1803 end subroutine refl10cm_wsm5
1804 !+---+-----------------------------------------------------------------+
1805 !-------------------------------------------------------------------
1806 SUBROUTINE wsm5init(den0,denr,dens,cl,cpv,allowed_to_read)
1807 !-------------------------------------------------------------------
1808 IMPLICIT NONE
1809 !-------------------------------------------------------------------
1810 !.... constants which may not be tunable
1811 REAL, INTENT(IN) :: den0,denr,dens,cl,cpv
1812 LOGICAL, INTENT(IN) :: allowed_to_read
1813 !
1814 pi = 4.*atan(1.)
1815 xlv1 = cl-cpv
1816 !
1817 qc0 = 4./3.*pi*denr*r0**3*xncr/den0 ! 0.419e-3 -- .61e-3
1818 qck1 = .104*9.8*peaut/(xncr*denr)**(1./3.)/xmyu*den0**(4./3.) ! 7.03
1819 pidnc = pi*denr/6. ! syb
1820 !
1821 bvtr1 = 1.+bvtr
1822 bvtr2 = 2.5+.5*bvtr
1823 bvtr3 = 3.+bvtr
1824 bvtr4 = 4.+bvtr
1825 g1pbr = rgmma(bvtr1)
1826 g3pbr = rgmma(bvtr3)
1827 g4pbr = rgmma(bvtr4) ! 17.837825
1828 g5pbro2 = rgmma(bvtr2) ! 1.8273
1829 pvtr = avtr*g4pbr/6.
1830 eacrr = 1.0
1831 pacrr = pi*n0r*avtr*g3pbr*.25*eacrr
1832 precr1 = 2.*pi*n0r*.78
1833 precr2 = 2.*pi*n0r*.31*avtr**.5*g5pbro2
1834 xmmax = (dimax/dicon)**2
1835 roqimax = 2.08e22*dimax**8
1836 !
1837 bvts1 = 1.+bvts
1838 bvts2 = 2.5+.5*bvts
1839 bvts3 = 3.+bvts
1840 bvts4 = 4.+bvts
1841 g1pbs = rgmma(bvts1) !.8875
1842 g3pbs = rgmma(bvts3)
1843 g4pbs = rgmma(bvts4) ! 12.0786
1844 g5pbso2 = rgmma(bvts2)
1845 pvts = avts*g4pbs/6.
1846 pacrs = pi*n0s*avts*g3pbs*.25
1847 precs1 = 4.*n0s*.65
1848 precs2 = 4.*n0s*.44*avts**.5*g5pbso2
1849 pidn0r = pi*denr*n0r
1850 pidn0s = pi*dens*n0s
1851 pacrc = pi*n0s*avts*g3pbs*.25*eacrc
1852 !
1853 rslopermax = 1./lamdarmax
1854 rslopesmax = 1./lamdasmax
1855 rsloperbmax = rslopermax ** bvtr
1856 rslopesbmax = rslopesmax ** bvts
1857 rsloper2max = rslopermax * rslopermax
1858 rslopes2max = rslopesmax * rslopesmax
1859 rsloper3max = rsloper2max * rslopermax
1860 rslopes3max = rslopes2max * rslopesmax
1861 !
1862 !+---+-----------------------------------------------------------------+
1863 !..Set these variables needed for computing radar reflectivity. These
1864 !.. get used within radar_init to create other variables used in the
1865 !.. radar module.
1866 xam_r = PI*denr/6.
1867 xbm_r = 3.
1868 xmu_r = 0.
1869 xam_s = PI*dens/6.
1870 xbm_s = 3.
1871 xmu_s = 0.
1872 xam_g = PI*dens/6. ! These 3 variables for graupel are set but unused.
1873 xbm_g = 3.
1874 xmu_g = 0.
1876 call radar_init
1877 !+---+-----------------------------------------------------------------+
1879 END SUBROUTINE wsm5init
1880 ! ...................................................................
1881 REAL FUNCTION rgmma(x)
1882 !-------------------------------------------------------------------
1883 IMPLICIT NONE
1884 !-------------------------------------------------------------------
1885 ! rgmma function: use infinite product form
1886 REAL :: euler
1887 PARAMETER (euler=0.577215664901532)
1888 REAL :: x, y
1889 INTEGER :: i
1890 if(x.eq.1.)then
1891 rgmma=0.
1892 else
1893 rgmma=x*exp(euler*x)
1894 do i=1,10000
1895 y=float(i)
1896 rgmma=rgmma*(1.000+x/y)*exp(-x/y)
1897 enddo
1898 rgmma=1./rgmma
1899 endif
1900 END FUNCTION rgmma
1902 !-----------------------------------------------------------------------
1903 subroutine effectRad_wsm5 (t, qc, qi, qs, rho, qmin, t0c, &
1904 re_qc, re_qi, re_qs, kts, kte, ii, jj)
1906 !-----------------------------------------------------------------------
1907 ! Compute radiation effective radii of cloud water, ice, and snow for
1908 ! single-moment microphysics.
1909 ! These are entirely consistent with microphysics assumptions, not
1910 ! constant or otherwise ad hoc as is internal to most radiation
1911 ! schemes.
1912 ! Coded and implemented by Soo ya Bae, KIAPS, January 2015.
1913 !-----------------------------------------------------------------------
1915 implicit none
1917 !..Sub arguments
1918 integer, intent(in) :: kts, kte, ii, jj
1919 real, intent(in) :: qmin
1920 real, intent(in) :: t0c
1921 real, dimension( kts:kte ), intent(in):: t
1922 real, dimension( kts:kte ), intent(in):: qc
1923 real, dimension( kts:kte ), intent(in):: qi
1924 real, dimension( kts:kte ), intent(in):: qs
1925 real, dimension( kts:kte ), intent(in):: rho
1926 real, dimension( kts:kte ), intent(inout):: re_qc
1927 real, dimension( kts:kte ), intent(inout):: re_qi
1928 real, dimension( kts:kte ), intent(inout):: re_qs
1929 !..Local variables
1930 integer:: i,k
1931 integer :: inu_c
1932 real, dimension( kts:kte ):: ni
1933 real, dimension( kts:kte ):: rqc
1934 real, dimension( kts:kte ):: rqi
1935 real, dimension( kts:kte ):: rni
1936 real, dimension( kts:kte ):: rqs
1937 real :: temp
1938 real :: lamdac
1939 real :: supcol, n0sfac, lamdas
1940 real :: diai ! diameter of ice in m
1941 logical :: has_qc, has_qi, has_qs
1942 !..Minimum microphys values
1943 real, parameter :: R1 = 1.E-12
1944 real, parameter :: R2 = 1.E-6
1945 !..Mass power law relations: mass = am*D**bm
1946 real, parameter :: bm_r = 3.0
1947 real, parameter :: obmr = 1.0/bm_r
1948 real, parameter :: nc0 = 3.E8
1949 !-----------------------------------------------------------------------
1950 has_qc = .false.
1951 has_qi = .false.
1952 has_qs = .false.
1954 do k = kts, kte
1955 ! for cloud
1956 rqc(k) = max(R1, qc(k)*rho(k))
1957 if (rqc(k).gt.R1) has_qc = .true.
1958 ! for ice
1959 rqi(k) = max(R1, qi(k)*rho(k))
1960 temp = (rho(k)*max(qi(k),qmin))
1961 temp = sqrt(sqrt(temp*temp*temp))
1962 ni(k) = min(max(5.38e7*temp,1.e3),1.e6)
1963 rni(k)= max(R2, ni(k)*rho(k))
1964 if (rqi(k).gt.R1 .and. rni(k).gt.R2) has_qi = .true.
1965 ! for snow
1966 rqs(k) = max(R1, qs(k)*rho(k))
1967 if (rqs(k).gt.R1) has_qs = .true.
1968 enddo
1970 if (has_qc) then
1971 do k=kts,kte
1972 if (rqc(k).le.R1) CYCLE
1973 lamdac = (pidnc*nc0/rqc(k))**obmr
1974 re_qc(k) = max(2.51E-6,min(1.5*(1.0/lamdac),50.E-6))
1975 enddo
1976 endif
1978 if (has_qi) then
1979 do k=kts,kte
1980 if (rqi(k).le.R1 .or. rni(k).le.R2) CYCLE
1981 diai = 11.9*sqrt(rqi(k)/ni(k))
1982 re_qi(k) = max(10.01E-6,min(0.75*0.163*diai,125.E-6))
1983 enddo
1984 endif
1986 if (has_qs) then
1987 do k=kts,kte
1988 if (rqs(k).le.R1) CYCLE
1989 supcol = t0c-t(k)
1990 n0sfac = max(min(exp(alpha*supcol),n0smax/n0s),1.)
1991 lamdas = sqrt(sqrt(pidn0s*n0sfac/rqs(k)))
1992 re_qs(k) = max(25.E-6,min(0.5*(1./lamdas), 999.E-6))
1993 enddo
1994 endif
1996 end subroutine effectRad_wsm5
1997 !-----------------------------------------------------------------------
1999 END MODULE module_mp_wsm5
2000 #endif