| blob | 085ef2f441c780136292112ccfc32f9c1e4027bb |
1 #if ( RWORDSIZE == 4 )
2 # define VREC vsrec
3 # define VSQRT vssqrt
4 #else
5 # define VREC vrec
6 # define VSQRT vsqrt
7 #endif
9 !Including inline expansion statistical function
10 MODULE module_mp_wsm5
11 !
12 !
13 REAL, PARAMETER, PRIVATE :: dtcldcr = 120. ! maximum time step for minor loops
14 REAL, PARAMETER, PRIVATE :: n0r = 8.e6 ! intercept parameter rain
15 REAL, PARAMETER, PRIVATE :: avtr = 841.9 ! a constant for terminal velocity of rain
16 REAL, PARAMETER, PRIVATE :: bvtr = 0.8 ! a constant for terminal velocity of rain
17 REAL, PARAMETER, PRIVATE :: r0 = .8e-5 ! 8 microm in contrast to 10 micro m
18 REAL, PARAMETER, PRIVATE :: peaut = .55 ! collection efficiency
19 REAL, PARAMETER, PRIVATE :: xncr = 3.e8 ! maritime cloud in contrast to 3.e8 in tc80
20 REAL, PARAMETER, PRIVATE :: xmyu = 1.718e-5 ! the dynamic viscosity kgm-1s-1
21 REAL, PARAMETER, PRIVATE :: avts = 11.72 ! a constant for terminal velocity of snow
22 REAL, PARAMETER, PRIVATE :: bvts = .41 ! a constant for terminal velocity of snow
23 REAL, PARAMETER, PRIVATE :: n0smax = 1.e11 ! maximum n0s (t=-90C unlimited)
24 REAL, PARAMETER, PRIVATE :: lamdarmax = 8.e4 ! limited maximum value for slope parameter of rain
25 REAL, PARAMETER, PRIVATE :: lamdasmax = 1.e5 ! limited maximum value for slope parameter of snow
26 REAL, PARAMETER, PRIVATE :: lamdagmax = 6.e4 ! limited maximum value for slope parameter of graupel
27 REAL, PARAMETER, PRIVATE :: dicon = 11.9 ! constant for the cloud-ice diamter
28 REAL, PARAMETER, PRIVATE :: dimax = 500.e-6 ! limited maximum value for the cloud-ice diamter
29 REAL, PARAMETER, PRIVATE :: n0s = 2.e6 ! temperature dependent intercept parameter snow
30 REAL, PARAMETER, PRIVATE :: alpha = .12 ! .122 exponen factor for n0s
31 REAL, PARAMETER, PRIVATE :: pfrz1 = 100. ! constant in Biggs freezing
32 REAL, PARAMETER, PRIVATE :: pfrz2 = 0.66 ! constant in Biggs freezing
33 REAL, PARAMETER, PRIVATE :: qcrmin = 1.e-9 ! minimun values for qr, qs, and qg
34 REAL, PARAMETER, PRIVATE :: eacrc = 1.0 ! Snow/cloud-water collection efficiency
35 REAL, SAVE :: &
36 qc0, qck1,bvtr1,bvtr2,bvtr3,bvtr4,g1pbr, &
37 g3pbr,g4pbr,g5pbro2,pvtr,eacrr,pacrr, &
38 precr1,precr2,xmmax,roqimax,bvts1, &
39 bvts2,bvts3,bvts4,g1pbs,g3pbs,g4pbs, &
40 g5pbso2,pvts,pacrs,precs1,precs2,pidn0r, &
41 pidn0s,xlv1,pacrc, &
42 rslopermax,rslopesmax,rslopegmax, &
43 rsloperbmax,rslopesbmax,rslopegbmax, &
44 rsloper2max,rslopes2max,rslopeg2max, &
45 rsloper3max,rslopes3max,rslopeg3max
46 !
47 ! Specifies code-inlining of fpvs function in WSM52D below. JM 20040507
48 !
49 CONTAINS
50 !===================================================================
51 !
52 SUBROUTINE wsm5(th, q, qc, qr, qi, qs &
53 ,den, pii, p, delz &
54 ,delt,g, cpd, cpv, rd, rv, t0c &
55 ,ep1, ep2, qmin &
56 ,XLS, XLV0, XLF0, den0, denr &
57 ,cliq,cice,psat &
58 ,rain, rainncv &
59 ,snow, snowncv &
60 ,sr &
61 ,ids,ide, jds,jde, kds,kde &
62 ,ims,ime, jms,jme, kms,kme &
63 ,its,ite, jts,jte, kts,kte &
64 )
65 #ifdef _OPENMP
66 use omp_lib
67 #endif
68 !-------------------------------------------------------------------
69 IMPLICIT NONE
70 !-------------------------------------------------------------------
71 !
72 ! This code is a 5-class mixed ice microphyiscs scheme (WSM5) of the WRF
73 ! Single-Moment MicroPhyiscs (WSMMP). The WSMMP assumes that ice nuclei
74 ! number concentration is a function of temperature, and seperate assumption
75 ! is developed, in which ice crystal number concentration is a function
76 ! of ice amount. A theoretical background of the ice-microphysics and related
77 ! processes in the WSMMPs are described in Hong et al. (2004).
78 ! Production terms in the WSM6 scheme are described in Hong and Lim (2006).
79 ! All units are in m.k.s. and source/sink terms in kgkg-1s-1.
80 !
81 ! WSM5 cloud scheme
82 !
83 ! Coded by Song-You Hong (Yonsei Univ.)
84 ! Jimy Dudhia (NCAR) and Shu-Hua Chen (UC Davis)
85 ! Summer 2002
86 !
87 ! Implemented by Song-You Hong (Yonsei Univ.) and Jimy Dudhia (NCAR)
88 ! Summer 2003
89 !
90 ! Reference) Hong, Dudhia, Chen (HDC, 2004) Mon. Wea. Rev.
91 ! Rutledge, Hobbs (RH83, 1983) J. Atmos. Sci.
92 ! Hong and Lim (HL, 2006) J. Korean Meteor. Soc.
93 !
94 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
95 ims,ime, jms,jme, kms,kme , &
96 its,ite, jts,jte, kts,kte
97 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
98 INTENT(INOUT) :: &
99 th, &
100 q, &
101 qc, &
102 qi, &
103 qr, &
104 qs
105 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
106 INTENT(IN ) :: &
107 den, &
108 pii, &
109 p, &
110 delz
111 REAL, INTENT(IN ) :: delt, &
112 g, &
113 rd, &
114 rv, &
115 t0c, &
116 den0, &
117 cpd, &
118 cpv, &
119 ep1, &
120 ep2, &
121 qmin, &
122 XLS, &
123 XLV0, &
124 XLF0, &
125 cliq, &
126 cice, &
127 psat, &
128 denr
129 REAL, DIMENSION( ims:ime , jms:jme ), &
130 INTENT(INOUT) :: rain, &
131 rainncv, &
132 sr
134 REAL, DIMENSION( ims:ime , jms:jme ), OPTIONAL, &
135 INTENT(INOUT) :: snow, &
136 snowncv
138 ! LOCAL VAR
139 REAL, DIMENSION( its:ite , kts:kte ) :: t
140 REAL, DIMENSION( its:ite , kts:kte, 2 ) :: qci, qrs
141 CHARACTER*256 :: emess
142 INTEGER :: mkx_test
143 INTEGER :: i,j,k
145 #ifdef _ACCEL_PROF
146 INTEGER :: l
147 real*8 wsm3_t(8,256), wsm5_t(8,256), t1, t2
148 common /wsm_times/ wsm3_t(8,256), wsm5_t(8,256)
149 #endif
152 !-------------------------------------------------------------------
154 #ifdef _ACCEL_PROF
155 call cpu_time(t1)
156 #endif
158 #ifndef RUN_ON_GPU
160 #ifdef _ACCEL
162 ! Need to send th, pii, qc, qi, qr, qs
163 ! Don't send t
165 CALL wsm52D(th, pii, q, qc, qr, qi, qs &
166 ,den &
167 ,p, delz &
168 ,delt,g, cpd, cpv, rd, rv, t0c &
169 ,ep1, ep2, qmin &
170 ,XLS, XLV0, XLF0, den0, denr &
171 ,cliq,cice,psat &
172 ,rain,rainncv &
173 ,sr &
174 ,ids,ide, jds,jde, kds,kde &
175 ,ims,ime, jms,jme, kms,kme &
176 ,its,ite, jts,jte, kts,kte &
177 ,snow,snowncv &
178 )
180 #else
182 DO j=jts,jte
183 DO k=kts,kte
184 DO i=its,ite
185 t(i,k)=th(i,k,j)*pii(i,k,j)
186 qci(i,k,1) = qc(i,k,j)
187 qci(i,k,2) = qi(i,k,j)
188 qrs(i,k,1) = qr(i,k,j)
189 qrs(i,k,2) = qs(i,k,j)
190 ENDDO
191 ENDDO
193 ! Sending array starting locations of optional variables may cause
194 ! troubles, so we explicitly change the call.
196 CALL wsm52D(t, q(ims,kms,j), qci, qrs &
197 ,den(ims,kms,j) &
198 ,p(ims,kms,j), delz(ims,kms,j) &
199 ,delt,g, cpd, cpv, rd, rv, t0c &
200 ,ep1, ep2, qmin &
201 ,XLS, XLV0, XLF0, den0, denr &
202 ,cliq,cice,psat &
203 ,j &
204 ,rain(ims,j),rainncv(ims,j) &
205 ,sr(ims,j) &
206 ,ids,ide, jds,jde, kds,kde &
207 ,ims,ime, jms,jme, kms,kme &
208 ,its,ite, jts,jte, kts,kte &
209 ,snow,snowncv &
210 )
212 DO K=kts,kte
213 DO I=its,ite
214 th(i,k,j)=t(i,k)/pii(i,k,j)
215 qc(i,k,j) = qci(i,k,1)
216 qi(i,k,j) = qci(i,k,2)
217 qr(i,k,j) = qrs(i,k,1)
218 qs(i,k,j) = qrs(i,k,2)
219 ENDDO
220 ENDDO
221 ENDDO
222 #endif
223 #else
224 CALL get_wsm5_gpu_levels ( mkx_test )
225 IF ( mkx_test .LT. kte ) THEN
226 WRITE(emess,*)'Number of levels compiled for GPU WSM5 too small. ', &
227 mkx_test,' < ',kte
228 CALL wrf_error_fatal(emess)
229 ENDIF
230 CALL wsm5_host ( &
231 th(its:ite,kts:kte,jts:jte), pii(its:ite,kts:kte,jts:jte) &
232 ,q(its:ite,kts:kte,jts:jte), qc(its:ite,kts:kte,jts:jte) &
233 ,qi(its:ite,kts:kte,jts:jte), qr(its:ite,kts:kte,jts:jte) &
234 ,qs(its:ite,kts:kte,jts:jte), den(its:ite,kts:kte,jts:jte) &
235 ,p(its:ite,kts:kte,jts:jte), delz(its:ite,kts:kte,jts:jte) &
236 ,delt &
237 ,rain(its:ite,jts:jte),rainncv(its:ite,jts:jte) &
238 ,snow(its:ite,jts:jte),snowncv(its:ite,jts:jte) &
239 ,sr(its:ite,jts:jte) &
240 ,its, ite, jts, jte, kts, kte &
241 ,its, ite, jts, jte, kts, kte &
242 ,its, ite, jts, jte, kts, kte &
243 )
244 #endif
247 #ifdef _ACCEL_PROF
248 call cpu_time(t2)
249 #ifdef _OPENMP
250 l = omp_get_thread_num() + 1
251 #else
252 l = 1
253 #endif
254 wsm5_t(1,l) = wsm5_t(1,l) + (t2 - t1)
255 #endif
257 END SUBROUTINE wsm5
260 #ifdef _ACCEL
262 !===================================================================
263 !
264 SUBROUTINE wsm52D(th, pii, q, qc, qr, qi, qqs, den, p, delz &
265 ,delt,g, cpd, cpv, rd, rv, t0c &
266 ,ep1, ep2, qmin &
267 ,XLS, XLV0, XLF0, den0, denr &
268 ,cliq,cice,psat &
269 ,rain,rainncv &
270 ,sr &
271 ,ids,ide, jds,jde, kds,kde &
272 ,ims,ime, jms,jme, kms,kme &
273 ,its,ite, jts,jte, kts,kte &
274 ,snow,snowncv &
275 )
276 !-------------------------------------------------------------------
277 IMPLICIT NONE
278 !-------------------------------------------------------------------
279 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
280 ims,ime, jms,jme, kms,kme , &
281 its,ite, jts,jte, kts,kte
282 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
283 INTENT(INOUT) :: &
284 th
285 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
286 INTENT(IN) :: &
287 pii
288 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
289 INTENT(INOUT) :: &
290 qc, &
291 qr, &
292 qi, &
293 qqs
294 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
295 INTENT(INOUT) :: &
296 q
297 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
298 INTENT(IN ) :: &
299 den, &
300 p, &
301 delz
302 REAL, INTENT(IN ) :: delt, &
303 g, &
304 cpd, &
305 cpv, &
306 t0c, &
307 den0, &
308 rd, &
309 rv, &
310 ep1, &
311 ep2, &
312 qmin, &
313 XLS, &
314 XLV0, &
315 XLF0, &
316 cliq, &
317 cice, &
318 psat, &
319 denr
320 REAL, DIMENSION( ims:ime, jms:jme ), &
321 INTENT(INOUT) :: rain, &
322 rainncv, &
323 sr
325 REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, &
326 INTENT(INOUT) :: snow, &
327 snowncv
329 ! LOCAL VAR
330 REAL, DIMENSION( its:ite , kts:kte , 2) :: &
331 rh, &
332 qs, &
333 rslope, &
334 rslope2, &
335 rslope3, &
336 rslopeb, &
337 falk, &
338 fall, &
339 work1
340 REAL, DIMENSION( its:ite , kts:kte, jts:jte ) :: &
341 t
342 REAL, DIMENSION( its:ite , kts:kte , 2 ) :: &
343 qci, &
344 qrs
345 REAL, DIMENSION( its:ite , kts:kte ) :: &
346 falkc, &
347 fallc, &
348 xl, &
349 cpm, &
350 denfac, &
351 xni, &
352 n0sfac, &
353 work2, &
354 work1c, &
355 work2c
357 REAL, DIMENSION( its:ite , kts:kte ) :: &
358 pigen, &
359 pidep, &
360 psdep, &
361 praut, &
362 psaut, &
363 prevp, &
364 psevp, &
365 pracw, &
366 psacw, &
367 psaci, &
368 pcond, &
369 psmlt
370 INTEGER :: &
371 mstep, &
372 numdt
373 REAL :: rmstep
374 REAL dtcldden, rdelz, rdtcld
375 LOGICAL :: flgcld
377 #define WSM_NO_CONDITIONAL_IN_VECTOR
378 #ifdef WSM_NO_CONDITIONAL_IN_VECTOR
379 REAL :: xal, xbl
380 #endif
382 REAL :: pi, &
383 cpmcal, xlcal, lamdar, lamdas, diffus, &
384 viscos, xka, venfac, conden, diffac, &
385 x, y, z, a, b, c, d, e, &
386 qdt, holdrr, holdrs, supcol, supcolt, pvt, &
387 coeres, supsat, dtcld, xmi, eacrs, satdt, &
388 vt2i,vt2s,acrfac, &
389 qimax, diameter, xni0, roqi0, &
390 fallsum, fallsum_qsi, xlwork2, factor, source, &
391 value, xlf, pfrzdtc, pfrzdtr, supice, holdc, holdci
392 ! variables for optimization
393 REAL, DIMENSION( its:ite ) :: tvec1
394 REAL :: temp
395 INTEGER :: i, j, k, &
396 iprt, latd, lond, loop, loops, ifsat, n
397 ! Temporaries used for inlining fpvs function
398 REAL :: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
399 REAL :: logtr
400 !
401 !=================================================================
402 ! compute internal functions
403 !
404 cpmcal(x) = cpd*(1.-max(x,qmin))+max(x,qmin)*cpv
405 xlcal(x) = xlv0-xlv1*(x-t0c)
406 !----------------------------------------------------------------
407 ! size distributions: (x=mixing ratio, y=air density):
408 ! valid for mixing ratio > 1.e-9 kg/kg.
409 !
410 ! Optimizatin : A**B => exp(log(A)*(B))
411 lamdar(x,y)= sqrt(sqrt(pidn0r/(x*y))) ! (pidn0r/(x*y))**.25
412 lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y))) ! (pidn0s*z/(x*y))**.25
413 !
414 !----------------------------------------------------------------
415 ! diffus: diffusion coefficient of the water vapor
416 ! viscos: kinematic viscosity(m2s-1)
417 ! diffus(x,y) = 8.794e-5 * exp(log(x)*(1.81)) / y ! 8.794e-5*x**1.81/y
418 ! viscos(x,y) = 1.496e-6 * (x*sqrt(x)) /(x+120.)/y ! 1.496e-6*x**1.5/(x+120.)/y
419 ! xka(x,y) = 1.414e3*viscos(x,y)*y
420 ! diffac(a,b,c,d,e) = d*a*a/(xka(c,d)*rv*c*c)+1./(e*diffus(c,b))
421 ! venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333))) &
422 ! /sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
423 ! conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+d*d/(rv*e)*c/(a*a))
424 !
425 !
426 pi = 4. * atan(1.)
427 !
428 !----------------------------------------------------------------
429 ! paddint 0 for negative values generated by dynamics
430 !
432 !
433 ! Moved outside of accelerator region
434 !
435 loops = max(nint(delt/dtcldcr),1)
436 dtcld = delt/loops
437 if(delt.le.dtcldcr) dtcld = delt
438 !
440 !....!$acc local(t) &
442 IF ( PRESENT (snowncv) .AND. PRESENT (snow)) THEN
444 !$acc region &
445 !$acc local(t) &
446 !$acc copyin(delz(:,:,:),p(:,:,:),den(:,:,:),pii(:,:,:)) &
447 !$acc copyout(snowncv(:,:),rainncv(:,:),sr(:,:)) &
448 !$acc copy(qqs(:,:,:),qr(:,:,:),qi(:,:,:),qc(:,:,:)) &
449 !$acc copy(th(:,:,:),q(:,:,:),snow(:,:),rain(:,:))
450 !$acc do &
451 !$acc private(rh,qs,rslope,rslope2,rslope3,rslopeb,falk,fall) &
452 !$acc private(work1,qci,qrs,falkc,fallc,xl,cpm,denfac,xni) &
453 !$acc private(n0sfac,work2,work1c,work2c,pigen,pidep,psdep) &
454 !$acc private(praut,psaut,prevp,psevp) &
455 !$acc private(pracw,psacw,psaci,pcond,psmlt) &
456 !$acc parallel
457 do j = jts, jte
458 !$acc do &
459 !$acc private(numdt,mstep) &
460 !$acc kernel vector
461 do i = its, ite
462 do k = kts, kte
463 t(i,k,j)=th(i,k,j)*pii(i,k,j)
464 qci(i,k,1) = max(qc(i,k,j),0.0)
465 qci(i,k,2) = max(qi(i,k,j),0.0)
466 qrs(i,k,1) = max(qr(i,k,j),0.0)
467 qrs(i,k,2) = max(qqs(i,k,j),0.0)
468 enddo
469 !
470 !----------------------------------------------------------------
471 ! latent heat for phase changes and heat capacity. neglect the
472 ! changes during microphysical process calculation
473 ! emanuel(1994)
474 !
475 do k = kts, kte
476 cpm(i,k) = cpmcal(q(i,k,j))
477 xl(i,k) = xlcal(t(i,k,j))
478 enddo
479 !
480 !----------------------------------------------------------------
481 ! compute the minor time steps.
482 !
483 ! loops = max(nint(delt/dtcldcr),1)
484 ! dtcld = delt/loops
485 ! if(delt.le.dtcldcr) dtcld = delt
486 !
487 do loop = 1,loops
488 !
489 !----------------------------------------------------------------
490 ! initialize the large scale variables
491 !
492 mstep = 1
493 flgcld = .true.
494 !
495 do k = kts, kte
496 denfac(i,k) = sqrt(den0/den(i,k,j))
497 enddo
498 ! do k = kts, kte
499 ! CALL VREC( tvec1(its), den(its,k,j), ite-its+1)
500 ! do i = its, ite
501 ! tvec1(i) = tvec1(i)*den0
502 ! enddo
503 ! CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
504 ! enddo
505 !
506 ! Inline expansion for fpvs
507 ! qs(i,k,1) = fpvs(t(i,k,j),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
508 ! qs(i,k,2) = fpvs(t(i,k,j),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
509 hsub = xls
510 hvap = xlv0
511 cvap = cpv
512 ttp=t0c+0.01
513 dldt=cvap-cliq
514 xa=-dldt/rv
515 xb=xa+hvap/(rv*ttp)
516 dldti=cvap-cice
517 xai=-dldti/rv
518 xbi=xai+hsub/(rv*ttp)
520 ! this is for compilers where the conditional inhibits vectorization
521 #ifdef WSM_NO_CONDITIONAL_IN_VECTOR
522 do k = kts, kte
523 if(t(i,k,j).lt.ttp) then
524 xal = xai
525 xbl = xbi
526 else
527 xal = xa
528 xbl = xb
529 endif
530 tr=ttp/t(i,k,j)
531 logtr=log(tr)
532 qs(i,k,1)=psat*exp(logtr*(xa)+xb*(1.-tr))
533 qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k,j) - qs(i,k,1))
534 qs(i,k,1) = max(qs(i,k,1),qmin)
535 rh(i,k,1) = max(q(i,k,j) / qs(i,k,1),qmin)
536 qs(i,k,2)=psat*exp(logtr*(xal)+xbl*(1.-tr))
537 qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k,j) - qs(i,k,2))
538 qs(i,k,2) = max(qs(i,k,2),qmin)
539 rh(i,k,2) = max(q(i,k,j) / qs(i,k,2),qmin)
540 enddo
541 #else
542 do k = kts, kte
543 tr=ttp/t(i,k,j)
544 logtr=log(tr)
545 qs(i,k,1)=psat*exp(logtr*(xa)+xb*(1.-tr))
546 qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k,j) - qs(i,k,1))
547 qs(i,k,1) = max(qs(i,k,1),qmin)
548 rh(i,k,1) = max(q(i,k,j) / qs(i,k,1),qmin)
549 if(t(i,k,j).lt.ttp) then
550 qs(i,k,2)=psat*exp(logtr*(xai)+xbi*(1.-tr))
551 else
552 qs(i,k,2)=psat*exp(logtr*(xa)+xb*(1.-tr))
553 endif
554 qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k,j) - qs(i,k,2))
555 qs(i,k,2) = max(qs(i,k,2),qmin)
556 rh(i,k,2) = max(q(i,k,j) / qs(i,k,2),qmin)
557 enddo
558 #endif
560 !
561 !----------------------------------------------------------------
562 ! initialize the variables for microphysical physics
563 !
564 !
565 do k = kts, kte
566 prevp(i,k) = 0.
567 psdep(i,k) = 0.
568 praut(i,k) = 0.
569 psaut(i,k) = 0.
570 pracw(i,k) = 0.
571 psaci(i,k) = 0.
572 psacw(i,k) = 0.
573 pigen(i,k) = 0.
574 pidep(i,k) = 0.
575 pcond(i,k) = 0.
576 psmlt(i,k) = 0.
577 psevp(i,k) = 0.
578 falk(i,k,1) = 0.
579 falk(i,k,2) = 0.
580 fall(i,k,1) = 0.
581 fall(i,k,2) = 0.
582 fallc(i,k) = 0.
583 falkc(i,k) = 0.
584 xni(i,k) = 1.e3
585 enddo
586 !
587 !----------------------------------------------------------------
588 ! compute the fallout term:
589 ! first, vertical terminal velosity for minor loops
590 !
591 do k = kts, kte
592 supcol = t0c-t(i,k,j)
593 !---------------------------------------------------------------
594 ! n0s: Intercept parameter for snow [m-4] [HDC 6]
595 !---------------------------------------------------------------
596 n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
597 if(qrs(i,k,1).le.qcrmin)then
598 rslope(i,k,1) = rslopermax
599 rslopeb(i,k,1) = rsloperbmax
600 rslope2(i,k,1) = rsloper2max
601 rslope3(i,k,1) = rsloper3max
602 else
603 rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k,j))
604 rslopeb(i,k,1) = exp(log(rslope(i,k,1))*(bvtr))
605 rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
606 rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
607 endif
608 if(qrs(i,k,2).le.qcrmin)then
609 rslope(i,k,2) = rslopesmax
610 rslopeb(i,k,2) = rslopesbmax
611 rslope2(i,k,2) = rslopes2max
612 rslope3(i,k,2) = rslopes3max
613 else
614 rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k,j),n0sfac(i,k))
615 rslopeb(i,k,2) = exp(log(rslope(i,k,2))*(bvts))
616 rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
617 rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
618 endif
619 !-------------------------------------------------------------
620 ! Ni: ice crystal number concentraiton [HDC 5c]
621 !-------------------------------------------------------------
622 ! xni(i,k) = min(max(5.38e7*(den(i,k,j) &
623 ! *max(qci(i,k,2),qmin))**0.75,1.e3),1.e6)
624 temp = (den(i,k,j)*max(qci(i,k,2),qmin))
625 temp = sqrt(sqrt(temp*temp*temp))
626 xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
627 enddo
628 !
629 numdt = 1
630 do k = kte, kts, -1
631 work1(i,k,1) = pvtr*rslopeb(i,k,1)*denfac(i,k)/delz(i,k,j)
632 work1(i,k,2) = pvts*rslopeb(i,k,2)*denfac(i,k)/delz(i,k,j)
633 numdt = max(nint(max(work1(i,k,1),work1(i,k,2))*dtcld+.5),1)
634 if(numdt.ge.mstep) mstep = numdt
635 enddo
636 rmstep = 1./mstep
637 !
638 do n = 1, mstep
639 k = kte
640 ! falk(i,k,1) = den(i,k,j)*qrs(i,k,1)*work1(i,k,1)/mstep(i)
641 ! falk(i,k,2) = den(i,k,j)*qrs(i,k,2)*work1(i,k,2)/mstep(i)
642 falk(i,k,1) = den(i,k,j)*qrs(i,k,1)*work1(i,k,1)*rmstep
643 falk(i,k,2) = den(i,k,j)*qrs(i,k,2)*work1(i,k,2)*rmstep
644 fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
645 fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
646 ! qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcld/den(i,k,j),0.)
647 ! qrs(i,k,2) = max(qrs(i,k,2)-falk(i,k,2)*dtcld/den(i,k,j),0.)
648 dtcldden = dtcld/den(i,k,j)
649 qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcldden,0.)
650 qrs(i,k,2) = max(qrs(i,k,2)-falk(i,k,2)*dtcldden,0.)
651 ! endif
652 do k = kte-1, kts, -1
653 falk(i,k,1) = den(i,k,j)*qrs(i,k,1)*work1(i,k,1)*rmstep
654 falk(i,k,2) = den(i,k,j)*qrs(i,k,2)*work1(i,k,2)*rmstep
655 fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
656 fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
657 dtcldden = dtcld/den(i,k,j)
658 rdelz = 1./delz(i,k,j)
659 qrs(i,k,1) = max(qrs(i,k,1)-(falk(i,k,1)-falk(i,k+1,1) &
660 *delz(i,k+1,j)*rdelz)*dtcldden,0.)
661 qrs(i,k,2) = max(qrs(i,k,2)-(falk(i,k,2)-falk(i,k+1,2) &
662 *delz(i,k+1,j)*rdelz)*dtcldden,0.)
663 enddo
664 do k = kte, kts, -1
665 if(t(i,k,j).gt.t0c.and.qrs(i,k,2).gt.0.) then
666 !----------------------------------------------------------------
667 ! psmlt: melting of snow [HL A33] [RH83 A25]
668 ! (T>T0: S->R)
669 !----------------------------------------------------------------
670 xlf = xlf0
671 ! work2(i,k)= venfac(p(i,k),t(i,k,j),den(i,k,j))
672 work2(i,k)= (exp(log(((1.496e-6*((t(i,k,j))*sqrt(t(i,k,j))) &
673 /((t(i,k,j))+120.)/(den(i,k,j)))/(8.794e-5 &
674 *exp(log(t(i,k,j))*(1.81))/p(i,k,j)))) &
675 *((.3333333)))/sqrt((1.496e-6*((t(i,k,j)) &
676 *sqrt(t(i,k,j)))/((t(i,k,j))+120.)/(den(i,k,j)))) &
677 *sqrt(sqrt(den0/(den(i,k,j)))))
678 coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
679 ! psmlt(i,k) = xka(t(i,k,j),den(i,k,j))/xlf*(t0c-t(i,k,j))*pi/2. &
680 ! *n0sfac(i,k)*(precs1*rslope2(i,k,2)+precs2 &
681 ! *work2(i,k)*coeres)
682 psmlt(i,k) = (1.414e3*(1.496e-6*((t(i,k,j))*sqrt(t(i,k,j))) &
683 /((t(i,k,j))+120.)/(den(i,k,j)) )*(den(i,k,j))) &
684 /xlf*(t0c-t(i,k,j))*pi/2. &
685 *n0sfac(i,k)*(precs1*rslope2(i,k,2)+precs2 &
686 *work2(i,k)*coeres)
687 psmlt(i,k) = min(max(psmlt(i,k)*dtcld/mstep, &
688 -qrs(i,k,2)/mstep),0.)
689 qrs(i,k,2) = qrs(i,k,2) + psmlt(i,k)
690 qrs(i,k,1) = qrs(i,k,1) - psmlt(i,k)
691 t(i,k,j) = t(i,k,j) + xlf/cpm(i,k)*psmlt(i,k)
692 endif
693 enddo
694 enddo
697 !---------------------------------------------------------------
698 ! Vice [ms-1] : fallout of ice crystal [HDC 5a]
699 !---------------------------------------------------------------
700 mstep = 1
701 numdt = 1
702 do k = kte, kts, -1
703 if(qci(i,k,2).le.0.) then
704 work2c(i,k) = 0.
705 else
706 xmi = den(i,k,j)*qci(i,k,2)/xni(i,k)
707 ! diameter = min(dicon * sqrt(xmi),dimax)
708 diameter = max(min(dicon * sqrt(xmi),dimax), 1.e-25)
709 work1c(i,k) = 1.49e4*exp(log(diameter)*(1.31))
710 work2c(i,k) = work1c(i,k)/delz(i,k,j)
711 endif
712 numdt = max(nint(work2c(i,k)*dtcld+.5),1)
713 if(numdt.ge.mstep) mstep = numdt
714 enddo
715 !
716 do n = 1, mstep
717 k = kte
718 falkc(i,k) = den(i,k,j)*qci(i,k,2)*work2c(i,k)/mstep
719 holdc = falkc(i,k)
720 fallc(i,k) = fallc(i,k)+falkc(i,k)
721 holdci = qci(i,k,2)
722 qci(i,k,2) = max(qci(i,k,2)-falkc(i,k)*dtcld/den(i,k,j),0.)
723 ! endif
724 do k = kte-1, kts, -1
725 falkc(i,k) = den(i,k,j)*qci(i,k,2)*work2c(i,k)/mstep
726 holdc = falkc(i,k)
727 fallc(i,k) = fallc(i,k)+falkc(i,k)
728 holdci = qci(i,k,2)
729 qci(i,k,2) = max(qci(i,k,2)-(falkc(i,k)-falkc(i,k+1) &
730 *delz(i,k+1,j)/delz(i,k,j))*dtcld/den(i,k,j),0.)
731 ! endif
732 enddo
733 enddo
734 !
735 !
736 !----------------------------------------------------------------
737 ! rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
738 !
739 fallsum = fall(i,1,1)+fall(i,1,2)+fallc(i,1)
740 fallsum_qsi = fall(i,1,2)+fallc(i,1)
741 rainncv(i,j) = 0.
742 if(fallsum.gt.0.) then
743 rainncv(i,j) = fallsum*delz(i,1,j)/denr*dtcld*1000.
744 rain(i,j) = fallsum*delz(i,1,j)/denr*dtcld*1000. + rain(i,j)
745 endif
746 snowncv(i,j) = 0.
747 if(fallsum_qsi.gt.0.) then
748 snowncv(i,j) = fallsum_qsi*delz(i,kts,j)/denr*dtcld*1000.
749 snow(i,j) = fallsum_qsi*delz(i,kts,j)/denr*dtcld*1000. + snow(i,j)
750 endif
751 sr(i,j) = 0.
752 if(fallsum.gt.0.)sr(i,j)=fallsum_qsi*delz(i,kts,j)/denr*dtcld*1000. &
753 /(rainncv(i,j)+1.e-12)
754 !
755 !---------------------------------------------------------------
756 ! pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
757 ! (T>T0: I->C)
758 !---------------------------------------------------------------
759 do k = kts, kte
760 supcol = t0c-t(i,k,j)
761 xlf = xls-xl(i,k)
762 if(supcol.lt.0.) xlf = xlf0
763 if(supcol.lt.0.and.qci(i,k,2).gt.0.) then
764 qci(i,k,1) = qci(i,k,1) + qci(i,k,2)
765 t(i,k,j) = t(i,k,j) - xlf/cpm(i,k)*qci(i,k,2)
766 qci(i,k,2) = 0.
767 endif
768 !---------------------------------------------------------------
769 ! pihmf: homogeneous freezing of cloud water below -40c [HL A45]
770 ! (T<-40C: C->I)
771 !---------------------------------------------------------------
772 if(supcol.gt.40..and.qci(i,k,1).gt.0.) then
773 qci(i,k,2) = qci(i,k,2) + qci(i,k,1)
774 t(i,k,j) = t(i,k,j) + xlf/cpm(i,k)*qci(i,k,1)
775 qci(i,k,1) = 0.
776 endif
777 !---------------------------------------------------------------
778 ! pihtf: heterogeneous freezing of cloud water [HL A44]
779 ! (T0>T>-40C: C->I)
780 !---------------------------------------------------------------
781 if(supcol.gt.0..and.qci(i,k,1).gt.0.) then
782 supcolt=min(supcol,50.)
783 ! pfrzdtc = min(pfrz1*(exp(pfrz2*supcol)-1.) &
784 ! *den(i,k,j)/denr/xncr*qci(i,k,1)**2*dtcld,qci(i,k,1))
785 pfrzdtc = min(pfrz1*(exp(pfrz2*supcolt)-1.) &
786 *den(i,k,j)/denr/xncr*qci(i,k,1)*qci(i,k,1)*dtcld,qci(i,k,1))
787 qci(i,k,2) = qci(i,k,2) + pfrzdtc
788 t(i,k,j) = t(i,k,j) + xlf/cpm(i,k)*pfrzdtc
789 qci(i,k,1) = qci(i,k,1)-pfrzdtc
790 endif
791 !---------------------------------------------------------------
792 ! psfrz: freezing of rain water [HL A20] [LFO 45]
793 ! (T<T0, R->S)
794 !---------------------------------------------------------------
795 if(supcol.gt.0..and.qrs(i,k,1).gt.0.) then
796 supcolt=min(supcol,50.)
797 ! pfrzdtr = min(20.*pi**2*pfrz1*n0r*denr/den(i,k,j) &
798 ! *(exp(pfrz2*supcol)-1.)*rslope(i,k,1)**7*dtcld, &
799 ! qrs(i,k,1))
800 temp = rslope(i,k,1)
801 temp = temp*temp*temp*temp*temp*temp*temp
802 pfrzdtr = min(20.*(pi*pi)*pfrz1*n0r*denr/den(i,k,j) &
803 *(exp(pfrz2*supcolt)-1.)*temp*dtcld, &
804 qrs(i,k,1))
805 qrs(i,k,2) = qrs(i,k,2) + pfrzdtr
806 t(i,k,j) = t(i,k,j) + xlf/cpm(i,k)*pfrzdtr
807 qrs(i,k,1) = qrs(i,k,1)-pfrzdtr
808 endif
809 enddo
810 !
811 !----------------------------------------------------------------
812 ! rsloper: reverse of the slope parameter of the rain(m)
813 ! xka: thermal conductivity of air(jm-1s-1k-1)
814 ! work1: the thermodynamic term in the denominator associated with
815 ! heat conduction and vapor diffusion
816 ! (ry88, y93, h85)
817 ! work2: parameter associated with the ventilation effects(y93)
818 !
819 do k = kts, kte
820 if(qrs(i,k,1).le.qcrmin)then
821 rslope(i,k,1) = rslopermax
822 rslopeb(i,k,1) = rsloperbmax
823 rslope2(i,k,1) = rsloper2max
824 rslope3(i,k,1) = rsloper3max
825 else
826 ! rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k,j))
827 rslope(i,k,1) = 1./(sqrt(sqrt(pidn0r/((qrs(i,k,1))*(den(i,k,j))))))
828 rslopeb(i,k,1) = exp(log(rslope(i,k,1))*(bvtr))
829 rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
830 rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
831 endif
832 if(qrs(i,k,2).le.qcrmin)then
833 rslope(i,k,2) = rslopesmax
834 rslopeb(i,k,2) = rslopesbmax
835 rslope2(i,k,2) = rslopes2max
836 rslope3(i,k,2) = rslopes3max
837 else
838 ! rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k,j),n0sfac(i,k))
839 rslope(i,k,2) = 1./(sqrt(sqrt(pidn0s*(n0sfac(i,k))/((qrs(i,k,2)) &
840 *(den(i,k,j))))))
841 rslopeb(i,k,2) = exp(log(rslope(i,k,2))*(bvts))
842 rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
843 rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
844 endif
845 enddo
846 !
847 do k = kts, kte
848 ! work1(i,k,1) = diffac(xl(i,k),p(i,k,j),t(i,k,j),den(i,k,j),qs(i,k,1))
849 work1(i,k,1) = ((((den(i,k,j))*(xl(i,k))*(xl(i,k)))*((t(i,k,j))+120.) &
850 *(den(i,k,j)))/(1.414e3*(1.496e-6*((t(i,k,j))*sqrt(t(i,k,j))))&
851 *(den(i,k,j))*(rv*(t(i,k,j))*(t(i,k,j))))) &
852 + p(i,k,j)/((qs(i,k,1))*(8.794e-5*exp(log(t(i,k,j))*(1.81))))
853 ! work1(i,k,2) = diffac(xls,p(i,k,j),t(i,k,j),den(i,k,j),qs(i,k,2))
854 work1(i,k,2) = ((((den(i,k,j))*(xls)*(xls))*((t(i,k,j))+120.)*(den(i,k,j)))&
855 /(1.414e3*(1.496e-6*((t(i,k,j))*sqrt(t(i,k,j))))*(den(i,k,j)) &
856 *(rv*(t(i,k,j))*(t(i,k,j)))) &
857 + p(i,k,j)/(qs(i,k,2)*(8.794e-5*exp(log(t(i,k,j))*(1.81)))))
858 ! work2(i,k) = venfac(p(i,k,j),t(i,k,j),den(i,k,j))
859 work2(i,k) = (exp(.3333333*log(((1.496e-6 * ((t(i,k,j))*sqrt(t(i,k,j)))) &
860 *p(i,k,j))/(((t(i,k,j))+120.)*den(i,k,j)*(8.794e-5 &
861 *exp(log(t(i,k,j))*(1.81))))))*sqrt(sqrt(den0/(den(i,k,j))))) &
862 /sqrt((1.496e-6*((t(i,k,j))*sqrt(t(i,k,j)))) &
863 /(((t(i,k,j))+120.)*den(i,k,j)))
864 enddo
865 !
866 !===============================================================
867 !
868 ! warm rain processes
869 !
870 ! - follows the processes in RH83 and LFO except for autoconcersion
871 !
872 !===============================================================
873 !
874 do k = kts, kte
875 supsat = max(q(i,k,j),qmin)-qs(i,k,1)
876 satdt = supsat/dtcld
877 !---------------------------------------------------------------
878 ! praut: auto conversion rate from cloud to rain [HDC 16]
879 ! (C->R)
880 !---------------------------------------------------------------
881 if(qci(i,k,1).gt.qc0) then
882 praut(i,k) = qck1*exp(log(qci(i,k,1))*((7./3.)))
883 praut(i,k) = min(praut(i,k),qci(i,k,1)/dtcld)
884 endif
885 !---------------------------------------------------------------
886 ! pracw: accretion of cloud water by rain [HL A40] [LFO 51]
887 ! (C->R)
888 !---------------------------------------------------------------
889 if(qrs(i,k,1).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
890 pracw(i,k) = min(pacrr*rslope3(i,k,1)*rslopeb(i,k,1) &
891 *qci(i,k,1)*denfac(i,k),qci(i,k,1)/dtcld)
892 endif
893 !---------------------------------------------------------------
894 ! prevp: evaporation/condensation rate of rain [HDC 14]
895 ! (V->R or R->V)
896 !---------------------------------------------------------------
897 if(qrs(i,k,1).gt.0.) then
898 coeres = rslope2(i,k,1)*sqrt(rslope(i,k,1)*rslopeb(i,k,1))
899 prevp(i,k) = (rh(i,k,1)-1.)*(precr1*rslope2(i,k,1) &
900 +precr2*work2(i,k)*coeres)/work1(i,k,1)
901 if(prevp(i,k).lt.0.) then
902 prevp(i,k) = max(prevp(i,k),-qrs(i,k,1)/dtcld)
903 prevp(i,k) = max(prevp(i,k),satdt/2)
904 else
905 prevp(i,k) = min(prevp(i,k),satdt/2)
906 endif
907 endif
908 enddo
909 !
910 !===============================================================
911 !
912 ! cold rain processes
913 !
914 ! - follows the revised ice microphysics processes in HDC
915 ! - the processes same as in RH83 and RH84 and LFO behave
916 ! following ice crystal hapits defined in HDC, inclduing
917 ! intercept parameter for snow (n0s), ice crystal number
918 ! concentration (ni), ice nuclei number concentration
919 ! (n0i), ice diameter (d)
920 !
921 !===============================================================
922 !
923 rdtcld = 1./dtcld
924 do k = kts, kte
925 supcol = t0c-t(i,k,j)
926 supsat = max(q(i,k,j),qmin)-qs(i,k,2)
927 satdt = supsat/dtcld
928 ifsat = 0
929 !-------------------------------------------------------------
930 ! Ni: ice crystal number concentraiton [HDC 5c]
931 !-------------------------------------------------------------
932 ! xni(i,k) = min(max(5.38e7*(den(i,k,j) &
933 ! *max(qci(i,k,2),qmin))**0.75,1.e3),1.e6)
934 temp = (den(i,k,j)*max(qci(i,k,2),qmin))
935 temp = sqrt(sqrt(temp*temp*temp))
936 xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
937 eacrs = exp(0.07*(-supcol))
938 !
939 if(supcol.gt.0) then
940 if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,2).gt.qmin) then
941 xmi = den(i,k,j)*qci(i,k,2)/xni(i,k)
942 diameter = min(dicon * sqrt(xmi),dimax)
943 vt2i = 1.49e4*diameter**1.31
944 vt2s = pvts*rslopeb(i,k,2)*denfac(i,k)
945 !-------------------------------------------------------------
946 ! psaci: Accretion of cloud ice by rain [HDC 10]
947 ! (T<T0: I->S)
948 !-------------------------------------------------------------
949 acrfac = 2.*rslope3(i,k,2)+2.*diameter*rslope2(i,k,2) &
950 +diameter**2*rslope(i,k,2)
951 psaci(i,k) = pi*qci(i,k,2)*eacrs*n0s*n0sfac(i,k) &
952 *abs(vt2s-vt2i)*acrfac/4.
953 endif
954 endif
955 !-------------------------------------------------------------
956 ! psacw: Accretion of cloud water by snow [HL A7] [LFO 24]
957 ! (T<T0: C->S, and T>=T0: C->R)
958 !-------------------------------------------------------------
959 if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
960 psacw(i,k) = min(pacrc*n0sfac(i,k)*rslope3(i,k,2) &
961 *rslopeb(i,k,2)*qci(i,k,1)*denfac(i,k) &
962 ! ,qci(i,k,1)/dtcld)
963 ,qci(i,k,1)*rdtcld)
964 endif
965 if(supcol .gt. 0) then
966 !-------------------------------------------------------------
967 ! pidep: Deposition/Sublimation rate of ice [HDC 9]
968 ! (T<T0: V->I or I->V)
969 !-------------------------------------------------------------
970 if(qci(i,k,2).gt.0.and.ifsat.ne.1) then
971 xmi = den(i,k,j)*qci(i,k,2)/xni(i,k)
972 diameter = dicon * sqrt(xmi)
973 pidep(i,k) = 4.*diameter*xni(i,k)*(rh(i,k,2)-1.)/work1(i,k,2)
974 supice = satdt-prevp(i,k)
975 if(pidep(i,k).lt.0.) then
976 ! pidep(i,k) = max(max(pidep(i,k),satdt/2),supice)
977 ! pidep(i,k) = max(pidep(i,k),-qci(i,k,2)/dtcld)
978 pidep(i,k) = max(max(pidep(i,k),satdt*.5),supice)
979 pidep(i,k) = max(pidep(i,k),-qci(i,k,2)*rdtcld)
980 else
981 ! pidep(i,k) = min(min(pidep(i,k),satdt/2),supice)
982 pidep(i,k) = min(min(pidep(i,k),satdt*.5),supice)
983 endif
984 if(abs(prevp(i,k)+pidep(i,k)).ge.abs(satdt)) ifsat = 1
985 endif
986 !-------------------------------------------------------------
987 ! psdep: deposition/sublimation rate of snow [HDC 14]
988 ! (V->S or S->V)
989 !-------------------------------------------------------------
990 if(qrs(i,k,2).gt.0..and.ifsat.ne.1) then
991 coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
992 psdep(i,k) = (rh(i,k,2)-1.)*n0sfac(i,k) &
993 *(precs1*rslope2(i,k,2)+precs2 &
994 *work2(i,k)*coeres)/work1(i,k,2)
995 supice = satdt-prevp(i,k)-pidep(i,k)
996 if(psdep(i,k).lt.0.) then
997 ! psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)/dtcld)
998 ! psdep(i,k) = max(max(psdep(i,k),satdt/2),supice)
999 psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)*rdtcld)
1000 psdep(i,k) = max(max(psdep(i,k),satdt*.5),supice)
1001 else
1002 ! psdep(i,k) = min(min(psdep(i,k),satdt/2),supice)
1003 psdep(i,k) = min(min(psdep(i,k),satdt*.5),supice)
1004 endif
1005 if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)).ge.abs(satdt)) &
1006 ifsat = 1
1007 endif
1008 !-------------------------------------------------------------
1009 ! pigen: generation(nucleation) of ice from vapor [HL A50] [HDC 7-8]
1010 ! (T<T0: V->I)
1011 !-------------------------------------------------------------
1012 if(supsat.gt.0.and.ifsat.ne.1) then
1013 supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)
1014 xni0 = 1.e3*exp(0.1*supcol)
1015 roqi0 = 4.92e-11*exp(log(xni0)*(1.33))
1016 pigen(i,k) = max(0.,(roqi0/den(i,k,j)-max(qci(i,k,2),0.)) &
1017 ! /dtcld)
1018 *rdtcld)
1019 pigen(i,k) = min(min(pigen(i,k),satdt),supice)
1020 endif
1021 !
1022 !-------------------------------------------------------------
1023 ! psaut: conversion(aggregation) of ice to snow [HDC 12]
1024 ! (T<T0: I->S)
1025 !-------------------------------------------------------------
1026 if(qci(i,k,2).gt.0.) then
1027 qimax = roqimax/den(i,k,j)
1028 ! psaut(i,k) = max(0.,(qci(i,k,2)-qimax)/dtcld)
1029 psaut(i,k) = max(0.,(qci(i,k,2)-qimax)*rdtcld)
1030 endif
1031 endif
1032 !-------------------------------------------------------------
1033 ! psevp: Evaporation of melting snow [HL A35] [RH83 A27]
1034 ! (T>T0: S->V)
1035 !-------------------------------------------------------------
1036 if(supcol.lt.0.) then
1037 if(qrs(i,k,2).gt.0..and.rh(i,k,1).lt.1.) &
1038 psevp(i,k) = psdep(i,k)*work1(i,k,2)/work1(i,k,1)
1039 ! psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)/dtcld),0.)
1040 psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)*rdtcld),0.)
1041 endif
1042 enddo
1043 !
1044 !
1045 !----------------------------------------------------------------
1046 ! check mass conservation of generation terms and feedback to the
1047 ! large scale
1048 !
1049 do k = kts, kte
1050 if(t(i,k,j).le.t0c) then
1051 !
1052 ! cloud water
1053 !
1054 value = max(qmin,qci(i,k,1))
1055 source = (praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
1056 if (source.gt.value) then
1057 factor = value/source
1058 praut(i,k) = praut(i,k)*factor
1059 pracw(i,k) = pracw(i,k)*factor
1060 psacw(i,k) = psacw(i,k)*factor
1061 endif
1062 !
1063 ! cloud ice
1064 !
1065 value = max(qmin,qci(i,k,2))
1066 source = (psaut(i,k)+psaci(i,k)-pigen(i,k)-pidep(i,k))*dtcld
1067 if (source.gt.value) then
1068 factor = value/source
1069 psaut(i,k) = psaut(i,k)*factor
1070 psaci(i,k) = psaci(i,k)*factor
1071 pigen(i,k) = pigen(i,k)*factor
1072 pidep(i,k) = pidep(i,k)*factor
1073 endif
1074 !
1075 ! rain
1076 !
1077 !
1078 value = max(qmin,qrs(i,k,1))
1079 source = (-praut(i,k)-pracw(i,k)-prevp(i,k))*dtcld
1080 if (source.gt.value) then
1081 factor = value/source
1082 praut(i,k) = praut(i,k)*factor
1083 pracw(i,k) = pracw(i,k)*factor
1084 prevp(i,k) = prevp(i,k)*factor
1085 endif
1086 !
1087 ! snow
1088 !
1089 value = max(qmin,qrs(i,k,2))
1090 source = (-psdep(i,k)-psaut(i,k)-psaci(i,k)-psacw(i,k))*dtcld
1091 if (source.gt.value) then
1092 factor = value/source
1093 psdep(i,k) = psdep(i,k)*factor
1094 psaut(i,k) = psaut(i,k)*factor
1095 psaci(i,k) = psaci(i,k)*factor
1096 psacw(i,k) = psacw(i,k)*factor
1097 endif
1098 !
1099 work2(i,k)=-(prevp(i,k)+psdep(i,k)+pigen(i,k)+pidep(i,k))
1100 ! update
1101 q(i,k,j) = q(i,k,j)+work2(i,k)*dtcld
1102 qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
1103 +psacw(i,k))*dtcld,0.)
1104 qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
1105 +prevp(i,k))*dtcld,0.)
1106 qci(i,k,2) = max(qci(i,k,2)-(psaut(i,k)+psaci(i,k) &
1107 -pigen(i,k)-pidep(i,k))*dtcld,0.)
1108 qrs(i,k,2) = max(qrs(i,k,2)+(psdep(i,k)+psaut(i,k) &
1109 +psaci(i,k)+psacw(i,k))*dtcld,0.)
1110 xlf = xls-xl(i,k)
1111 xlwork2 = -xls*(psdep(i,k)+pidep(i,k)+pigen(i,k)) &
1112 -xl(i,k)*prevp(i,k)-xlf*psacw(i,k)
1113 t(i,k,j) = t(i,k,j)-xlwork2/cpm(i,k)*dtcld
1114 else
1115 !
1116 ! cloud water
1117 !
1118 value = max(qmin,qci(i,k,1))
1119 source=(praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
1120 if (source.gt.value) then
1121 factor = value/source
1122 praut(i,k) = praut(i,k)*factor
1123 pracw(i,k) = pracw(i,k)*factor
1124 psacw(i,k) = psacw(i,k)*factor
1125 endif
1126 !
1127 ! rain
1128 !
1129 value = max(qmin,qrs(i,k,1))
1130 source = (-praut(i,k)-pracw(i,k)-prevp(i,k)-psacw(i,k))*dtcld
1131 if (source.gt.value) then
1132 factor = value/source
1133 praut(i,k) = praut(i,k)*factor
1134 pracw(i,k) = pracw(i,k)*factor
1135 prevp(i,k) = prevp(i,k)*factor
1136 psacw(i,k) = psacw(i,k)*factor
1137 endif
1138 !
1139 ! snow
1140 !
1141 value = max(qcrmin,qrs(i,k,2))
1142 source=(-psevp(i,k))*dtcld
1143 if (source.gt.value) then
1144 factor = value/source
1145 psevp(i,k) = psevp(i,k)*factor
1146 endif
1147 work2(i,k)=-(prevp(i,k)+psevp(i,k))
1148 ! update
1149 q(i,k,j) = q(i,k,j)+work2(i,k)*dtcld
1150 qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
1151 +psacw(i,k))*dtcld,0.)
1152 qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
1153 +prevp(i,k) +psacw(i,k))*dtcld,0.)
1154 qrs(i,k,2) = max(qrs(i,k,2)+psevp(i,k)*dtcld,0.)
1155 xlf = xls-xl(i,k)
1156 xlwork2 = -xl(i,k)*(prevp(i,k)+psevp(i,k))
1157 t(i,k,j) = t(i,k,j)-xlwork2/cpm(i,k)*dtcld
1158 endif
1159 enddo
1160 !
1161 ! Inline expansion for fpvs
1162 ! qs(i,k,1) = fpvs(t(i,k,j),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1163 ! qs(i,k,2) = fpvs(t(i,k,j),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1164 hsub = xls
1165 hvap = xlv0
1166 cvap = cpv
1167 ttp=t0c+0.01
1168 dldt=cvap-cliq
1169 xa=-dldt/rv
1170 xb=xa+hvap/(rv*ttp)
1171 dldti=cvap-cice
1172 xai=-dldti/rv
1173 xbi=xai+hsub/(rv*ttp)
1174 do k = kts, kte
1175 tr=ttp/t(i,k,j)
1176 logtr = log(tr)
1177 qs(i,k,1)=psat*exp(logtr*(xa)+xb*(1.-tr))
1178 qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k,j) - qs(i,k,1))
1179 qs(i,k,1) = max(qs(i,k,1),qmin)
1180 enddo
1181 !
1182 !----------------------------------------------------------------
1183 ! pcond: condensational/evaporational rate of cloud water [HL A46] [RH83 A6]
1184 ! if there exists additional water vapor condensated/if
1185 ! evaporation of cloud water is not enough to remove subsaturation
1186 !
1187 do k = kts, kte
1188 ! work1(i,k,1) = conden(t(i,k,j),q(i,k,j),qs(i,k,1),xl(i,k),cpm(i,k))
1189 work1(i,k,1) = ((max(q(i,k,j),qmin)-(qs(i,k,1)))/(1.+(xl(i,k)) &
1190 *(xl(i,k))/(rv*(cpm(i,k)))*(qs(i,k,1)) &
1191 /((t(i,k,j))*(t(i,k,j)))))
1192 work2(i,k) = qci(i,k,1)+work1(i,k,1)
1193 pcond(i,k) = min(max(work1(i,k,1)/dtcld,0.),max(q(i,k,j),0.)/dtcld)
1194 if(qci(i,k,1).gt.0..and.work1(i,k,1).lt.0.) &
1195 pcond(i,k) = max(work1(i,k,1),-qci(i,k,1))/dtcld
1196 q(i,k,j) = q(i,k,j)-pcond(i,k)*dtcld
1197 qci(i,k,1) = max(qci(i,k,1)+pcond(i,k)*dtcld,0.)
1198 t(i,k,j) = t(i,k,j)+pcond(i,k)*xl(i,k)/cpm(i,k)*dtcld
1199 enddo
1200 !
1201 !
1202 !----------------------------------------------------------------
1203 ! padding for small values
1204 !
1205 do k = kts, kte
1206 if(qci(i,k,1).le.qmin) qci(i,k,1) = 0.0
1207 if(qci(i,k,2).le.qmin) qci(i,k,2) = 0.0
1208 enddo
1209 enddo ! big loops
1211 DO K=kts,kte
1212 th(i,k,j)=t(i,k,j)/pii(i,k,j)
1213 qc(i,k,j) = qci(i,k,1)
1214 qi(i,k,j) = qci(i,k,2)
1215 qr(i,k,j) = qrs(i,k,1)
1216 qqs(i,k,j) = qrs(i,k,2)
1217 ENDDO
1221 ENDDO ! i loop
1222 enddo ! j loop
1223 !$acc end region
1225 ELSE
1227 !
1228 ! Moved outside of accelerator region
1229 !
1230 loops = max(nint(delt/dtcldcr),1)
1231 dtcld = delt/loops
1232 if(delt.le.dtcldcr) dtcld = delt
1234 !$acc region &
1235 !$acc local(t) &
1236 !$acc copyin(delz(:,:,:),p(:,:,:),den(:,:,:),pii(:,:,:)) &
1237 !$acc copyout(rainncv(:,:),sr(:,:)) &
1238 !$acc copy(qqs(:,:,:),qr(:,:,:),qi(:,:,:),qc(:,:,:)) &
1239 !$acc copy(th(:,:,:),q(:,:,:),rain(:,:))
1240 !$acc do &
1241 !$acc private(rh,qs,rslope,rslope2,rslope3,rslopeb,falk,fall) &
1242 !$acc private(work1,qci,qrs,falkc,fallc,xl,cpm,denfac,xni) &
1243 !$acc private(n0sfac,work2,work1c,work2c,pigen,pidep,psdep) &
1244 !$acc private(praut,psaut,prevp,psevp) &
1245 !$acc private(pracw,psacw,psaci,pcond,psmlt) &
1246 !$acc parallel
1247 do j = jts, jte
1248 !$acc do &
1249 !$acc private(numdt,mstep) &
1250 !$acc kernel vector
1251 do i = its, ite
1252 do k = kts, kte
1253 t(i,k,j)=th(i,k,j)*pii(i,k,j)
1254 qci(i,k,1) = max(qc(i,k,j),0.0)
1255 qci(i,k,2) = max(qi(i,k,j),0.0)
1256 qrs(i,k,1) = max(qr(i,k,j),0.0)
1257 qrs(i,k,2) = max(qqs(i,k,j),0.0)
1258 enddo
1259 !
1260 !----------------------------------------------------------------
1261 ! latent heat for phase changes and heat capacity. neglect the
1262 ! changes during microphysical process calculation
1263 ! emanuel(1994)
1264 !
1265 do k = kts, kte
1266 cpm(i,k) = cpmcal(q(i,k,j))
1267 xl(i,k) = xlcal(t(i,k,j))
1268 enddo
1269 !
1270 !----------------------------------------------------------------
1271 ! compute the minor time steps.
1272 !
1273 ! loops = max(nint(delt/dtcldcr),1)
1274 ! dtcld = delt/loops
1275 ! if(delt.le.dtcldcr) dtcld = delt
1276 !
1277 do loop = 1,loops
1278 !
1279 !----------------------------------------------------------------
1280 ! initialize the large scale variables
1281 !
1282 mstep = 1
1283 flgcld = .true.
1284 !
1285 do k = kts, kte
1286 denfac(i,k) = sqrt(den0/den(i,k,j))
1287 enddo
1288 ! do k = kts, kte
1289 ! CALL VREC( tvec1(its), den(its,k,j), ite-its+1)
1290 ! do i = its, ite
1291 ! tvec1(i) = tvec1(i)*den0
1292 ! enddo
1293 ! CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
1294 ! enddo
1295 !
1296 ! Inline expansion for fpvs
1297 ! qs(i,k,1) = fpvs(t(i,k,j),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1298 ! qs(i,k,2) = fpvs(t(i,k,j),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1299 hsub = xls
1300 hvap = xlv0
1301 cvap = cpv
1302 ttp=t0c+0.01
1303 dldt=cvap-cliq
1304 xa=-dldt/rv
1305 xb=xa+hvap/(rv*ttp)
1306 dldti=cvap-cice
1307 xai=-dldti/rv
1308 xbi=xai+hsub/(rv*ttp)
1310 ! this is for compilers where the conditional inhibits vectorization
1311 #ifdef WSM_NO_CONDITIONAL_IN_VECTOR
1312 do k = kts, kte
1313 if(t(i,k,j).lt.ttp) then
1314 xal = xai
1315 xbl = xbi
1316 else
1317 xal = xa
1318 xbl = xb
1319 endif
1320 tr=ttp/t(i,k,j)
1321 logtr=log(tr)
1322 qs(i,k,1)=psat*exp(logtr*(xa)+xb*(1.-tr))
1323 qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k,j) - qs(i,k,1))
1324 qs(i,k,1) = max(qs(i,k,1),qmin)
1325 rh(i,k,1) = max(q(i,k,j) / qs(i,k,1),qmin)
1326 qs(i,k,2)=psat*exp(logtr*(xal)+xbl*(1.-tr))
1327 qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k,j) - qs(i,k,2))
1328 qs(i,k,2) = max(qs(i,k,2),qmin)
1329 rh(i,k,2) = max(q(i,k,j) / qs(i,k,2),qmin)
1330 enddo
1331 #else
1332 do k = kts, kte
1333 tr=ttp/t(i,k,j)
1334 logtr=log(tr)
1335 qs(i,k,1)=psat*exp(logtr*(xa)+xb*(1.-tr))
1336 qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k,j) - qs(i,k,1))
1337 qs(i,k,1) = max(qs(i,k,1),qmin)
1338 rh(i,k,1) = max(q(i,k,j) / qs(i,k,1),qmin)
1339 if(t(i,k,j).lt.ttp) then
1340 qs(i,k,2)=psat*exp(logtr*(xai)+xbi*(1.-tr))
1341 else
1342 qs(i,k,2)=psat*exp(logtr*(xa)+xb*(1.-tr))
1343 endif
1344 qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k,j) - qs(i,k,2))
1345 qs(i,k,2) = max(qs(i,k,2),qmin)
1346 rh(i,k,2) = max(q(i,k,j) / qs(i,k,2),qmin)
1347 enddo
1348 #endif
1349 !
1350 !----------------------------------------------------------------
1351 ! initialize the variables for microphysical physics
1352 !
1353 !
1354 do k = kts, kte
1355 prevp(i,k) = 0.
1356 psdep(i,k) = 0.
1357 praut(i,k) = 0.
1358 psaut(i,k) = 0.
1359 pracw(i,k) = 0.
1360 psaci(i,k) = 0.
1361 psacw(i,k) = 0.
1362 pigen(i,k) = 0.
1363 pidep(i,k) = 0.
1364 pcond(i,k) = 0.
1365 psmlt(i,k) = 0.
1366 psevp(i,k) = 0.
1367 falk(i,k,1) = 0.
1368 falk(i,k,2) = 0.
1369 fall(i,k,1) = 0.
1370 fall(i,k,2) = 0.
1371 fallc(i,k) = 0.
1372 falkc(i,k) = 0.
1373 xni(i,k) = 1.e3
1374 enddo
1375 !
1376 !----------------------------------------------------------------
1377 ! compute the fallout term:
1378 ! first, vertical terminal velosity for minor loops
1379 !
1380 do k = kts, kte
1381 supcol = t0c-t(i,k,j)
1382 !---------------------------------------------------------------
1383 ! n0s: Intercept parameter for snow [m-4] [HDC 6]
1384 !---------------------------------------------------------------
1385 n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
1386 if(qrs(i,k,1).le.qcrmin)then
1387 rslope(i,k,1) = rslopermax
1388 rslopeb(i,k,1) = rsloperbmax
1389 rslope2(i,k,1) = rsloper2max
1390 rslope3(i,k,1) = rsloper3max
1391 else
1392 rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k,j))
1393 rslopeb(i,k,1) = exp(log(rslope(i,k,1))*(bvtr))
1394 rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
1395 rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
1396 endif
1397 if(qrs(i,k,2).le.qcrmin)then
1398 rslope(i,k,2) = rslopesmax
1399 rslopeb(i,k,2) = rslopesbmax
1400 rslope2(i,k,2) = rslopes2max
1401 rslope3(i,k,2) = rslopes3max
1402 else
1403 rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k,j),n0sfac(i,k))
1404 rslopeb(i,k,2) = exp(log(rslope(i,k,2))*(bvts))
1405 rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
1406 rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
1407 endif
1408 !-------------------------------------------------------------
1409 ! Ni: ice crystal number concentraiton [HDC 5c]
1410 !-------------------------------------------------------------
1411 ! xni(i,k) = min(max(5.38e7*(den(i,k,j) &
1412 ! *max(qci(i,k,2),qmin))**0.75,1.e3),1.e6)
1413 temp = (den(i,k,j)*max(qci(i,k,2),qmin))
1414 temp = sqrt(sqrt(temp*temp*temp))
1415 xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
1416 enddo
1417 !
1418 numdt = 1
1419 do k = kte, kts, -1
1420 work1(i,k,1) = pvtr*rslopeb(i,k,1)*denfac(i,k)/delz(i,k,j)
1421 work1(i,k,2) = pvts*rslopeb(i,k,2)*denfac(i,k)/delz(i,k,j)
1422 numdt = max(nint(max(work1(i,k,1),work1(i,k,2))*dtcld+.5),1)
1423 if(numdt.ge.mstep) mstep = numdt
1424 enddo
1425 rmstep = 1./mstep
1426 !
1427 do n = 1, mstep
1428 k = kte
1429 ! falk(i,k,1) = den(i,k,j)*qrs(i,k,1)*work1(i,k,1)/mstep(i)
1430 ! falk(i,k,2) = den(i,k,j)*qrs(i,k,2)*work1(i,k,2)/mstep(i)
1431 falk(i,k,1) = den(i,k,j)*qrs(i,k,1)*work1(i,k,1)*rmstep
1432 falk(i,k,2) = den(i,k,j)*qrs(i,k,2)*work1(i,k,2)*rmstep
1433 fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
1434 fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
1435 ! qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcld/den(i,k,j),0.)
1436 ! qrs(i,k,2) = max(qrs(i,k,2)-falk(i,k,2)*dtcld/den(i,k,j),0.)
1437 dtcldden = dtcld/den(i,k,j)
1438 qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcldden,0.)
1439 qrs(i,k,2) = max(qrs(i,k,2)-falk(i,k,2)*dtcldden,0.)
1440 ! endif
1441 do k = kte-1, kts, -1
1442 falk(i,k,1) = den(i,k,j)*qrs(i,k,1)*work1(i,k,1)*rmstep
1443 falk(i,k,2) = den(i,k,j)*qrs(i,k,2)*work1(i,k,2)*rmstep
1444 fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
1445 fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
1446 dtcldden = dtcld/den(i,k,j)
1447 rdelz = 1./delz(i,k,j)
1448 qrs(i,k,1) = max(qrs(i,k,1)-(falk(i,k,1)-falk(i,k+1,1) &
1449 *delz(i,k+1,j)*rdelz)*dtcldden,0.)
1450 qrs(i,k,2) = max(qrs(i,k,2)-(falk(i,k,2)-falk(i,k+1,2) &
1451 *delz(i,k+1,j)*rdelz)*dtcldden,0.)
1452 enddo
1453 do k = kte, kts, -1
1454 if(t(i,k,j).gt.t0c.and.qrs(i,k,2).gt.0.) then
1455 !----------------------------------------------------------------
1456 ! psmlt: melting of snow [HL A33] [RH83 A25]
1457 ! (T>T0: S->R)
1458 !----------------------------------------------------------------
1459 xlf = xlf0
1460 ! work2(i,k)= venfac(p(i,k),t(i,k,j),den(i,k,j))
1461 work2(i,k)= (exp(log(((1.496e-6*((t(i,k,j))*sqrt(t(i,k,j))) &
1462 /((t(i,k,j))+120.)/(den(i,k,j)))/(8.794e-5 &
1463 *exp(log(t(i,k,j))*(1.81))/p(i,k,j)))) &
1464 *((.3333333)))/sqrt((1.496e-6*((t(i,k,j)) &
1465 *sqrt(t(i,k,j)))/((t(i,k,j))+120.)/(den(i,k,j)))) &
1466 *sqrt(sqrt(den0/(den(i,k,j)))))
1467 coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
1468 ! psmlt(i,k) = xka(t(i,k,j),den(i,k,j))/xlf*(t0c-t(i,k,j))*pi/2. &
1469 ! *n0sfac(i,k)*(precs1*rslope2(i,k,2)+precs2 &
1470 ! *work2(i,k)*coeres)
1471 psmlt(i,k) = (1.414e3*(1.496e-6*((t(i,k,j))*sqrt(t(i,k,j))) &
1472 /((t(i,k,j))+120.)/(den(i,k,j)) )*(den(i,k,j))) &
1473 /xlf*(t0c-t(i,k,j))*pi/2. &
1474 *n0sfac(i,k)*(precs1*rslope2(i,k,2)+precs2 &
1475 *work2(i,k)*coeres)
1476 psmlt(i,k) = min(max(psmlt(i,k)*dtcld/mstep, &
1477 -qrs(i,k,2)/mstep),0.)
1478 qrs(i,k,2) = qrs(i,k,2) + psmlt(i,k)
1479 qrs(i,k,1) = qrs(i,k,1) - psmlt(i,k)
1480 t(i,k,j) = t(i,k,j) + xlf/cpm(i,k)*psmlt(i,k)
1481 endif
1482 enddo
1483 enddo
1484 !---------------------------------------------------------------
1485 ! Vice [ms-1] : fallout of ice crystal [HDC 5a]
1486 !---------------------------------------------------------------
1487 mstep = 1
1488 numdt = 1
1489 do k = kte, kts, -1
1490 if(qci(i,k,2).le.0.) then
1491 work2c(i,k) = 0.
1492 else
1493 xmi = den(i,k,j)*qci(i,k,2)/xni(i,k)
1494 ! diameter = min(dicon * sqrt(xmi),dimax)
1495 diameter = max(min(dicon * sqrt(xmi),dimax), 1.e-25)
1496 work1c(i,k) = 1.49e4*exp(log(diameter)*(1.31))
1497 work2c(i,k) = work1c(i,k)/delz(i,k,j)
1498 endif
1499 numdt = max(nint(work2c(i,k)*dtcld+.5),1)
1500 if(numdt.ge.mstep) mstep = numdt
1501 enddo
1502 !
1503 do n = 1, mstep
1504 k = kte
1505 falkc(i,k) = den(i,k,j)*qci(i,k,2)*work2c(i,k)/mstep
1506 holdc = falkc(i,k)
1507 fallc(i,k) = fallc(i,k)+falkc(i,k)
1508 holdci = qci(i,k,2)
1509 qci(i,k,2) = max(qci(i,k,2)-falkc(i,k)*dtcld/den(i,k,j),0.)
1510 do k = kte-1, kts, -1
1511 falkc(i,k) = den(i,k,j)*qci(i,k,2)*work2c(i,k)/mstep
1512 holdc = falkc(i,k)
1513 fallc(i,k) = fallc(i,k)+falkc(i,k)
1514 holdci = qci(i,k,2)
1515 qci(i,k,2) = max(qci(i,k,2)-(falkc(i,k)-falkc(i,k+1) &
1516 *delz(i,k+1,j)/delz(i,k,j))*dtcld/den(i,k,j),0.)
1517 enddo
1518 enddo
1519 !
1520 !
1521 !----------------------------------------------------------------
1522 ! rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
1523 !
1524 fallsum = fall(i,1,1)+fall(i,1,2)+fallc(i,1)
1525 fallsum_qsi = fall(i,1,2)+fallc(i,1)
1526 rainncv(i,j) = 0.
1527 if(fallsum.gt.0.) then
1528 rainncv(i,j) = fallsum*delz(i,1,j)/denr*dtcld*1000.
1529 rain(i,j) = fallsum*delz(i,1,j)/denr*dtcld*1000. + rain(i,j)
1530 endif
1531 sr(i,j) = 0.
1532 if(fallsum.gt.0.)sr(i,j)=fallsum_qsi*delz(i,kts,j)/denr*dtcld*1000. &
1533 /(rainncv(i,j)+1.e-12)
1534 !
1535 !---------------------------------------------------------------
1536 ! pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
1537 ! (T>T0: I->C)
1538 !---------------------------------------------------------------
1539 do k = kts, kte
1540 supcol = t0c-t(i,k,j)
1541 xlf = xls-xl(i,k)
1542 if(supcol.lt.0.) xlf = xlf0
1543 if(supcol.lt.0.and.qci(i,k,2).gt.0.) then
1544 qci(i,k,1) = qci(i,k,1) + qci(i,k,2)
1545 t(i,k,j) = t(i,k,j) - xlf/cpm(i,k)*qci(i,k,2)
1546 qci(i,k,2) = 0.
1547 endif
1548 !---------------------------------------------------------------
1549 ! pihmf: homogeneous freezing of cloud water below -40c [HL A45]
1550 ! (T<-40C: C->I)
1551 !---------------------------------------------------------------
1552 if(supcol.gt.40..and.qci(i,k,1).gt.0.) then
1553 qci(i,k,2) = qci(i,k,2) + qci(i,k,1)
1554 t(i,k,j) = t(i,k,j) + xlf/cpm(i,k)*qci(i,k,1)
1555 qci(i,k,1) = 0.
1556 endif
1557 !---------------------------------------------------------------
1558 ! pihtf: heterogeneous freezing of cloud water [HL A44]
1559 ! (T0>T>-40C: C->I)
1560 !---------------------------------------------------------------
1561 if(supcol.gt.0..and.qci(i,k,1).gt.0.) then
1562 supcolt=min(supcol,50.)
1563 ! pfrzdtc = min(pfrz1*(exp(pfrz2*supcol)-1.) &
1564 ! *den(i,k,j)/denr/xncr*qci(i,k,1)**2*dtcld,qci(i,k,1))
1565 pfrzdtc = min(pfrz1*(exp(pfrz2*supcolt)-1.) &
1566 *den(i,k,j)/denr/xncr*qci(i,k,1)*qci(i,k,1)*dtcld,qci(i,k,1))
1567 qci(i,k,2) = qci(i,k,2) + pfrzdtc
1568 t(i,k,j) = t(i,k,j) + xlf/cpm(i,k)*pfrzdtc
1569 qci(i,k,1) = qci(i,k,1)-pfrzdtc
1570 endif
1571 !---------------------------------------------------------------
1572 ! psfrz: freezing of rain water [HL A20] [LFO 45]
1573 ! (T<T0, R->S)
1574 !---------------------------------------------------------------
1575 if(supcol.gt.0..and.qrs(i,k,1).gt.0.) then
1576 supcolt=min(supcol,50.)
1577 ! pfrzdtr = min(20.*pi**2*pfrz1*n0r*denr/den(i,k,j) &
1578 ! *(exp(pfrz2*supcol)-1.)*rslope(i,k,1)**7*dtcld, &
1579 ! qrs(i,k,1))
1580 temp = rslope(i,k,1)
1581 temp = temp*temp*temp*temp*temp*temp*temp
1582 pfrzdtr = min(20.*(pi*pi)*pfrz1*n0r*denr/den(i,k,j) &
1583 *(exp(pfrz2*supcolt)-1.)*temp*dtcld, &
1584 qrs(i,k,1))
1585 qrs(i,k,2) = qrs(i,k,2) + pfrzdtr
1586 t(i,k,j) = t(i,k,j) + xlf/cpm(i,k)*pfrzdtr
1587 qrs(i,k,1) = qrs(i,k,1)-pfrzdtr
1588 endif
1589 enddo
1590 !
1591 !----------------------------------------------------------------
1592 ! rsloper: reverse of the slope parameter of the rain(m)
1593 ! xka: thermal conductivity of air(jm-1s-1k-1)
1594 ! work1: the thermodynamic term in the denominator associated with
1595 ! heat conduction and vapor diffusion
1596 ! (ry88, y93, h85)
1597 ! work2: parameter associated with the ventilation effects(y93)
1598 !
1599 do k = kts, kte
1600 if(qrs(i,k,1).le.qcrmin)then
1601 rslope(i,k,1) = rslopermax
1602 rslopeb(i,k,1) = rsloperbmax
1603 rslope2(i,k,1) = rsloper2max
1604 rslope3(i,k,1) = rsloper3max
1605 else
1606 ! rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k,j))
1607 rslope(i,k,1) = 1./(sqrt(sqrt(pidn0r/((qrs(i,k,1))*(den(i,k,j))))))
1608 rslopeb(i,k,1) = exp(log(rslope(i,k,1))*(bvtr))
1609 rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
1610 rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
1611 endif
1612 if(qrs(i,k,2).le.qcrmin)then
1613 rslope(i,k,2) = rslopesmax
1614 rslopeb(i,k,2) = rslopesbmax
1615 rslope2(i,k,2) = rslopes2max
1616 rslope3(i,k,2) = rslopes3max
1617 else
1618 ! rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k,j),n0sfac(i,k))
1619 rslope(i,k,2) = 1./(sqrt(sqrt(pidn0s*(n0sfac(i,k))/((qrs(i,k,2)) &
1620 *(den(i,k,j))))))
1621 rslopeb(i,k,2) = exp(log(rslope(i,k,2))*(bvts))
1622 rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
1623 rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
1624 endif
1625 enddo
1626 !
1627 do k = kts, kte
1628 ! work1(i,k,1) = diffac(xl(i,k),p(i,k,j),t(i,k,j),den(i,k,j),qs(i,k,1))
1629 work1(i,k,1) = ((((den(i,k,j))*(xl(i,k))*(xl(i,k)))*((t(i,k,j))+120.) &
1630 *(den(i,k,j)))/(1.414e3*(1.496e-6*((t(i,k,j))*sqrt(t(i,k,j))))&
1631 *(den(i,k,j))*(rv*(t(i,k,j))*(t(i,k,j))))) &
1632 + p(i,k,j)/((qs(i,k,1))*(8.794e-5*exp(log(t(i,k,j))*(1.81))))
1633 ! work1(i,k,2) = diffac(xls,p(i,k,j),t(i,k,j),den(i,k,j),qs(i,k,2))
1634 work1(i,k,2) = ((((den(i,k,j))*(xls)*(xls))*((t(i,k,j))+120.)*(den(i,k,j)))&
1635 /(1.414e3*(1.496e-6*((t(i,k,j))*sqrt(t(i,k,j))))*(den(i,k,j)) &
1636 *(rv*(t(i,k,j))*(t(i,k,j)))) &
1637 + p(i,k,j)/(qs(i,k,2)*(8.794e-5*exp(log(t(i,k,j))*(1.81)))))
1638 ! work2(i,k) = venfac(p(i,k,j),t(i,k,j),den(i,k,j))
1639 work2(i,k) = (exp(.3333333*log(((1.496e-6 * ((t(i,k,j))*sqrt(t(i,k,j)))) &
1640 *p(i,k,j))/(((t(i,k,j))+120.)*den(i,k,j)*(8.794e-5 &
1641 *exp(log(t(i,k,j))*(1.81))))))*sqrt(sqrt(den0/(den(i,k,j))))) &
1642 /sqrt((1.496e-6*((t(i,k,j))*sqrt(t(i,k,j)))) &
1643 /(((t(i,k,j))+120.)*den(i,k,j)))
1644 enddo
1645 !
1646 !===============================================================
1647 !
1648 ! warm rain processes
1649 !
1650 ! - follows the processes in RH83 and LFO except for autoconcersion
1651 !
1652 !===============================================================
1653 !
1654 do k = kts, kte
1655 supsat = max(q(i,k,j),qmin)-qs(i,k,1)
1656 satdt = supsat/dtcld
1657 !---------------------------------------------------------------
1658 ! praut: auto conversion rate from cloud to rain [HDC 16]
1659 ! (C->R)
1660 !---------------------------------------------------------------
1661 if(qci(i,k,1).gt.qc0) then
1662 praut(i,k) = qck1*exp(log(qci(i,k,1))*((7./3.)))
1663 praut(i,k) = min(praut(i,k),qci(i,k,1)/dtcld)
1664 endif
1665 !---------------------------------------------------------------
1666 ! pracw: accretion of cloud water by rain [HL A40] [LFO 51]
1667 ! (C->R)
1668 !---------------------------------------------------------------
1669 if(qrs(i,k,1).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
1670 pracw(i,k) = min(pacrr*rslope3(i,k,1)*rslopeb(i,k,1) &
1671 *qci(i,k,1)*denfac(i,k),qci(i,k,1)/dtcld)
1672 endif
1673 !---------------------------------------------------------------
1674 ! prevp: evaporation/condensation rate of rain [HDC 14]
1675 ! (V->R or R->V)
1676 !---------------------------------------------------------------
1677 if(qrs(i,k,1).gt.0.) then
1678 coeres = rslope2(i,k,1)*sqrt(rslope(i,k,1)*rslopeb(i,k,1))
1679 prevp(i,k) = (rh(i,k,1)-1.)*(precr1*rslope2(i,k,1) &
1680 +precr2*work2(i,k)*coeres)/work1(i,k,1)
1681 if(prevp(i,k).lt.0.) then
1682 prevp(i,k) = max(prevp(i,k),-qrs(i,k,1)/dtcld)
1683 prevp(i,k) = max(prevp(i,k),satdt/2)
1684 else
1685 prevp(i,k) = min(prevp(i,k),satdt/2)
1686 endif
1687 endif
1688 enddo
1689 !
1690 !===============================================================
1691 !
1692 ! cold rain processes
1693 !
1694 ! - follows the revised ice microphysics processes in HDC
1695 ! - the processes same as in RH83 and RH84 and LFO behave
1696 ! following ice crystal hapits defined in HDC, inclduing
1697 ! intercept parameter for snow (n0s), ice crystal number
1698 ! concentration (ni), ice nuclei number concentration
1699 ! (n0i), ice diameter (d)
1700 !
1701 !===============================================================
1702 !
1703 rdtcld = 1./dtcld
1704 do k = kts, kte
1705 supcol = t0c-t(i,k,j)
1706 supsat = max(q(i,k,j),qmin)-qs(i,k,2)
1707 satdt = supsat/dtcld
1708 ifsat = 0
1709 !-------------------------------------------------------------
1710 ! Ni: ice crystal number concentraiton [HDC 5c]
1711 !-------------------------------------------------------------
1712 ! xni(i,k) = min(max(5.38e7*(den(i,k,j) &
1713 ! *max(qci(i,k,2),qmin))**0.75,1.e3),1.e6)
1714 temp = (den(i,k,j)*max(qci(i,k,2),qmin))
1715 temp = sqrt(sqrt(temp*temp*temp))
1716 xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
1717 eacrs = exp(0.07*(-supcol))
1718 !
1719 if(supcol.gt.0) then
1720 if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,2).gt.qmin) then
1721 xmi = den(i,k,j)*qci(i,k,2)/xni(i,k)
1722 diameter = min(dicon * sqrt(xmi),dimax)
1723 vt2i = 1.49e4*diameter**1.31
1724 vt2s = pvts*rslopeb(i,k,2)*denfac(i,k)
1725 !-------------------------------------------------------------
1726 ! psaci: Accretion of cloud ice by rain [HDC 10]
1727 ! (T<T0: I->S)
1728 !-------------------------------------------------------------
1729 acrfac = 2.*rslope3(i,k,2)+2.*diameter*rslope2(i,k,2) &
1730 +diameter**2*rslope(i,k,2)
1731 psaci(i,k) = pi*qci(i,k,2)*eacrs*n0s*n0sfac(i,k) &
1732 *abs(vt2s-vt2i)*acrfac/4.
1733 endif
1734 endif
1735 !-------------------------------------------------------------
1736 ! psacw: Accretion of cloud water by snow [HL A7] [LFO 24]
1737 ! (T<T0: C->S, and T>=T0: C->R)
1738 !-------------------------------------------------------------
1739 if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
1740 psacw(i,k) = min(pacrc*n0sfac(i,k)*rslope3(i,k,2) &
1741 *rslopeb(i,k,2)*qci(i,k,1)*denfac(i,k) &
1742 ! ,qci(i,k,1)/dtcld)
1743 ,qci(i,k,1)*rdtcld)
1744 endif
1745 if(supcol .gt. 0) then
1746 !-------------------------------------------------------------
1747 ! pidep: Deposition/Sublimation rate of ice [HDC 9]
1748 ! (T<T0: V->I or I->V)
1749 !-------------------------------------------------------------
1750 if(qci(i,k,2).gt.0.and.ifsat.ne.1) then
1751 xmi = den(i,k,j)*qci(i,k,2)/xni(i,k)
1752 diameter = dicon * sqrt(xmi)
1753 pidep(i,k) = 4.*diameter*xni(i,k)*(rh(i,k,2)-1.)/work1(i,k,2)
1754 supice = satdt-prevp(i,k)
1755 if(pidep(i,k).lt.0.) then
1756 ! pidep(i,k) = max(max(pidep(i,k),satdt/2),supice)
1757 ! pidep(i,k) = max(pidep(i,k),-qci(i,k,2)/dtcld)
1758 pidep(i,k) = max(max(pidep(i,k),satdt*.5),supice)
1759 pidep(i,k) = max(pidep(i,k),-qci(i,k,2)*rdtcld)
1760 else
1761 ! pidep(i,k) = min(min(pidep(i,k),satdt/2),supice)
1762 pidep(i,k) = min(min(pidep(i,k),satdt*.5),supice)
1763 endif
1764 if(abs(prevp(i,k)+pidep(i,k)).ge.abs(satdt)) ifsat = 1
1765 endif
1766 !-------------------------------------------------------------
1767 ! psdep: deposition/sublimation rate of snow [HDC 14]
1768 ! (V->S or S->V)
1769 !-------------------------------------------------------------
1770 if(qrs(i,k,2).gt.0..and.ifsat.ne.1) then
1771 coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
1772 psdep(i,k) = (rh(i,k,2)-1.)*n0sfac(i,k) &
1773 *(precs1*rslope2(i,k,2)+precs2 &
1774 *work2(i,k)*coeres)/work1(i,k,2)
1775 supice = satdt-prevp(i,k)-pidep(i,k)
1776 if(psdep(i,k).lt.0.) then
1777 ! psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)/dtcld)
1778 ! psdep(i,k) = max(max(psdep(i,k),satdt/2),supice)
1779 psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)*rdtcld)
1780 psdep(i,k) = max(max(psdep(i,k),satdt*.5),supice)
1781 else
1782 ! psdep(i,k) = min(min(psdep(i,k),satdt/2),supice)
1783 psdep(i,k) = min(min(psdep(i,k),satdt*.5),supice)
1784 endif
1785 if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)).ge.abs(satdt)) &
1786 ifsat = 1
1787 endif
1788 !-------------------------------------------------------------
1789 ! pigen: generation(nucleation) of ice from vapor [HL A50] [HDC 7-8]
1790 ! (T<T0: V->I)
1791 !-------------------------------------------------------------
1792 if(supsat.gt.0.and.ifsat.ne.1) then
1793 supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)
1794 xni0 = 1.e3*exp(0.1*supcol)
1795 roqi0 = 4.92e-11*exp(log(xni0)*(1.33))
1796 pigen(i,k) = max(0.,(roqi0/den(i,k,j)-max(qci(i,k,2),0.)) &
1797 ! /dtcld)
1798 *rdtcld)
1799 pigen(i,k) = min(min(pigen(i,k),satdt),supice)
1800 endif
1801 !
1802 !-------------------------------------------------------------
1803 ! psaut: conversion(aggregation) of ice to snow [HDC 12]
1804 ! (T<T0: I->S)
1805 !-------------------------------------------------------------
1806 if(qci(i,k,2).gt.0.) then
1807 qimax = roqimax/den(i,k,j)
1808 ! psaut(i,k) = max(0.,(qci(i,k,2)-qimax)/dtcld)
1809 psaut(i,k) = max(0.,(qci(i,k,2)-qimax)*rdtcld)
1810 endif
1811 endif
1812 !-------------------------------------------------------------
1813 ! psevp: Evaporation of melting snow [HL A35] [RH83 A27]
1814 ! (T>T0: S->V)
1815 !-------------------------------------------------------------
1816 if(supcol.lt.0.) then
1817 if(qrs(i,k,2).gt.0..and.rh(i,k,1).lt.1.) &
1818 psevp(i,k) = psdep(i,k)*work1(i,k,2)/work1(i,k,1)
1819 ! psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)/dtcld),0.)
1820 psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)*rdtcld),0.)
1821 endif
1822 enddo
1823 !
1824 !
1825 !----------------------------------------------------------------
1826 ! check mass conservation of generation terms and feedback to the
1827 ! large scale
1828 !
1829 do k = kts, kte
1830 if(t(i,k,j).le.t0c) then
1831 !
1832 ! cloud water
1833 !
1834 value = max(qmin,qci(i,k,1))
1835 source = (praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
1836 if (source.gt.value) then
1837 factor = value/source
1838 praut(i,k) = praut(i,k)*factor
1839 pracw(i,k) = pracw(i,k)*factor
1840 psacw(i,k) = psacw(i,k)*factor
1841 endif
1842 !
1843 ! cloud ice
1844 !
1845 value = max(qmin,qci(i,k,2))
1846 source = (psaut(i,k)+psaci(i,k)-pigen(i,k)-pidep(i,k))*dtcld
1847 if (source.gt.value) then
1848 factor = value/source
1849 psaut(i,k) = psaut(i,k)*factor
1850 psaci(i,k) = psaci(i,k)*factor
1851 pigen(i,k) = pigen(i,k)*factor
1852 pidep(i,k) = pidep(i,k)*factor
1853 endif
1854 !
1855 ! rain
1856 !
1857 !
1858 value = max(qmin,qrs(i,k,1))
1859 source = (-praut(i,k)-pracw(i,k)-prevp(i,k))*dtcld
1860 if (source.gt.value) then
1861 factor = value/source
1862 praut(i,k) = praut(i,k)*factor
1863 pracw(i,k) = pracw(i,k)*factor
1864 prevp(i,k) = prevp(i,k)*factor
1865 endif
1866 !
1867 ! snow
1868 !
1869 value = max(qmin,qrs(i,k,2))
1870 source = (-psdep(i,k)-psaut(i,k)-psaci(i,k)-psacw(i,k))*dtcld
1871 if (source.gt.value) then
1872 factor = value/source
1873 psdep(i,k) = psdep(i,k)*factor
1874 psaut(i,k) = psaut(i,k)*factor
1875 psaci(i,k) = psaci(i,k)*factor
1876 psacw(i,k) = psacw(i,k)*factor
1877 endif
1878 !
1879 work2(i,k)=-(prevp(i,k)+psdep(i,k)+pigen(i,k)+pidep(i,k))
1880 ! update
1881 q(i,k,j) = q(i,k,j)+work2(i,k)*dtcld
1882 qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
1883 +psacw(i,k))*dtcld,0.)
1884 qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
1885 +prevp(i,k))*dtcld,0.)
1886 qci(i,k,2) = max(qci(i,k,2)-(psaut(i,k)+psaci(i,k) &
1887 -pigen(i,k)-pidep(i,k))*dtcld,0.)
1888 qrs(i,k,2) = max(qrs(i,k,2)+(psdep(i,k)+psaut(i,k) &
1889 +psaci(i,k)+psacw(i,k))*dtcld,0.)
1890 xlf = xls-xl(i,k)
1891 xlwork2 = -xls*(psdep(i,k)+pidep(i,k)+pigen(i,k)) &
1892 -xl(i,k)*prevp(i,k)-xlf*psacw(i,k)
1893 t(i,k,j) = t(i,k,j)-xlwork2/cpm(i,k)*dtcld
1894 else
1895 !
1896 ! cloud water
1897 !
1898 value = max(qmin,qci(i,k,1))
1899 source=(praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
1900 if (source.gt.value) then
1901 factor = value/source
1902 praut(i,k) = praut(i,k)*factor
1903 pracw(i,k) = pracw(i,k)*factor
1904 psacw(i,k) = psacw(i,k)*factor
1905 endif
1906 !
1907 ! rain
1908 !
1909 value = max(qmin,qrs(i,k,1))
1910 source = (-praut(i,k)-pracw(i,k)-prevp(i,k)-psacw(i,k))*dtcld
1911 if (source.gt.value) then
1912 factor = value/source
1913 praut(i,k) = praut(i,k)*factor
1914 pracw(i,k) = pracw(i,k)*factor
1915 prevp(i,k) = prevp(i,k)*factor
1916 psacw(i,k) = psacw(i,k)*factor
1917 endif
1918 !
1919 ! snow
1920 !
1921 value = max(qcrmin,qrs(i,k,2))
1922 source=(-psevp(i,k))*dtcld
1923 if (source.gt.value) then
1924 factor = value/source
1925 psevp(i,k) = psevp(i,k)*factor
1926 endif
1927 work2(i,k)=-(prevp(i,k)+psevp(i,k))
1928 ! update
1929 q(i,k,j) = q(i,k,j)+work2(i,k)*dtcld
1930 qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
1931 +psacw(i,k))*dtcld,0.)
1932 qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
1933 +prevp(i,k) +psacw(i,k))*dtcld,0.)
1934 qrs(i,k,2) = max(qrs(i,k,2)+psevp(i,k)*dtcld,0.)
1935 xlf = xls-xl(i,k)
1936 xlwork2 = -xl(i,k)*(prevp(i,k)+psevp(i,k))
1937 t(i,k,j) = t(i,k,j)-xlwork2/cpm(i,k)*dtcld
1938 endif
1939 enddo
1940 !
1941 ! Inline expansion for fpvs
1942 ! qs(i,k,1) = fpvs(t(i,k,j),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1943 ! qs(i,k,2) = fpvs(t(i,k,j),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
1944 hsub = xls
1945 hvap = xlv0
1946 cvap = cpv
1947 ttp=t0c+0.01
1948 dldt=cvap-cliq
1949 xa=-dldt/rv
1950 xb=xa+hvap/(rv*ttp)
1951 dldti=cvap-cice
1952 xai=-dldti/rv
1953 xbi=xai+hsub/(rv*ttp)
1954 do k = kts, kte
1955 tr=ttp/t(i,k,j)
1956 logtr = log(tr)
1957 qs(i,k,1)=psat*exp(logtr*(xa)+xb*(1.-tr))
1958 qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k,j) - qs(i,k,1))
1959 qs(i,k,1) = max(qs(i,k,1),qmin)
1960 enddo
1961 !
1962 !----------------------------------------------------------------
1963 ! pcond: condensational/evaporational rate of cloud water [HL A46] [RH83 A6]
1964 ! if there exists additional water vapor condensated/if
1965 ! evaporation of cloud water is not enough to remove subsaturation
1966 !
1967 do k = kts, kte
1968 ! work1(i,k,1) = conden(t(i,k,j),q(i,k,j),qs(i,k,1),xl(i,k),cpm(i,k))
1969 work1(i,k,1) = ((max(q(i,k,j),qmin)-(qs(i,k,1)))/(1.+(xl(i,k)) &
1970 *(xl(i,k))/(rv*(cpm(i,k)))*(qs(i,k,1)) &
1971 /((t(i,k,j))*(t(i,k,j)))))
1972 work2(i,k) = qci(i,k,1)+work1(i,k,1)
1973 pcond(i,k) = min(max(work1(i,k,1)/dtcld,0.),max(q(i,k,j),0.)/dtcld)
1974 if(qci(i,k,1).gt.0..and.work1(i,k,1).lt.0.) &
1975 pcond(i,k) = max(work1(i,k,1),-qci(i,k,1))/dtcld
1976 q(i,k,j) = q(i,k,j)-pcond(i,k)*dtcld
1977 qci(i,k,1) = max(qci(i,k,1)+pcond(i,k)*dtcld,0.)
1978 t(i,k,j) = t(i,k,j)+pcond(i,k)*xl(i,k)/cpm(i,k)*dtcld
1979 enddo
1980 !
1981 !
1982 !----------------------------------------------------------------
1983 ! padding for small values
1984 !
1985 do k = kts, kte
1986 if(qci(i,k,1).le.qmin) qci(i,k,1) = 0.0
1987 if(qci(i,k,2).le.qmin) qci(i,k,2) = 0.0
1988 enddo
1989 enddo ! big loops
1991 DO K=kts,kte
1992 th(i,k,j)=t(i,k,j)/pii(i,k,j)
1993 qc(i,k,j) = qci(i,k,1)
1994 qi(i,k,j) = qci(i,k,2)
1995 qr(i,k,j) = qrs(i,k,1)
1996 qqs(i,k,j) = qrs(i,k,2)
1997 ENDDO
1999 ENDDO ! i loop
2000 enddo ! j loop
2001 !$acc end region
2003 ENDIF
2005 END SUBROUTINE wsm52d
2008 #else
2010 !===================================================================
2011 !
2012 SUBROUTINE wsm52D(t, q, qci, qrs, den, p, delz &
2013 ,delt,g, cpd, cpv, rd, rv, t0c &
2014 ,ep1, ep2, qmin &
2015 ,XLS, XLV0, XLF0, den0, denr &
2016 ,cliq,cice,psat &
2017 ,lat &
2018 ,rain,rainncv &
2019 ,sr &
2020 ,ids,ide, jds,jde, kds,kde &
2021 ,ims,ime, jms,jme, kms,kme &
2022 ,its,ite, jts,jte, kts,kte &
2023 ,snow,snowncv &
2024 )
2025 !-------------------------------------------------------------------
2026 IMPLICIT NONE
2027 !-------------------------------------------------------------------
2028 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
2029 ims,ime, jms,jme, kms,kme , &
2030 its,ite, jts,jte, kts,kte, &
2031 lat
2032 REAL, DIMENSION( its:ite , kts:kte ), &
2033 INTENT(INOUT) :: &
2034 t
2035 REAL, DIMENSION( its:ite , kts:kte, 2 ), &
2036 INTENT(INOUT) :: &
2037 qci, &
2038 qrs
2039 REAL, DIMENSION( ims:ime , kms:kme ), &
2040 INTENT(INOUT) :: &
2041 q
2042 REAL, DIMENSION( ims:ime , kms:kme ), &
2043 INTENT(IN ) :: &
2044 den, &
2045 p, &
2046 delz
2047 REAL, INTENT(IN ) :: delt, &
2048 g, &
2049 cpd, &
2050 cpv, &
2051 t0c, &
2052 den0, &
2053 rd, &
2054 rv, &
2055 ep1, &
2056 ep2, &
2057 qmin, &
2058 XLS, &
2059 XLV0, &
2060 XLF0, &
2061 cliq, &
2062 cice, &
2063 psat, &
2064 denr
2065 REAL, DIMENSION( ims:ime ), &
2066 INTENT(INOUT) :: rain, &
2067 rainncv, &
2068 sr
2070 REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, &
2071 INTENT(INOUT) :: snow, &
2072 snowncv
2074 ! LOCAL VAR
2075 REAL, DIMENSION( its:ite , kts:kte , 2) :: &
2076 rh, &
2077 qs, &
2078 rslope, &
2079 rslope2, &
2080 rslope3, &
2081 rslopeb, &
2082 falk, &
2083 fall, &
2084 work1
2086 REAL, DIMENSION( its:ite , kts:kte ) :: &
2087 falkc, &
2088 fallc, &
2089 xl, &
2090 cpm, &
2091 denfac, &
2092 xni, &
2093 n0sfac, &
2094 work2, &
2095 work1c, &
2096 work2c
2098 REAL, DIMENSION( its:ite , kts:kte ) :: &
2099 pigen, &
2100 pidep, &
2101 psdep, &
2102 praut, &
2103 psaut, &
2104 prevp, &
2105 psevp, &
2106 pracw, &
2107 psacw, &
2108 psaci, &
2109 pcond, &
2110 psmlt
2111 INTEGER, DIMENSION( its:ite ) :: &
2112 mstep, &
2113 numdt
2114 REAL, DIMENSION(its:ite) :: rmstep
2115 REAL dtcldden, rdelz, rdtcld
2116 LOGICAL, DIMENSION( its:ite ) :: flgcld
2118 #define WSM_NO_CONDITIONAL_IN_VECTOR
2119 #ifdef WSM_NO_CONDITIONAL_IN_VECTOR
2120 REAL, DIMENSION(its:ite) :: xal, xbl
2121 #endif
2123 REAL :: pi, &
2124 cpmcal, xlcal, lamdar, lamdas, diffus, &
2125 viscos, xka, venfac, conden, diffac, &
2126 x, y, z, a, b, c, d, e, &
2127 qdt, holdrr, holdrs, supcol, supcolt, pvt, &
2128 coeres, supsat, dtcld, xmi, eacrs, satdt, &
2129 vt2i,vt2s,acrfac, &
2130 qimax, diameter, xni0, roqi0, &
2131 fallsum, fallsum_qsi, xlwork2, factor, source, &
2132 value, xlf, pfrzdtc, pfrzdtr, supice, holdc, holdci
2133 ! variables for optimization
2134 REAL, DIMENSION( its:ite ) :: tvec1
2135 REAL :: temp
2136 INTEGER :: i, j, k, mstepmax, &
2137 iprt, latd, lond, loop, loops, ifsat, n
2138 ! Temporaries used for inlining fpvs function
2139 REAL :: dldti, xb, xai, tr, xbi, xa, hvap, cvap, hsub, dldt, ttp
2140 REAL :: logtr
2141 !
2142 !=================================================================
2143 ! compute internal functions
2144 !
2145 cpmcal(x) = cpd*(1.-max(x,qmin))+max(x,qmin)*cpv
2146 xlcal(x) = xlv0-xlv1*(x-t0c)
2147 !----------------------------------------------------------------
2148 ! size distributions: (x=mixing ratio, y=air density):
2149 ! valid for mixing ratio > 1.e-9 kg/kg.
2150 !
2151 ! Optimizatin : A**B => exp(log(A)*(B))
2152 lamdar(x,y)= sqrt(sqrt(pidn0r/(x*y))) ! (pidn0r/(x*y))**.25
2153 lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y))) ! (pidn0s*z/(x*y))**.25
2154 !
2155 !----------------------------------------------------------------
2156 ! diffus: diffusion coefficient of the water vapor
2157 ! viscos: kinematic viscosity(m2s-1)
2158 ! diffus(x,y) = 8.794e-5 * exp(log(x)*(1.81)) / y ! 8.794e-5*x**1.81/y
2159 ! viscos(x,y) = 1.496e-6 * (x*sqrt(x)) /(x+120.)/y ! 1.496e-6*x**1.5/(x+120.)/y
2160 ! xka(x,y) = 1.414e3*viscos(x,y)*y
2161 ! diffac(a,b,c,d,e) = d*a*a/(xka(c,d)*rv*c*c)+1./(e*diffus(c,b))
2162 ! venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333))) &
2163 ! /sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
2164 ! conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+d*d/(rv*e)*c/(a*a))
2165 !
2166 !
2167 pi = 4. * atan(1.)
2168 !
2169 !----------------------------------------------------------------
2170 ! paddint 0 for negative values generated by dynamics
2171 !
2172 do k = kts, kte
2173 do i = its, ite
2174 qci(i,k,1) = max(qci(i,k,1),0.0)
2175 qrs(i,k,1) = max(qrs(i,k,1),0.0)
2176 qci(i,k,2) = max(qci(i,k,2),0.0)
2177 qrs(i,k,2) = max(qrs(i,k,2),0.0)
2178 enddo
2179 enddo
2180 !
2181 !----------------------------------------------------------------
2182 ! latent heat for phase changes and heat capacity. neglect the
2183 ! changes during microphysical process calculation
2184 ! emanuel(1994)
2185 !
2186 do k = kts, kte
2187 do i = its, ite
2188 cpm(i,k) = cpmcal(q(i,k))
2189 xl(i,k) = xlcal(t(i,k))
2190 enddo
2191 enddo
2192 !
2193 !----------------------------------------------------------------
2194 ! compute the minor time steps.
2195 !
2196 loops = max(nint(delt/dtcldcr),1)
2197 dtcld = delt/loops
2198 if(delt.le.dtcldcr) dtcld = delt
2199 !
2200 do loop = 1,loops
2201 !
2202 !----------------------------------------------------------------
2203 ! initialize the large scale variables
2204 !
2205 do i = its, ite
2206 mstep(i) = 1
2207 flgcld(i) = .true.
2208 enddo
2209 !
2210 ! do k = kts, kte
2211 ! do i = its, ite
2212 ! denfac(i,k) = sqrt(den0/den(i,k))
2213 ! enddo
2214 ! enddo
2215 do k = kts, kte
2216 CALL VREC( tvec1(its), den(its,k), ite-its+1)
2217 do i = its, ite
2218 tvec1(i) = tvec1(i)*den0
2219 enddo
2220 CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
2221 enddo
2222 !
2223 ! Inline expansion for fpvs
2224 ! qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
2225 ! qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
2226 hsub = xls
2227 hvap = xlv0
2228 cvap = cpv
2229 ttp=t0c+0.01
2230 dldt=cvap-cliq
2231 xa=-dldt/rv
2232 xb=xa+hvap/(rv*ttp)
2233 dldti=cvap-cice
2234 xai=-dldti/rv
2235 xbi=xai+hsub/(rv*ttp)
2237 ! this is for compilers where the conditional inhibits vectorization
2238 #ifdef WSM_NO_CONDITIONAL_IN_VECTOR
2239 do k = kts, kte
2240 do i = its, ite
2241 if(t(i,k).lt.ttp) then
2242 xal(i) = xai
2243 xbl(i) = xbi
2244 else
2245 xal(i) = xa
2246 xbl(i) = xb
2247 endif
2248 enddo
2249 do i = its, ite
2250 tr=ttp/t(i,k)
2251 logtr=log(tr)
2252 qs(i,k,1)=psat*exp(logtr*(xa)+xb*(1.-tr))
2253 qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
2254 qs(i,k,1) = max(qs(i,k,1),qmin)
2255 rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
2256 qs(i,k,2)=psat*exp(logtr*(xal(i))+xbl(i)*(1.-tr))
2257 qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
2258 qs(i,k,2) = max(qs(i,k,2),qmin)
2259 rh(i,k,2) = max(q(i,k) / qs(i,k,2),qmin)
2260 enddo
2261 enddo
2262 #else
2263 do k = kts, kte
2264 do i = its, ite
2265 tr=ttp/t(i,k)
2266 logtr=log(tr)
2267 qs(i,k,1)=psat*exp(logtr*(xa)+xb*(1.-tr))
2268 qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
2269 qs(i,k,1) = max(qs(i,k,1),qmin)
2270 rh(i,k,1) = max(q(i,k) / qs(i,k,1),qmin)
2271 if(t(i,k).lt.ttp) then
2272 qs(i,k,2)=psat*exp(logtr*(xai)+xbi*(1.-tr))
2273 else
2274 qs(i,k,2)=psat*exp(logtr*(xa)+xb*(1.-tr))
2275 endif
2276 qs(i,k,2) = ep2 * qs(i,k,2) / (p(i,k) - qs(i,k,2))
2277 qs(i,k,2) = max(qs(i,k,2),qmin)
2278 rh(i,k,2) = max(q(i,k) / qs(i,k,2),qmin)
2279 enddo
2280 enddo
2281 #endif
2282 !
2283 !----------------------------------------------------------------
2284 ! initialize the variables for microphysical physics
2285 !
2286 !
2287 do k = kts, kte
2288 do i = its, ite
2289 prevp(i,k) = 0.
2290 psdep(i,k) = 0.
2291 praut(i,k) = 0.
2292 psaut(i,k) = 0.
2293 pracw(i,k) = 0.
2294 psaci(i,k) = 0.
2295 psacw(i,k) = 0.
2296 pigen(i,k) = 0.
2297 pidep(i,k) = 0.
2298 pcond(i,k) = 0.
2299 psmlt(i,k) = 0.
2300 psevp(i,k) = 0.
2301 falk(i,k,1) = 0.
2302 falk(i,k,2) = 0.
2303 fall(i,k,1) = 0.
2304 fall(i,k,2) = 0.
2305 fallc(i,k) = 0.
2306 falkc(i,k) = 0.
2307 xni(i,k) = 1.e3
2308 enddo
2309 enddo
2310 !
2311 !----------------------------------------------------------------
2312 ! compute the fallout term:
2313 ! first, vertical terminal velosity for minor loops
2314 !
2315 do k = kts, kte
2316 do i = its, ite
2317 supcol = t0c-t(i,k)
2318 !---------------------------------------------------------------
2319 ! n0s: Intercept parameter for snow [m-4] [HDC 6]
2320 !---------------------------------------------------------------
2321 n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
2322 if(qrs(i,k,1).le.qcrmin)then
2323 rslope(i,k,1) = rslopermax
2324 rslopeb(i,k,1) = rsloperbmax
2325 rslope2(i,k,1) = rsloper2max
2326 rslope3(i,k,1) = rsloper3max
2327 else
2328 rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k))
2329 rslopeb(i,k,1) = exp(log(rslope(i,k,1))*(bvtr))
2330 rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
2331 rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
2332 endif
2333 if(qrs(i,k,2).le.qcrmin)then
2334 rslope(i,k,2) = rslopesmax
2335 rslopeb(i,k,2) = rslopesbmax
2336 rslope2(i,k,2) = rslopes2max
2337 rslope3(i,k,2) = rslopes3max
2338 else
2339 rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k))
2340 rslopeb(i,k,2) = exp(log(rslope(i,k,2))*(bvts))
2341 rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
2342 rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
2343 endif
2344 !-------------------------------------------------------------
2345 ! Ni: ice crystal number concentraiton [HDC 5c]
2346 !-------------------------------------------------------------
2347 ! xni(i,k) = min(max(5.38e7*(den(i,k) &
2348 ! *max(qci(i,k,2),qmin))**0.75,1.e3),1.e6)
2349 temp = (den(i,k)*max(qci(i,k,2),qmin))
2350 temp = sqrt(sqrt(temp*temp*temp))
2351 xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
2352 enddo
2353 enddo
2354 !
2355 mstepmax = 1
2356 numdt = 1
2357 do k = kte, kts, -1
2358 do i = its, ite
2359 work1(i,k,1) = pvtr*rslopeb(i,k,1)*denfac(i,k)/delz(i,k)
2360 work1(i,k,2) = pvts*rslopeb(i,k,2)*denfac(i,k)/delz(i,k)
2361 numdt(i) = max(nint(max(work1(i,k,1),work1(i,k,2))*dtcld+.5),1)
2362 if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
2363 enddo
2364 enddo
2365 do i = its, ite
2366 if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
2367 rmstep(i) = 1./mstep(i)
2368 enddo
2369 !
2370 do n = 1, mstepmax
2371 k = kte
2372 do i = its, ite
2373 if(n.le.mstep(i)) then
2374 ! falk(i,k,1) = den(i,k)*qrs(i,k,1)*work1(i,k,1)/mstep(i)
2375 ! falk(i,k,2) = den(i,k)*qrs(i,k,2)*work1(i,k,2)/mstep(i)
2376 falk(i,k,1) = den(i,k)*qrs(i,k,1)*work1(i,k,1)*rmstep(i)
2377 falk(i,k,2) = den(i,k)*qrs(i,k,2)*work1(i,k,2)*rmstep(i)
2378 fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
2379 fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
2380 ! qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcld/den(i,k),0.)
2381 ! qrs(i,k,2) = max(qrs(i,k,2)-falk(i,k,2)*dtcld/den(i,k),0.)
2382 dtcldden = dtcld/den(i,k)
2383 qrs(i,k,1) = max(qrs(i,k,1)-falk(i,k,1)*dtcldden,0.)
2384 qrs(i,k,2) = max(qrs(i,k,2)-falk(i,k,2)*dtcldden,0.)
2385 endif
2386 enddo
2387 do k = kte-1, kts, -1
2388 do i = its, ite
2389 if(n.le.mstep(i)) then
2390 falk(i,k,1) = den(i,k)*qrs(i,k,1)*work1(i,k,1)*rmstep(i)
2391 falk(i,k,2) = den(i,k)*qrs(i,k,2)*work1(i,k,2)*rmstep(i)
2392 fall(i,k,1) = fall(i,k,1)+falk(i,k,1)
2393 fall(i,k,2) = fall(i,k,2)+falk(i,k,2)
2394 dtcldden = dtcld/den(i,k)
2395 rdelz = 1./delz(i,k)
2396 qrs(i,k,1) = max(qrs(i,k,1)-(falk(i,k,1)-falk(i,k+1,1) &
2397 *delz(i,k+1)*rdelz)*dtcldden,0.)
2398 qrs(i,k,2) = max(qrs(i,k,2)-(falk(i,k,2)-falk(i,k+1,2) &
2399 *delz(i,k+1)*rdelz)*dtcldden,0.)
2400 endif
2401 enddo
2402 enddo
2403 do k = kte, kts, -1
2404 do i = its, ite
2405 if(n.le.mstep(i)) then
2406 if(t(i,k).gt.t0c.and.qrs(i,k,2).gt.0.) then
2407 !----------------------------------------------------------------
2408 ! psmlt: melting of snow [HL A33] [RH83 A25]
2409 ! (T>T0: S->R)
2410 !----------------------------------------------------------------
2411 xlf = xlf0
2412 ! work2(i,k)= venfac(p(i,k),t(i,k),den(i,k))
2413 work2(i,k)= (exp(log(((1.496e-6*((t(i,k))*sqrt(t(i,k))) &
2414 /((t(i,k))+120.)/(den(i,k)))/(8.794e-5 &
2415 *exp(log(t(i,k))*(1.81))/p(i,k)))) &
2416 *((.3333333)))/sqrt((1.496e-6*((t(i,k)) &
2417 *sqrt(t(i,k)))/((t(i,k))+120.)/(den(i,k)))) &
2418 *sqrt(sqrt(den0/(den(i,k)))))
2419 coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
2420 ! psmlt(i,k) = xka(t(i,k),den(i,k))/xlf*(t0c-t(i,k))*pi/2. &
2421 ! *n0sfac(i,k)*(precs1*rslope2(i,k,2)+precs2 &
2422 ! *work2(i,k)*coeres)
2423 psmlt(i,k) = (1.414e3*(1.496e-6*((t(i,k))*sqrt(t(i,k))) &
2424 /((t(i,k))+120.)/(den(i,k)) )*(den(i,k))) &
2425 /xlf*(t0c-t(i,k))*pi/2. &
2426 *n0sfac(i,k)*(precs1*rslope2(i,k,2)+precs2 &
2427 *work2(i,k)*coeres)
2428 psmlt(i,k) = min(max(psmlt(i,k)*dtcld/mstep(i), &
2429 -qrs(i,k,2)/mstep(i)),0.)
2430 qrs(i,k,2) = qrs(i,k,2) + psmlt(i,k)
2431 qrs(i,k,1) = qrs(i,k,1) - psmlt(i,k)
2432 t(i,k) = t(i,k) + xlf/cpm(i,k)*psmlt(i,k)
2433 endif
2434 endif
2435 enddo
2436 enddo
2437 enddo
2438 !---------------------------------------------------------------
2439 ! Vice [ms-1] : fallout of ice crystal [HDC 5a]
2440 !---------------------------------------------------------------
2441 mstepmax = 1
2442 mstep = 1
2443 numdt = 1
2444 do k = kte, kts, -1
2445 do i = its, ite
2446 if(qci(i,k,2).le.0.) then
2447 work2c(i,k) = 0.
2448 else
2449 xmi = den(i,k)*qci(i,k,2)/xni(i,k)
2450 ! diameter = min(dicon * sqrt(xmi),dimax)
2451 diameter = max(min(dicon * sqrt(xmi),dimax), 1.e-25)
2452 work1c(i,k) = 1.49e4*exp(log(diameter)*(1.31))
2453 work2c(i,k) = work1c(i,k)/delz(i,k)
2454 endif
2455 numdt(i) = max(nint(work2c(i,k)*dtcld+.5),1)
2456 if(numdt(i).ge.mstep(i)) mstep(i) = numdt(i)
2457 enddo
2458 enddo
2459 do i = its, ite
2460 if(mstepmax.le.mstep(i)) mstepmax = mstep(i)
2461 enddo
2462 !
2463 do n = 1, mstepmax
2464 k = kte
2465 do i = its, ite
2466 if(n.le.mstep(i)) then
2467 falkc(i,k) = den(i,k)*qci(i,k,2)*work2c(i,k)/mstep(i)
2468 holdc = falkc(i,k)
2469 fallc(i,k) = fallc(i,k)+falkc(i,k)
2470 holdci = qci(i,k,2)
2471 qci(i,k,2) = max(qci(i,k,2)-falkc(i,k)*dtcld/den(i,k),0.)
2472 endif
2473 enddo
2474 do k = kte-1, kts, -1
2475 do i = its, ite
2476 if(n.le.mstep(i)) then
2477 falkc(i,k) = den(i,k)*qci(i,k,2)*work2c(i,k)/mstep(i)
2478 holdc = falkc(i,k)
2479 fallc(i,k) = fallc(i,k)+falkc(i,k)
2480 holdci = qci(i,k,2)
2481 qci(i,k,2) = max(qci(i,k,2)-(falkc(i,k)-falkc(i,k+1) &
2482 *delz(i,k+1)/delz(i,k))*dtcld/den(i,k),0.)
2483 endif
2484 enddo
2485 enddo
2486 enddo
2487 !
2488 !
2489 !----------------------------------------------------------------
2490 ! rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
2491 !
2492 do i = its, ite
2493 fallsum = fall(i,1,1)+fall(i,1,2)+fallc(i,1)
2494 fallsum_qsi = fall(i,1,2)+fallc(i,1)
2495 rainncv(i) = 0.
2496 if(fallsum.gt.0.) then
2497 rainncv(i) = fallsum*delz(i,1)/denr*dtcld*1000.
2498 rain(i) = fallsum*delz(i,1)/denr*dtcld*1000. + rain(i)
2499 endif
2500 IF ( PRESENT (snowncv) .AND. PRESENT (snow)) THEN
2501 snowncv(i,lat) = 0.
2502 if(fallsum_qsi.gt.0.) then
2503 snowncv(i,lat) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000.
2504 snow(i,lat) = fallsum_qsi*delz(i,kts)/denr*dtcld*1000. + snow(i,lat)
2505 endif
2506 ENDIF
2507 sr(i) = 0.
2508 if(fallsum.gt.0.)sr(i)=fallsum_qsi*delz(i,kts)/denr*dtcld*1000. &
2509 /(rainncv(i)+1.e-12)
2510 enddo
2511 !
2512 !---------------------------------------------------------------
2513 ! pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
2514 ! (T>T0: I->C)
2515 !---------------------------------------------------------------
2516 do k = kts, kte
2517 do i = its, ite
2518 supcol = t0c-t(i,k)
2519 xlf = xls-xl(i,k)
2520 if(supcol.lt.0.) xlf = xlf0
2521 if(supcol.lt.0.and.qci(i,k,2).gt.0.) then
2522 qci(i,k,1) = qci(i,k,1) + qci(i,k,2)
2523 t(i,k) = t(i,k) - xlf/cpm(i,k)*qci(i,k,2)
2524 qci(i,k,2) = 0.
2525 endif
2526 !---------------------------------------------------------------
2527 ! pihmf: homogeneous freezing of cloud water below -40c [HL A45]
2528 ! (T<-40C: C->I)
2529 !---------------------------------------------------------------
2530 if(supcol.gt.40..and.qci(i,k,1).gt.0.) then
2531 qci(i,k,2) = qci(i,k,2) + qci(i,k,1)
2532 t(i,k) = t(i,k) + xlf/cpm(i,k)*qci(i,k,1)
2533 qci(i,k,1) = 0.
2534 endif
2535 !---------------------------------------------------------------
2536 ! pihtf: heterogeneous freezing of cloud water [HL A44]
2537 ! (T0>T>-40C: C->I)
2538 !---------------------------------------------------------------
2539 if(supcol.gt.0..and.qci(i,k,1).gt.0.) then
2540 supcolt=min(supcol,50.)
2541 ! pfrzdtc = min(pfrz1*(exp(pfrz2*supcol)-1.) &
2542 ! *den(i,k)/denr/xncr*qci(i,k,1)**2*dtcld,qci(i,k,1))
2543 pfrzdtc = min(pfrz1*(exp(pfrz2*supcolt)-1.) &
2544 *den(i,k)/denr/xncr*qci(i,k,1)*qci(i,k,1)*dtcld,qci(i,k,1))
2545 qci(i,k,2) = qci(i,k,2) + pfrzdtc
2546 t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtc
2547 qci(i,k,1) = qci(i,k,1)-pfrzdtc
2548 endif
2549 !---------------------------------------------------------------
2550 ! psfrz: freezing of rain water [HL A20] [LFO 45]
2551 ! (T<T0, R->S)
2552 !---------------------------------------------------------------
2553 if(supcol.gt.0..and.qrs(i,k,1).gt.0.) then
2554 supcolt=min(supcol,50.)
2555 ! pfrzdtr = min(20.*pi**2*pfrz1*n0r*denr/den(i,k) &
2556 ! *(exp(pfrz2*supcol)-1.)*rslope(i,k,1)**7*dtcld, &
2557 ! qrs(i,k,1))
2558 temp = rslope(i,k,1)
2559 temp = temp*temp*temp*temp*temp*temp*temp
2560 pfrzdtr = min(20.*(pi*pi)*pfrz1*n0r*denr/den(i,k) &
2561 *(exp(pfrz2*supcolt)-1.)*temp*dtcld, &
2562 qrs(i,k,1))
2563 qrs(i,k,2) = qrs(i,k,2) + pfrzdtr
2564 t(i,k) = t(i,k) + xlf/cpm(i,k)*pfrzdtr
2565 qrs(i,k,1) = qrs(i,k,1)-pfrzdtr
2566 endif
2567 enddo
2568 enddo
2569 !
2570 !----------------------------------------------------------------
2571 ! rsloper: reverse of the slope parameter of the rain(m)
2572 ! xka: thermal conductivity of air(jm-1s-1k-1)
2573 ! work1: the thermodynamic term in the denominator associated with
2574 ! heat conduction and vapor diffusion
2575 ! (ry88, y93, h85)
2576 ! work2: parameter associated with the ventilation effects(y93)
2577 !
2578 do k = kts, kte
2579 do i = its, ite
2580 if(qrs(i,k,1).le.qcrmin)then
2581 rslope(i,k,1) = rslopermax
2582 rslopeb(i,k,1) = rsloperbmax
2583 rslope2(i,k,1) = rsloper2max
2584 rslope3(i,k,1) = rsloper3max
2585 else
2586 ! rslope(i,k,1) = 1./lamdar(qrs(i,k,1),den(i,k))
2587 rslope(i,k,1) = 1./(sqrt(sqrt(pidn0r/((qrs(i,k,1))*(den(i,k))))))
2588 rslopeb(i,k,1) = exp(log(rslope(i,k,1))*(bvtr))
2589 rslope2(i,k,1) = rslope(i,k,1)*rslope(i,k,1)
2590 rslope3(i,k,1) = rslope2(i,k,1)*rslope(i,k,1)
2591 endif
2592 if(qrs(i,k,2).le.qcrmin)then
2593 rslope(i,k,2) = rslopesmax
2594 rslopeb(i,k,2) = rslopesbmax
2595 rslope2(i,k,2) = rslopes2max
2596 rslope3(i,k,2) = rslopes3max
2597 else
2598 ! rslope(i,k,2) = 1./lamdas(qrs(i,k,2),den(i,k),n0sfac(i,k))
2599 rslope(i,k,2) = 1./(sqrt(sqrt(pidn0s*(n0sfac(i,k))/((qrs(i,k,2)) &
2600 *(den(i,k))))))
2601 rslopeb(i,k,2) = exp(log(rslope(i,k,2))*(bvts))
2602 rslope2(i,k,2) = rslope(i,k,2)*rslope(i,k,2)
2603 rslope3(i,k,2) = rslope2(i,k,2)*rslope(i,k,2)
2604 endif
2605 enddo
2606 enddo
2607 !
2608 do k = kts, kte
2609 do i = its, ite
2610 ! work1(i,k,1) = diffac(xl(i,k),p(i,k),t(i,k),den(i,k),qs(i,k,1))
2611 work1(i,k,1) = ((((den(i,k))*(xl(i,k))*(xl(i,k)))*((t(i,k))+120.) &
2612 *(den(i,k)))/(1.414e3*(1.496e-6*((t(i,k))*sqrt(t(i,k))))&
2613 *(den(i,k))*(rv*(t(i,k))*(t(i,k))))) &
2614 + p(i,k)/((qs(i,k,1))*(8.794e-5*exp(log(t(i,k))*(1.81))))
2615 ! work1(i,k,2) = diffac(xls,p(i,k),t(i,k),den(i,k),qs(i,k,2))
2616 work1(i,k,2) = ((((den(i,k))*(xls)*(xls))*((t(i,k))+120.)*(den(i,k)))&
2617 /(1.414e3*(1.496e-6*((t(i,k))*sqrt(t(i,k))))*(den(i,k)) &
2618 *(rv*(t(i,k))*(t(i,k)))) &
2619 + p(i,k)/(qs(i,k,2)*(8.794e-5*exp(log(t(i,k))*(1.81)))))
2620 ! work2(i,k) = venfac(p(i,k),t(i,k),den(i,k))
2621 work2(i,k) = (exp(.3333333*log(((1.496e-6 * ((t(i,k))*sqrt(t(i,k)))) &
2622 *p(i,k))/(((t(i,k))+120.)*den(i,k)*(8.794e-5 &
2623 *exp(log(t(i,k))*(1.81))))))*sqrt(sqrt(den0/(den(i,k))))) &
2624 /sqrt((1.496e-6*((t(i,k))*sqrt(t(i,k)))) &
2625 /(((t(i,k))+120.)*den(i,k)))
2626 enddo
2627 enddo
2628 !
2629 !===============================================================
2630 !
2631 ! warm rain processes
2632 !
2633 ! - follows the processes in RH83 and LFO except for autoconcersion
2634 !
2635 !===============================================================
2636 !
2637 do k = kts, kte
2638 do i = its, ite
2639 supsat = max(q(i,k),qmin)-qs(i,k,1)
2640 satdt = supsat/dtcld
2641 !---------------------------------------------------------------
2642 ! praut: auto conversion rate from cloud to rain [HDC 16]
2643 ! (C->R)
2644 !---------------------------------------------------------------
2645 if(qci(i,k,1).gt.qc0) then
2646 praut(i,k) = qck1*exp(log(qci(i,k,1))*((7./3.)))
2647 praut(i,k) = min(praut(i,k),qci(i,k,1)/dtcld)
2648 endif
2649 !---------------------------------------------------------------
2650 ! pracw: accretion of cloud water by rain [HL A40] [LFO 51]
2651 ! (C->R)
2652 !---------------------------------------------------------------
2653 if(qrs(i,k,1).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
2654 pracw(i,k) = min(pacrr*rslope3(i,k,1)*rslopeb(i,k,1) &
2655 *qci(i,k,1)*denfac(i,k),qci(i,k,1)/dtcld)
2656 endif
2657 !---------------------------------------------------------------
2658 ! prevp: evaporation/condensation rate of rain [HDC 14]
2659 ! (V->R or R->V)
2660 !---------------------------------------------------------------
2661 if(qrs(i,k,1).gt.0.) then
2662 coeres = rslope2(i,k,1)*sqrt(rslope(i,k,1)*rslopeb(i,k,1))
2663 prevp(i,k) = (rh(i,k,1)-1.)*(precr1*rslope2(i,k,1) &
2664 +precr2*work2(i,k)*coeres)/work1(i,k,1)
2665 if(prevp(i,k).lt.0.) then
2666 prevp(i,k) = max(prevp(i,k),-qrs(i,k,1)/dtcld)
2667 prevp(i,k) = max(prevp(i,k),satdt/2)
2668 else
2669 prevp(i,k) = min(prevp(i,k),satdt/2)
2670 endif
2671 endif
2672 enddo
2673 enddo
2674 !
2675 !===============================================================
2676 !
2677 ! cold rain processes
2678 !
2679 ! - follows the revised ice microphysics processes in HDC
2680 ! - the processes same as in RH83 and RH84 and LFO behave
2681 ! following ice crystal hapits defined in HDC, inclduing
2682 ! intercept parameter for snow (n0s), ice crystal number
2683 ! concentration (ni), ice nuclei number concentration
2684 ! (n0i), ice diameter (d)
2685 !
2686 !===============================================================
2687 !
2688 rdtcld = 1./dtcld
2689 do k = kts, kte
2690 do i = its, ite
2691 supcol = t0c-t(i,k)
2692 supsat = max(q(i,k),qmin)-qs(i,k,2)
2693 satdt = supsat/dtcld
2694 ifsat = 0
2695 !-------------------------------------------------------------
2696 ! Ni: ice crystal number concentraiton [HDC 5c]
2697 !-------------------------------------------------------------
2698 ! xni(i,k) = min(max(5.38e7*(den(i,k) &
2699 ! *max(qci(i,k,2),qmin))**0.75,1.e3),1.e6)
2700 temp = (den(i,k)*max(qci(i,k,2),qmin))
2701 temp = sqrt(sqrt(temp*temp*temp))
2702 xni(i,k) = min(max(5.38e7*temp,1.e3),1.e6)
2703 eacrs = exp(0.07*(-supcol))
2704 !
2705 if(supcol.gt.0) then
2706 if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,2).gt.qmin) then
2707 xmi = den(i,k)*qci(i,k,2)/xni(i,k)
2708 diameter = min(dicon * sqrt(xmi),dimax)
2709 vt2i = 1.49e4*diameter**1.31
2710 vt2s = pvts*rslopeb(i,k,2)*denfac(i,k)
2711 !-------------------------------------------------------------
2712 ! psaci: Accretion of cloud ice by rain [HDC 10]
2713 ! (T<T0: I->S)
2714 !-------------------------------------------------------------
2715 acrfac = 2.*rslope3(i,k,2)+2.*diameter*rslope2(i,k,2) &
2716 +diameter**2*rslope(i,k,2)
2717 psaci(i,k) = pi*qci(i,k,2)*eacrs*n0s*n0sfac(i,k) &
2718 *abs(vt2s-vt2i)*acrfac/4.
2719 endif
2720 endif
2721 !-------------------------------------------------------------
2722 ! psacw: Accretion of cloud water by snow [HL A7] [LFO 24]
2723 ! (T<T0: C->S, and T>=T0: C->R)
2724 !-------------------------------------------------------------
2725 if(qrs(i,k,2).gt.qcrmin.and.qci(i,k,1).gt.qmin) then
2726 psacw(i,k) = min(pacrc*n0sfac(i,k)*rslope3(i,k,2) &
2727 *rslopeb(i,k,2)*qci(i,k,1)*denfac(i,k) &
2728 ! ,qci(i,k,1)/dtcld)
2729 ,qci(i,k,1)*rdtcld)
2730 endif
2731 if(supcol .gt. 0) then
2732 !-------------------------------------------------------------
2733 ! pidep: Deposition/Sublimation rate of ice [HDC 9]
2734 ! (T<T0: V->I or I->V)
2735 !-------------------------------------------------------------
2736 if(qci(i,k,2).gt.0.and.ifsat.ne.1) then
2737 xmi = den(i,k)*qci(i,k,2)/xni(i,k)
2738 diameter = dicon * sqrt(xmi)
2739 pidep(i,k) = 4.*diameter*xni(i,k)*(rh(i,k,2)-1.)/work1(i,k,2)
2740 supice = satdt-prevp(i,k)
2741 if(pidep(i,k).lt.0.) then
2742 ! pidep(i,k) = max(max(pidep(i,k),satdt/2),supice)
2743 ! pidep(i,k) = max(pidep(i,k),-qci(i,k,2)/dtcld)
2744 pidep(i,k) = max(max(pidep(i,k),satdt*.5),supice)
2745 pidep(i,k) = max(pidep(i,k),-qci(i,k,2)*rdtcld)
2746 else
2747 ! pidep(i,k) = min(min(pidep(i,k),satdt/2),supice)
2748 pidep(i,k) = min(min(pidep(i,k),satdt*.5),supice)
2749 endif
2750 if(abs(prevp(i,k)+pidep(i,k)).ge.abs(satdt)) ifsat = 1
2751 endif
2752 !-------------------------------------------------------------
2753 ! psdep: deposition/sublimation rate of snow [HDC 14]
2754 ! (V->S or S->V)
2755 !-------------------------------------------------------------
2756 if(qrs(i,k,2).gt.0..and.ifsat.ne.1) then
2757 coeres = rslope2(i,k,2)*sqrt(rslope(i,k,2)*rslopeb(i,k,2))
2758 psdep(i,k) = (rh(i,k,2)-1.)*n0sfac(i,k) &
2759 *(precs1*rslope2(i,k,2)+precs2 &
2760 *work2(i,k)*coeres)/work1(i,k,2)
2761 supice = satdt-prevp(i,k)-pidep(i,k)
2762 if(psdep(i,k).lt.0.) then
2763 ! psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)/dtcld)
2764 ! psdep(i,k) = max(max(psdep(i,k),satdt/2),supice)
2765 psdep(i,k) = max(psdep(i,k),-qrs(i,k,2)*rdtcld)
2766 psdep(i,k) = max(max(psdep(i,k),satdt*.5),supice)
2767 else
2768 ! psdep(i,k) = min(min(psdep(i,k),satdt/2),supice)
2769 psdep(i,k) = min(min(psdep(i,k),satdt*.5),supice)
2770 endif
2771 if(abs(prevp(i,k)+pidep(i,k)+psdep(i,k)).ge.abs(satdt)) &
2772 ifsat = 1
2773 endif
2774 !-------------------------------------------------------------
2775 ! pigen: generation(nucleation) of ice from vapor [HL A50] [HDC 7-8]
2776 ! (T<T0: V->I)
2777 !-------------------------------------------------------------
2778 if(supsat.gt.0.and.ifsat.ne.1) then
2779 supice = satdt-prevp(i,k)-pidep(i,k)-psdep(i,k)
2780 xni0 = 1.e3*exp(0.1*supcol)
2781 roqi0 = 4.92e-11*exp(log(xni0)*(1.33))
2782 pigen(i,k) = max(0.,(roqi0/den(i,k)-max(qci(i,k,2),0.)) &
2783 ! /dtcld)
2784 *rdtcld)
2785 pigen(i,k) = min(min(pigen(i,k),satdt),supice)
2786 endif
2787 !
2788 !-------------------------------------------------------------
2789 ! psaut: conversion(aggregation) of ice to snow [HDC 12]
2790 ! (T<T0: I->S)
2791 !-------------------------------------------------------------
2792 if(qci(i,k,2).gt.0.) then
2793 qimax = roqimax/den(i,k)
2794 ! psaut(i,k) = max(0.,(qci(i,k,2)-qimax)/dtcld)
2795 psaut(i,k) = max(0.,(qci(i,k,2)-qimax)*rdtcld)
2796 endif
2797 endif
2798 !-------------------------------------------------------------
2799 ! psevp: Evaporation of melting snow [HL A35] [RH83 A27]
2800 ! (T>T0: S->V)
2801 !-------------------------------------------------------------
2802 if(supcol.lt.0.) then
2803 if(qrs(i,k,2).gt.0..and.rh(i,k,1).lt.1.) &
2804 psevp(i,k) = psdep(i,k)*work1(i,k,2)/work1(i,k,1)
2805 ! psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)/dtcld),0.)
2806 psevp(i,k) = min(max(psevp(i,k),-qrs(i,k,2)*rdtcld),0.)
2807 endif
2808 enddo
2809 enddo
2810 !
2811 !
2812 !----------------------------------------------------------------
2813 ! check mass conservation of generation terms and feedback to the
2814 ! large scale
2815 !
2816 do k = kts, kte
2817 do i = its, ite
2818 if(t(i,k).le.t0c) then
2819 !
2820 ! cloud water
2821 !
2822 value = max(qmin,qci(i,k,1))
2823 source = (praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
2824 if (source.gt.value) then
2825 factor = value/source
2826 praut(i,k) = praut(i,k)*factor
2827 pracw(i,k) = pracw(i,k)*factor
2828 psacw(i,k) = psacw(i,k)*factor
2829 endif
2830 !
2831 ! cloud ice
2832 !
2833 value = max(qmin,qci(i,k,2))
2834 source = (psaut(i,k)+psaci(i,k)-pigen(i,k)-pidep(i,k))*dtcld
2835 if (source.gt.value) then
2836 factor = value/source
2837 psaut(i,k) = psaut(i,k)*factor
2838 psaci(i,k) = psaci(i,k)*factor
2839 pigen(i,k) = pigen(i,k)*factor
2840 pidep(i,k) = pidep(i,k)*factor
2841 endif
2842 !
2843 ! rain
2844 !
2845 !
2846 value = max(qmin,qrs(i,k,1))
2847 source = (-praut(i,k)-pracw(i,k)-prevp(i,k))*dtcld
2848 if (source.gt.value) then
2849 factor = value/source
2850 praut(i,k) = praut(i,k)*factor
2851 pracw(i,k) = pracw(i,k)*factor
2852 prevp(i,k) = prevp(i,k)*factor
2853 endif
2854 !
2855 ! snow
2856 !
2857 value = max(qmin,qrs(i,k,2))
2858 source = (-psdep(i,k)-psaut(i,k)-psaci(i,k)-psacw(i,k))*dtcld
2859 if (source.gt.value) then
2860 factor = value/source
2861 psdep(i,k) = psdep(i,k)*factor
2862 psaut(i,k) = psaut(i,k)*factor
2863 psaci(i,k) = psaci(i,k)*factor
2864 psacw(i,k) = psacw(i,k)*factor
2865 endif
2866 !
2867 work2(i,k)=-(prevp(i,k)+psdep(i,k)+pigen(i,k)+pidep(i,k))
2868 ! update
2869 q(i,k) = q(i,k)+work2(i,k)*dtcld
2870 qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
2871 +psacw(i,k))*dtcld,0.)
2872 qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
2873 +prevp(i,k))*dtcld,0.)
2874 qci(i,k,2) = max(qci(i,k,2)-(psaut(i,k)+psaci(i,k) &
2875 -pigen(i,k)-pidep(i,k))*dtcld,0.)
2876 qrs(i,k,2) = max(qrs(i,k,2)+(psdep(i,k)+psaut(i,k) &
2877 +psaci(i,k)+psacw(i,k))*dtcld,0.)
2878 xlf = xls-xl(i,k)
2879 xlwork2 = -xls*(psdep(i,k)+pidep(i,k)+pigen(i,k)) &
2880 -xl(i,k)*prevp(i,k)-xlf*psacw(i,k)
2881 t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
2882 else
2883 !
2884 ! cloud water
2885 !
2886 value = max(qmin,qci(i,k,1))
2887 source=(praut(i,k)+pracw(i,k)+psacw(i,k))*dtcld
2888 if (source.gt.value) then
2889 factor = value/source
2890 praut(i,k) = praut(i,k)*factor
2891 pracw(i,k) = pracw(i,k)*factor
2892 psacw(i,k) = psacw(i,k)*factor
2893 endif
2894 !
2895 ! rain
2896 !
2897 value = max(qmin,qrs(i,k,1))
2898 source = (-praut(i,k)-pracw(i,k)-prevp(i,k)-psacw(i,k))*dtcld
2899 if (source.gt.value) then
2900 factor = value/source
2901 praut(i,k) = praut(i,k)*factor
2902 pracw(i,k) = pracw(i,k)*factor
2903 prevp(i,k) = prevp(i,k)*factor
2904 psacw(i,k) = psacw(i,k)*factor
2905 endif
2906 !
2907 ! snow
2908 !
2909 value = max(qcrmin,qrs(i,k,2))
2910 source=(-psevp(i,k))*dtcld
2911 if (source.gt.value) then
2912 factor = value/source
2913 psevp(i,k) = psevp(i,k)*factor
2914 endif
2915 work2(i,k)=-(prevp(i,k)+psevp(i,k))
2916 ! update
2917 q(i,k) = q(i,k)+work2(i,k)*dtcld
2918 qci(i,k,1) = max(qci(i,k,1)-(praut(i,k)+pracw(i,k) &
2919 +psacw(i,k))*dtcld,0.)
2920 qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,k)+pracw(i,k) &
2921 +prevp(i,k) +psacw(i,k))*dtcld,0.)
2922 qrs(i,k,2) = max(qrs(i,k,2)+psevp(i,k)*dtcld,0.)
2923 xlf = xls-xl(i,k)
2924 xlwork2 = -xl(i,k)*(prevp(i,k)+psevp(i,k))
2925 t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
2926 endif
2927 enddo
2928 enddo
2929 !
2930 ! Inline expansion for fpvs
2931 ! qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
2932 ! qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
2933 hsub = xls
2934 hvap = xlv0
2935 cvap = cpv
2936 ttp=t0c+0.01
2937 dldt=cvap-cliq
2938 xa=-dldt/rv
2939 xb=xa+hvap/(rv*ttp)
2940 dldti=cvap-cice
2941 xai=-dldti/rv
2942 xbi=xai+hsub/(rv*ttp)
2943 do k = kts, kte
2944 do i = its, ite
2945 tr=ttp/t(i,k)
2946 logtr = log(tr)
2947 qs(i,k,1)=psat*exp(logtr*(xa)+xb*(1.-tr))
2948 qs(i,k,1) = ep2 * qs(i,k,1) / (p(i,k) - qs(i,k,1))
2949 qs(i,k,1) = max(qs(i,k,1),qmin)
2950 enddo
2951 enddo
2952 !
2953 !----------------------------------------------------------------
2954 ! pcond: condensational/evaporational rate of cloud water [HL A46] [RH83 A6]
2955 ! if there exists additional water vapor condensated/if
2956 ! evaporation of cloud water is not enough to remove subsaturation
2957 !
2958 do k = kts, kte
2959 do i = its, ite
2960 ! work1(i,k,1) = conden(t(i,k),q(i,k),qs(i,k,1),xl(i,k),cpm(i,k))
2961 work1(i,k,1) = ((max(q(i,k),qmin)-(qs(i,k,1)))/(1.+(xl(i,k)) &
2962 *(xl(i,k))/(rv*(cpm(i,k)))*(qs(i,k,1)) &
2963 /((t(i,k))*(t(i,k)))))
2964 work2(i,k) = qci(i,k,1)+work1(i,k,1)
2965 pcond(i,k) = min(max(work1(i,k,1)/dtcld,0.),max(q(i,k),0.)/dtcld)
2966 if(qci(i,k,1).gt.0..and.work1(i,k,1).lt.0.) &
2967 pcond(i,k) = max(work1(i,k,1),-qci(i,k,1))/dtcld
2968 q(i,k) = q(i,k)-pcond(i,k)*dtcld
2969 qci(i,k,1) = max(qci(i,k,1)+pcond(i,k)*dtcld,0.)
2970 t(i,k) = t(i,k)+pcond(i,k)*xl(i,k)/cpm(i,k)*dtcld
2971 enddo
2972 enddo
2973 !
2974 !
2975 !----------------------------------------------------------------
2976 ! padding for small values
2977 !
2978 do k = kts, kte
2979 do i = its, ite
2980 if(qci(i,k,1).le.qmin) qci(i,k,1) = 0.0
2981 if(qci(i,k,2).le.qmin) qci(i,k,2) = 0.0
2982 enddo
2983 enddo
2984 enddo ! big loops
2985 END SUBROUTINE wsm52d
2987 #endif
2989 ! ...................................................................
2990 REAL FUNCTION rgmma(x)
2991 !-------------------------------------------------------------------
2992 IMPLICIT NONE
2993 !-------------------------------------------------------------------
2994 ! rgmma function: use infinite product form
2995 REAL :: euler
2996 PARAMETER (euler=0.577215664901532)
2997 REAL :: x, y
2998 INTEGER :: i
2999 if(x.eq.1.)then
3000 rgmma=0.
3001 else
3002 rgmma=x*exp(euler*x)
3003 do i=1,10000
3004 y=float(i)
3005 rgmma=rgmma*(1.000+x/y)*exp(-x/y)
3006 enddo
3007 rgmma=1./rgmma
3008 endif
3009 END FUNCTION rgmma
3010 !
3011 !--------------------------------------------------------------------------
3012 REAL FUNCTION fpvs(t,ice,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c)
3013 !--------------------------------------------------------------------------
3014 IMPLICIT NONE
3015 !--------------------------------------------------------------------------
3016 REAL t,rd,rv,cvap,cliq,cice,hvap,hsub,psat,t0c,dldt,xa,xb,dldti, &
3017 xai,xbi,ttp,tr
3018 INTEGER ice
3019 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3020 ttp=t0c+0.01
3021 dldt=cvap-cliq
3022 xa=-dldt/rv
3023 xb=xa+hvap/(rv*ttp)
3024 dldti=cvap-cice
3025 xai=-dldti/rv
3026 xbi=xai+hsub/(rv*ttp)
3027 tr=ttp/t
3028 if(t.lt.ttp.and.ice.eq.1) then
3029 fpvs=psat*exp(log(tr)*(xai))*exp(xbi*(1.-tr))
3030 else
3031 fpvs=psat*exp(log(tr)*(xa))*exp(xb*(1.-tr))
3032 endif
3033 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3034 END FUNCTION fpvs
3035 !-------------------------------------------------------------------
3036 SUBROUTINE wsm5init(den0,denr,dens,cl,cpv,allowed_to_read)
3037 !-------------------------------------------------------------------
3038 IMPLICIT NONE
3039 !-------------------------------------------------------------------
3040 !.... constants which may not be tunable
3041 REAL, INTENT(IN) :: den0,denr,dens,cl,cpv
3042 LOGICAL, INTENT(IN) :: allowed_to_read
3043 REAL :: pi
3044 !
3045 pi = 4.*atan(1.)
3046 xlv1 = cl-cpv
3047 !
3048 qc0 = 4./3.*pi*denr*r0**3*xncr/den0 ! 0.419e-3 -- .61e-3
3049 qck1 = .104*9.8*peaut/(xncr*denr)**(1./3.)/xmyu*den0**(4./3.) ! 7.03
3050 !
3051 bvtr1 = 1.+bvtr
3052 bvtr2 = 2.5+.5*bvtr
3053 bvtr3 = 3.+bvtr
3054 bvtr4 = 4.+bvtr
3055 g1pbr = rgmma(bvtr1)
3056 g3pbr = rgmma(bvtr3)
3057 g4pbr = rgmma(bvtr4) ! 17.837825
3058 g5pbro2 = rgmma(bvtr2) ! 1.8273
3059 pvtr = avtr*g4pbr/6.
3060 eacrr = 1.0
3061 pacrr = pi*n0r*avtr*g3pbr*.25*eacrr
3062 precr1 = 2.*pi*n0r*.78
3063 precr2 = 2.*pi*n0r*.31*avtr**.5*g5pbro2
3064 xmmax = (dimax/dicon)**2
3065 roqimax = 2.08e22*dimax**8
3066 !
3067 bvts1 = 1.+bvts
3068 bvts2 = 2.5+.5*bvts
3069 bvts3 = 3.+bvts
3070 bvts4 = 4.+bvts
3071 g1pbs = rgmma(bvts1) !.8875
3072 g3pbs = rgmma(bvts3)
3073 g4pbs = rgmma(bvts4) ! 12.0786
3074 g5pbso2 = rgmma(bvts2)
3075 pvts = avts*g4pbs/6.
3076 pacrs = pi*n0s*avts*g3pbs*.25
3077 precs1 = 4.*n0s*.65
3078 precs2 = 4.*n0s*.44*avts**.5*g5pbso2
3079 pidn0r = pi*denr*n0r
3080 pidn0s = pi*dens*n0s
3081 pacrc = pi*n0s*avts*g3pbs*.25*eacrc
3082 !
3083 rslopermax = 1./lamdarmax
3084 rslopesmax = 1./lamdasmax
3085 rsloperbmax = rslopermax ** bvtr
3086 rslopesbmax = rslopesmax ** bvts
3087 rsloper2max = rslopermax * rslopermax
3088 rslopes2max = rslopesmax * rslopesmax
3089 rsloper3max = rsloper2max * rslopermax
3090 rslopes3max = rslopes2max * rslopesmax
3091 !
3092 END SUBROUTINE wsm5init
3093 END MODULE module_mp_wsm5