?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Update version info for release v4.6.1 (#2122)
[WRF.git] / phys / mic-wsm5-3-5-code.h
blob87149c0d5c434a53cf8cd5c1af9080e50e5ab5d5
1 #define WSM_NO_CONDITIONAL_IN_VECTOR
2
3 SUBROUTINE wsm52D(t, q &
4 ,qci, qrs, den, p, delz &
5 ,delt,g, cpd, cpv, rd, rv, t0c &
6 ,ep1, ep2, qmin &
7 ,XLS, XLV0, XLF0, den0, denr &
8 ,cliq,cice,psat &
9 ,lon,lat &
10 ,rain,rainncv &
11 ,sr &
12 ,snow,snowncv &
13 ,nx0,nk0,irestrict &
14 ,doit,kts,kte )
15
16 !-------------------------------------------------------------------
17 IMPLICIT NONE
18 !-------------------------------------------------------------------
19 !
20 ! This code is a 5-class mixed ice microphyiscs scheme (WSM5) of the
21 ! Single-Moment MicroPhyiscs (WSMMP). The WSMMP assumes that ice nuclei
22 ! number concentration is a function of temperature, and seperate assumption
23 ! is developed, in which ice crystal number concentration is a function
24 ! of ice amount. A theoretical background of the ice-microphysics and related
25 ! processes in the WSMMPs are described in Hong et al. (2004).
26 ! Production terms in the WSM6 scheme are described in Hong and Lim (2006).
27 ! All units are in m.k.s. and source/sink terms in kgkg-1s-1.
28 !
29 ! WSM5 cloud scheme
30 !
31 ! Coded by Song-You Hong (Yonsei Univ.)
32 ! Jimy Dudhia (NCAR) and Shu-Hua Chen (UC Davis)
33 ! Summer 2002
34 !
35 ! Implemented by Song-You Hong (Yonsei Univ.) and Jimy Dudhia (NCAR)
36 ! Summer 2003
37 !
38 ! Reference) Hong, Dudhia, Chen (HDC, 2004) Mon. Wea. Rev.
39 ! Rutledge, Hobbs (RH83, 1983) J. Atmos. Sci.
40 ! Hong and Lim (HL, 2006) J. Korean Meteor. Soc.
41 !
42 INTEGER, INTENT(IN ) :: nx0,nk0,irestrict &
43 ,lon,lat,kts,kte
44
45 #define nx CHUNK
46 #define its 1
47 #define ite nx
48 #define ims 1
49 #define ime nx
50 #define kms kts
51 #define kme kte
52
53 REAL, DIMENSION( its:ite , kts:kte ), &
54 INTENT(INOUT) :: &
55 t
56 REAL, DIMENSION( its:ite , kts:kte, 2 ), &
57 INTENT(INOUT) :: &
58 qci, &
59 qrs
60 REAL, DIMENSION( ims:ime , kms:kme ), &
61 INTENT(INOUT) :: &
62 q
63 REAL, DIMENSION( ims:ime , kms:kme ), &
64 INTENT(IN ) :: &
65 den, &
66 p, &
67 delz
68 REAL, INTENT(IN ) :: delt, &
69 g, &
70 cpd, &
71 cpv, &
72 t0c, &
73 den0, &
74 rd, &
75 rv, &
76 ep1, &
77 ep2, &
78 qmin, &
79 XLS, &
80 XLV0, &
81 XLF0, &
82 cliq, &
83 cice, &
84 psat, &
85 denr
86 REAL, DIMENSION( ims:ime ), &
87 INTENT(INOUT) :: rain, &
88 rainncv, &
89 sr
90 REAL, DIMENSION( ims:ime ), OPTIONAL, &
91 INTENT(INOUT) :: snow, &
92 snowncv
93
94 LOGICAL, INTENT(IN) :: doit ! added for conformity with standalone driver
95 ! LOCAL VAR
96 REAL, DIMENSION( its:ite , kts:kte , 2) :: &
97 rh, &
98 qs, &
99 falk, &
100 fall
101
102 REAL, DIMENSION( its:ite , 1 , 2) :: &
103 work1
104 REAL, DIMENSION( its:ite , 1 ) :: &
105 work2
106
107 REAL, DIMENSION( its:ite , 2) :: &
108 rslope_v, &
109 rslope2_v, &
110 rslope3_v, &
111 rslopeb_v
112
113
114
115 REAL, DIMENSION( its:ite , kts:kte ) :: &
116 falkc, &
117 fallc, &
118 xl, &
119 cpm, &
120 denfac, &
121 xni, &
122 denqrs, &
123 denqci, &
124 n0sfac, &
125 workrs, &
126 work1c, &
127 work2c
128 REAL, DIMENSION( its:ite ) :: &
129 delqrs
130 REAL, DIMENSION( its:ite , 1 ) :: &
131 pigen, &
132 pidep, &
133 psdep, &
134 praut, &
135 psaut, &
136 prevp, &
137 psevp, &
138 pracw, &
139 psacw, &
140 psaci, &
141 pcond, &
142 psmlt
143
144 INTEGER, DIMENSION( its:ite ) :: &
145 mstep, &
146 numdt
147 REAL, DIMENSION(its:ite) :: tstepsnow
148 REAL, DIMENSION(its:ite) :: rmstep
149 REAL dtcldden, rdelz, rdtcld
150 LOGICAL, DIMENSION( its:ite ) :: flgcld
151 #ifdef WSM_NO_CONDITIONAL_IN_VECTOR
152 REAL, DIMENSION(its:ite) :: xal, xbl
153 #endif
154 REAL :: &
155 cpmcal, xlcal, diffus, &
156 viscos, xka, venfac, conden, diffac, &
157 x, y, z, a, b, c, d, e, &
158 qdt, holdrr, holdrs, supcol, supcolt, pvt, &
159 coeres, supsat, dtcld, xmi, eacrs, satdt, &
160 vt2i,vt2s,acrfac, &
161 qimax, diameter, xni0, roqi0, &
162 xlwork2, factor, source, &
163 value, xlf, pfrzdtc, pfrzdtr, supice, holdc, holdci
164
165 REAL, DIMENSION(CHUNK) :: fallsum, fallsum_qsi ! convert to vector
166 REAL, DIMENSION(CHUNK) :: supsat_v,satdt_v,coeres_v
167
168 ! variables for optimization
169 REAL, DIMENSION( its:ite ) :: tvec1
170 REAL :: temp
171 INTEGER :: i, j, k, mstepmax, &
172 iprt, latd, lond, loop, loops, ifsat, n, idim, kdim
173 ! Temporaries used for inlining fpvs function, and other vector stuff
174 REAL :: dldti, xb, xai, xbi, xa, hvap, cvap, hsub, dldt, ttp
175 REAL :: tr_v(its:ite),logtr_v(its:ite),supcol_v(its:ite),supcolt_v(its:ite),xlf_v(its:ite),temp_v(its:ite)
176 REAL :: diameter_v(CHUNK),supice_v(CHUNK)
177 INTEGER :: ifsat_v(CHUNK)
178 ! mask variable
179 LOGICAL*4 :: lmask(CHUNK)
180 !
181
182 # define LAMDAR(x,y) sqrt(sqrt(pidn0r/((x)*(y))))
183 # define LAMDAS(x,y,z) sqrt(sqrt(pidn0s*(z)/((x)*(y))))
184
185
186 !
187 !=================================================================
188 ! compute internal functions
189 !
190 #define CPMCAL(x) (cpd*(1.-max(x,qmin))+max(x,qmin)*cpv)
191 #define XLCAL(x) (xlv0-xlv1*((x)-t0c))
192 cpmcal(x) = cpd*(1.-max(x,qmin))+max(x,qmin)*cpv
193 xlcal(x) = xlv0-xlv1*(x-t0c)
194 !----------------------------------------------------------------
195 ! diffus: diffusion coefficient of the water vapor
196 ! viscos: kinematic viscosity(m2s-1)
197 ! diffus(x,y) = 8.794e-5 * exp(log(x)*(1.81)) / y ! 8.794e-5*x**1.81/y
198 ! viscos(x,y) = 1.496e-6 * (x*sqrt(x)) /(x+120.)/y ! 1.496e-6*x**1.5/(x+120.)/y
199 ! xka(x,y) = 1.414e3*viscos(x,y)*y
200 ! diffac(a,b,c,d,e) = d*a*a/(xka(c,d)*rv*c*c)+1./(e*diffus(c,b))
201 ! venfac(a,b,c) = exp(log((viscos(b,c)/diffus(b,a)))*((.3333333))) &
202 ! /sqrt(viscos(b,c))*sqrt(sqrt(den0/c))
203 ! conden(a,b,c,d,e) = (max(b,qmin)-c)/(1.+d*d/(rv*e)*c/(a*a))
204 !
205 !----------------------------------------------------------------
206 !DIR$ ASSUME_ALIGNED t:64
207 !DIR$ ASSUME_ALIGNED qci:64
208 !DIR$ ASSUME_ALIGNED qrs:64
209 !DIR$ ASSUME_ALIGNED q:64
210 !DIR$ ASSUME_ALIGNED den:64
211 !DIR$ ASSUME_ALIGNED p:64
212 !DIR$ ASSUME_ALIGNED delz:64
213 !DIR$ ASSUME_ALIGNED rain:64
214 !DIR$ ASSUME_ALIGNED rainncv:64
215 !DIR$ ASSUME_ALIGNED sr:64
216 !DIR$ ASSUME_ALIGNED snow:64
217 !DIR$ ASSUME_ALIGNED snowncv:64
218
219 if ( irestrict .le. 0 .OR. .NOT. doit ) return
220 idim = ite-its+1
221 kdim = kte-kts+1
222 lmask = .FALSE.
223 do i = its, min(irestrict,ite)
224 lmask(i) = .TRUE.
225 enddo
226 !
227 !----------------------------------------------------------------
228 ! paddint 0 for negative values generated by dynamics
229 !
230 do k = kts, kte
231 WHERE(lmask)
232 qci(:,k,1) = max(qci(:,k,1),0.0)
233 qrs(:,k,1) = max(qrs(:,k,1),0.0)
234 qci(:,k,2) = max(qci(:,k,2),0.0)
235 qrs(:,k,2) = max(qrs(:,k,2),0.0)
236 ENDWHERE
237 enddo
238 !
239 !----------------------------------------------------------------
240 ! latent heat for phase changes and heat capacity. neglect the
241 ! changes during microphysical process calculation
242 ! emanuel(1994)
243 !
244 do k = kts, kte
245 WHERE(lmask)
246 cpm(:,k) = CPMCAL(q(:,k))
247 xl(:,k) = XLCAL(t(:,k))
248 ENDWHERE
249 enddo
250 !
251 !----------------------------------------------------------------
252 ! initialize the surface rain, snow
253 !
254 WHERE(lmask)
255 rainncv(:) = 0.
256 sr(:) = 0.
257 tstepsnow(:) = 0.
258 ENDWHERE
259 IF(PRESENT (snowncv) .AND. PRESENT (snow)) THEN
260 WHERE( lmask ) snowncv(:) = 0.
261 ENDIF
262 !
263 !----------------------------------------------------------------
264 ! compute the minor time steps.
265 !
266 loops = max(nint(delt/dtcldcr),1)
267 dtcld = delt/loops
268 if(delt.le.dtcldcr) dtcld = delt
269 !
270 do loop = 1,loops
271 !
272 !----------------------------------------------------------------
273 ! initialize the large scale variables
274 !
275 WHERE (lmask)
276 mstep(:) = 1
277 flgcld(:) = .true.
278 ENDWHERE
279
280 ! this seems to be needed for the standalone version to agree with its reference
281 do k = kts, kte
282 CALL VREC( tvec1(its), den(its,k), ite-its+1)
283 WHERE( lmask )
284 tvec1(:) = tvec1(:)*den0
285 ENDWHERE
286 CALL VSQRT( denfac(its,k), tvec1(its), ite-its+1)
287 enddo
288
289 !
290 ! Inline expansion for fpvs
291 ! qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
292 ! qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
293 hsub = xls
294 hvap = xlv0
295 cvap = cpv
296 ttp=t0c+0.01
297 dldt=cvap-cliq
298 xa=-dldt/rv
299 xb=xa+hvap/(rv*ttp)
300 dldti=cvap-cice
301 xai=-dldti/rv
302 xbi=xai+hsub/(rv*ttp)
303 ! this is for compilers where the conditional inhibits vectorization
304 #ifdef WSM_NO_CONDITIONAL_IN_VECTOR
305 do k = kts, kte
306 WHERE( lmask )
307 WHERE(t(:,k).lt.ttp)
308 xal(:) = xai
309 xbl(:) = xbi
310 ELSEWHERE
311 xal(:) = xa
312 xbl(:) = xb
313 ENDWHERE
314 tr_v=ttp/t(:,k)
315 logtr_v=log(tr_v)
316 qs(:,k,1)=psat*exp(logtr_v*(xa)+xb*(1.-tr_v))
317 qs(:,k,1) = min(qs(:,k,1),0.99*p(:,k))
318 qs(:,k,1) = ep2 * qs(:,k,1) / (p(:,k) - qs(:,k,1))
319 qs(:,k,1) = max(qs(:,k,1),qmin)
320 rh(:,k,1) = max(q(:,k) / qs(:,k,1),qmin)
321 qs(:,k,2)=psat*exp(logtr_v*(xal(:))+xbl(:)*(1.-tr_v))
322 qs(:,k,2) = min(qs(:,k,2),0.99*p(:,k))
323 qs(:,k,2) = ep2 * qs(:,k,2) / (p(:,k) - qs(:,k,2))
324 qs(:,k,2) = max(qs(:,k,2),qmin)
325 rh(:,k,2) = max(q(:,k) / qs(:,k,2),qmin)
326 ENDWHERE
327 enddo
328 #else
329 Bad --- XEON_OPTIMIZED VERSION NEEDS WSM5_NO_CONDITIONAL_IN_VECTOR
330 #endif
331
332 !
333 !----------------------------------------------------------------
334 ! initialize the variables for microphysical physics
335 !
336 !
337 do k = kts, kte
338 do i = its, min(irestrict,ite)
339 falk(i,k,1) = 0.
340 falk(i,k,2) = 0.
341 fall(i,k,1) = 0.
342 fall(i,k,2) = 0.
343 fallc(i,k) = 0.
344 falkc(i,k) = 0.
345 xni(i,k) = 1.e3
346 enddo
347 enddo
348 !-------------------------------------------------------------
349 ! Ni: ice crystal number concentraiton [HDC 5c]
350 !-------------------------------------------------------------
351 !jm note reuse of tr_v as temporary
352 do k = kts, kte
353 WHERE( lmask )
354 tr_v = (den(:,k)*max(qci(:,k,2),qmin))
355 tr_v = sqrt(sqrt(tr_v*tr_v*tr_v))
356 xni(:,k) = min(max(5.38e7*tr_v,1.e3),1.e6)
357 ENDWHERE
358 enddo
359 !
360 !----------------------------------------------------------------
361 ! compute the fallout term:
362 ! first, vertical terminal velosity for minor loops
363 !----------------------------------------------------------------
364 !
365 #if 1
366 do k = kts, kte
367 WHERE(lmask)
368 !---------------------------------------------------------------
369 ! n0s: Intercept parameter for snow [m-4] [HDC 6]
370 !---------------------------------------------------------------
371 WHERE (qrs(:,k,1).le.0.0)
372 workrs(:,k) = 0.0
373 ELSEWHERE (qrs(:,k,1).le.qcrmin)
374 workrs(:,k) = pvtr*rsloperbmax*denfac(:,k)
375 ELSEWHERE
376 ! (rslopeb ( rslope ) )
377 workrs(:,k) = pvtr*( exp(log( 1./LAMDAR(qrs(:,k,1),den(:,k)) )*(bvtr)) )*denfac(:,k)
378 ENDWHERE
379 ENDWHERE
380 enddo
381
382 qrs(:,:,1) = den*qrs(:,:,1)
383 call nislfv_rain_plm(idim,kdim,den,denfac,t,delz,workrs,qrs(:,:,1), &
384 delqrs,dtcld,1,1,irestrict,lon,lat,.true.,1)
385 fall(:,:,1) = qrs(:,:,1)*workrs/delz
386 qrs(:,:,1) = max(qrs(:,:,1)/den,0.)
387 fall(:,1,1) = delqrs/delz(:,1)/dtcld
388
389 do k = kts, kte
390 WHERE(lmask)
391 WHERE (qrs(:,k,2).le.0.0)
392 workrs(:,k) = 0.0
393 ELSEWHERE (qrs(:,k,2).le.qcrmin)
394 workrs(:,k) = pvts*rslopesbmax*denfac(:,k)
395 ELSEWHERE
396 workrs(:,k) = pvts*(exp(log(1./ &
397 LAMDAS(qrs(:,k,2),den(:,k), max(min(exp(alpha*(t0c-t(:,k))),n0smax/n0s),1.)) &
398 ) *(bvts)))*denfac(:,k)
399 ENDWHERE
400 ENDWHERE
401 enddo
402 qrs(:,:,2) = den*qrs(:,:,2)
403 call nislfv_rain_plm(idim,kdim,den,denfac,t,delz,workrs,qrs(:,:,2), &
404 delqrs,dtcld,2,1,irestrict,lon,lat,.false.,2)
405 fall(:,:,2) = qrs(:,:,2)*workrs/delz
406 qrs(:,:,2) = max(qrs(:,:,2)/den,0.)
407 fall(:,1,2) = delqrs/delz(:,1)/dtcld
408 #else
409 WHERE( qrs(:,:,1) .le. 0.0 )
410 workr = 0.0
411 ELSEWHERE
412 workr = work1(:,:,1)
413 ENDWHERE
414 denqrs1 = den*qrs(:,:,1)
415 call nislfv_rain_plm(idim,kdim,den,denfac,t,delz,workr,denqrs1, &
416 delqrs1,dtcld,1,1,irestrict,lon,lat,.true.,1)
417 qrs(:,:,1) = max(denqrs1/den,0.)
418 fall(:,:,1) = denqrs1*workr/delz
419 fall(:,1,1) = delqrs1/delz(:,1)/dtcld
420 WHERE( qrs(:,:,2) .le. 0.0 )
421 works = 0.0
422 ELSEWHERE
423 works = work1(:,:,2)
424 ENDWHERE
425 denqrs2 = den*qrs(:,:,2)
426 call nislfv_rain_plm(idim,kdim,den,denfac,t,delz,works,denqrs2, &
427 delqrs2,dtcld,2,1,irestrict,lon,lat,.false.,2)
428 qrs(:,:,2) = max(denqrs2/den,0.)
429 fall(:,:,2) = denqrs2*works/delz
430 fall(:,1,2) = delqrs2/delz(:,1)/dtcld
431 #endif
432
433 !note reuse of tr_v as temporary for coeres
434 xlf = xlf0
435 do k = kte, kts, -1
436 psmlt = 0.
437 WHERE(lmask)
438 supcol_v = t0c-t(:,k)
439 n0sfac(:,k) = max(min(exp(alpha*supcol_v),n0smax/n0s),1.)
440 WHERE(qrs(:,k,2).le.qcrmin)
441 rslope_v(:,2) = rslopesmax
442 rslopeb_v(:,2) = rslopesbmax
443 rslope2_v(:,2) = rslopes2max
444 ELSEWHERE
445 rslope_v(:,2) = 1./LAMDAS(qrs(:,k,2),den(:,k),max(min(exp(alpha*(t0c-t(:,k))),n0smax/n0s),1.))
446 rslopeb_v(:,2) = exp(log(rslope_v(:,2))*(bvts))
447 rslope2_v(:,2) = rslope_v(:,2)*rslope_v(:,2)
448 ENDWHERE
449 WHERE (t(:,k).gt.t0c.and.qrs(:,k,2).gt.0.)
450 !----------------------------------------------------------------
451 ! psmlt: melting of snow [HL A33] [RH83 A25]
452 ! (T>T0: S->R)
453 !----------------------------------------------------------------
454 work2(:,1)= (exp(log(((1.496e-6*((t(:,k))*sqrt(t(:,k))) &
455 /((t(:,k))+120.)/(den(:,k)))/(8.794e-5 &
456 *exp(log(t(:,k))*(1.81))/p(:,k)))) &
457 *((.3333333)))/sqrt((1.496e-6*((t(:,k)) &
458 *sqrt(t(:,k)))/((t(:,k))+120.)/(den(:,k)))) &
459 *sqrt(sqrt(den0/(den(:,k)))))
460 tr_v = rslope2_v(:,2)*sqrt(rslope_v(:,2)*rslopeb_v(:,2))
461 psmlt(:,1) = (1.414e3*(1.496e-6*((t(:,k))*sqrt(t(:,k))) &
462 /((t(:,k))+120.)/(den(:,k)) )*(den(:,k))) &
463 /xlf*(t0c-t(:,k))*pi/2. &
464 *n0sfac(:,k)*(precs1*rslope2_v(:,2)+precs2 &
465 *work2(:,1)*tr_v)
466 psmlt(:,1) = min(max(psmlt(:,1)*dtcld/mstep(:), &
467 -qrs(:,k,2)/mstep(:)),0.)
468 qrs(:,k,2) = qrs(:,k,2) + psmlt(:,1)
469 qrs(:,k,1) = qrs(:,k,1) - psmlt(:,1)
470 t(:,k) = t(:,k) + xlf/cpm(:,k)*psmlt(:,1)
471 ENDWHERE
472 ENDWHERE
473 enddo
474 !---------------------------------------------------------------
475 ! Vice [ms-1] : fallout of ice crystal [HDC 5a]
476 !---------------------------------------------------------------
477 work1c = 0.
478 WHERE (qci(:,:,2).gt.0.)
479 work1c = &
480 1.49e4*exp( &
481 log( &
482 max(min(dicon * sqrt(den*qci(:,:,2)/xni),dimax), 1.e-25) &
483 )*(1.31) &
484 )
485 ENDWHERE
486 denqci = den*qci(:,:,2)
487 call nislfv_rain_plm(idim,kdim,den,denfac,t,delz,work1c,denqci, &
488 delqrs,dtcld,1,0,irestrict,lon,lat,.false.,3)
489 do k = kts, kte
490 WHERE(lmask)
491 qci(:,k,2) = max(denqci(:,k)/den(:,k),0.)
492 ENDWHERE
493 enddo
494 WHERE(lmask)
495 fallc(:,1) = delqrs(:)/delz(:,1)/dtcld
496 ENDWHERE
497 !
498 !----------------------------------------------------------------
499 ! rain (unit is mm/sec;kgm-2s-1: /1000*delt ===> m)==> mm for wrf
500 !
501 fallsum = fall(:,1,1)+fall(:,1,2)+fallc(:,1)
502 fallsum_qsi = fall(:,1,2)+fallc(:,1)
503 WHERE (lmask .and. fallsum.gt.0.)
504 rainncv = fallsum*delz(:,1)/denr*dtcld*1000. + rainncv
505 rain = fallsum*delz(:,1)/denr*dtcld*1000. + rain
506 ENDWHERE
507 WHERE (lmask .and. fallsum_qsi.gt.0.)
508 tstepsnow = fallsum_qsi*delz(:,kts)/denr*dtcld*1000. &
509 +tstepsnow
510 snowncv = fallsum_qsi*delz(:,kts)/denr*dtcld*1000. &
511 +snowncv
512 snow = fallsum_qsi*delz(:,kts)/denr*dtcld*1000. + snow
513 ENDWHERE
514 WHERE (lmask.and.fallsum.gt.0.)sr=tstepsnow/(rainncv+1.e-12)
515 !
516 !---------------------------------------------------------------
517 ! pimlt: instantaneous melting of cloud ice [HL A47] [RH83 A28]
518 ! (T>T0: I->C)
519 !---------------------------------------------------------------
520 do k = kts, kte
521 WHERE(lmask)
522 supcol_v = t0c-t(:,k)
523 xlf_v = xlf0
524 WHERE(supcol_v.ge.0.) xlf_v = xls-xl(:,k)
525 WHERE(supcol_v.lt.0..and.qci(:,k,2).gt.0.)
526 qci(:,k,1) = qci(:,k,1) + qci(:,k,2)
527 t(:,k) = t(:,k) - xlf_v/cpm(:,k)*qci(:,k,2)
528 qci(:,k,2) = 0.
529 ENDWHERE
530 !---------------------------------------------------------------
531 ! pihmf: homogeneous freezing of cloud water below -40c [HL A45]
532 ! (T<-40C: C->I)
533 !---------------------------------------------------------------
534 WHERE(supcol_v.gt.40..and.qci(:,k,1).gt.0.)
535 qci(:,k,2) = qci(:,k,2) + qci(:,k,1)
536 t(:,k) = t(:,k) + xlf_v/cpm(:,k)*qci(:,k,1)
537 qci(:,k,1) = 0.
538 ENDWHERE
539 !---------------------------------------------------------------
540 ! pihtf: heterogeneous freezing of cloud water [HL A44]
541 ! (T0>T>-40C: C->I)
542 !---------------------------------------------------------------
543 !jm note use of tr_v for temporary
544 WHERE(supcol_v.gt.0..and.qci(:,k,1).gt.0)
545 supcolt_v=min(supcol_v,50.)
546 tr_v = min(pfrz1*(exp(pfrz2*supcolt_v)-1.) &
547 *den(:,k)/denr/xncr*qci(:,k,1)*qci(:,k,1)*dtcld,qci(:,k,1))
548 qci(:,k,2) = qci(:,k,2) + tr_v
549 t(:,k) = t(:,k) + xlf_v/cpm(:,k)*tr_v
550 qci(:,k,1) = qci(:,k,1)-tr_v
551 ENDWHERE
552 !---------------------------------------------------------------
553 ! psfrz: freezing of rain water [HL A20] [LFO 45]
554 ! (T<T0, R->S)
555 !---------------------------------------------------------------
556 !jm note use of tr_v
557 WHERE(supcol_v.gt.0..and.qrs(:,k,1).gt.0)
558 supcolt_v=min(supcol_v,50.)
559 WHERE ( qrs(:,k,1).le.qcrmin)
560 temp_v = (rslopermax)
561 ELSEWHERE
562 temp_v = (1./LAMDAR(qrs(:,k,1),den(:,k)))
563 ENDWHERE
564 temp_v = temp_v*temp_v*temp_v*temp_v*temp_v*temp_v*temp_v
565 tr_v = min(20.*(pi*pi)*pfrz1*n0r*denr/den(:,k) &
566 *(exp(pfrz2*supcolt_v)-1.)*temp_v*dtcld, &
567 qrs(:,k,1))
568 qrs(:,k,2) = qrs(:,k,2) + tr_v
569 t(:,k) = t(:,k) + xlf_v/cpm(:,k)*tr_v
570 qrs(:,k,1) = qrs(:,k,1)-tr_v
571 ENDWHERE
572 ENDWHERE
573 enddo
574 !
575 !----------------------------------------------------------------
576 ! update the slope parameters for microphysics computation
577 !
578 !------------------------------------------------------------------
579 ! work1: the thermodynamic term in the denominator associated with
580 ! heat conduction and vapor diffusion
581 ! (ry88, y93, h85)
582 ! work2: parameter associated with the ventilation effects(y93)
583 !
584 rdtcld = 1./dtcld
585 ! big fused k-loop
586 do k = kts, kte
587 do i = its, min(irestrict,ite)
588 prevp(i,1) = 0.
589 psdep(i,1) = 0.
590 praut(i,1) = 0.
591 psaut(i,1) = 0.
592 pracw(i,1) = 0.
593 psaci(i,1) = 0.
594 psacw(i,1) = 0.
595 pigen(i,1) = 0.
596 pidep(i,1) = 0.
597 pcond(i,1) = 0.
598 psevp(i,1) = 0.
599 enddo
600 WHERE(lmask)
601 work1(:,1,1) = (((((den(:,k))*(xl(:,k))*(xl(:,k)))*((t(:,k))+120.) & ! was work1
602 *(den(:,k)))/(1.414e3*(1.496e-6*((t(:,k))*sqrt(t(:,k)))) &
603 *(den(:,k))*(rv*(t(:,k))*(t(:,k))))) &
604 +p(:,k)/((qs(:,k,1))*(8.794e-5*exp(log(t(:,k))*(1.81)))))
605 work1(:,1,2) = ((((den(:,k))*(xls)*(xls))*((t(:,k))+120.)*(den(:,k)))&
606 /(1.414e3*(1.496e-6*((t(:,k))*sqrt(t(:,k))))*(den(:,k)) &
607 *(rv*(t(:,k))*(t(:,k)))) &
608 + p(:,k)/(qs(:,k,2)*(8.794e-5*exp(log(t(:,k))*(1.81)))))
609 work2(:,1) = (exp(.3333333*log(((1.496e-6 * ((t(:,k))*sqrt(t(:,k)))) &
610 *p(:,k))/(((t(:,k))+120.)*den(:,k)*(8.794e-5 &
611 *exp(log(t(:,k))*(1.81))))))*sqrt(sqrt(den0/(den(:,k))))) &
612 /sqrt((1.496e-6*((t(:,k))*sqrt(t(:,k)))) &
613 /(((t(:,k))+120.)*den(:,k)))
614 ENDWHERE
615 !
616 !===============================================================
617 !
618 ! warm rain processes
619 !
620 ! - follows the processes in RH83 and LFO except for autoconcersion
621 !
622 !===============================================================
623 !
624 ! do k = kts, kte
625 WHERE(lmask)
626 WHERE(qrs(:,k,1).le.qcrmin)
627 rslope_v(:,1) = rslopermax
628 rslopeb_v(:,1) = rsloperbmax
629 rslope2_v(:,1) = rsloper2max
630 rslope3_v(:,1) = rsloper3max
631 ELSEWHERE
632 rslope_v(:,1) = 1./LAMDAR(qrs(:,k,1),den(:,k))
633 rslopeb_v(:,1) = exp(log(rslope_v(:,1))*(bvtr))
634 rslope2_v(:,1) = rslope_v(:,1)*rslope_v(:,1)
635 rslope3_v(:,1) = rslope2_v(:,1)*rslope_v(:,1)
636 ENDWHERE
637 supsat_v = max(q(:,k),qmin)-qs(:,k,1)
638 satdt_v = supsat_v/dtcld
639 !---------------------------------------------------------------
640 ! praut: auto conversion rate from cloud to rain [HDC 16]
641 ! (C->R)
642 !---------------------------------------------------------------
643 WHERE (qci(:,k,1).gt.qc0)
644 praut(:,1) = qck1*exp(log(qci(:,k,1))*((7./3.)))
645 praut(:,1) = min(praut(:,1),qci(:,k,1)/dtcld)
646 ENDWHERE
647 !---------------------------------------------------------------
648 ! pracw: accretion of cloud water by rain [HL A40] [LFO 51]
649 ! (C->R)
650 !---------------------------------------------------------------
651
652 WHERE (qrs(:,k,1).gt.qcrmin.and.qci(:,k,1).gt.qmin)
653 pracw(:,1) = min(pacrr*rslope3_v(:,1)*rslopeb_v(:,1) &
654 *qci(:,k,1)*denfac(:,k),qci(:,k,1)/dtcld)
655 ENDWHERE
656 !---------------------------------------------------------------
657 ! prevp: evaporation/condensation rate of rain [HDC 14]
658 ! (V->R or R->V)
659 !---------------------------------------------------------------
660 WHERE(qrs(:,k,1).gt.0.)
661 coeres_v = rslope2_v(:,1)*sqrt(rslope_v(:,1)*rslopeb_v(:,1))
662 prevp(:,1) = (rh(:,k,1)-1.)*(precr1*rslope2_v(:,1) &
663 +precr2*work2(:,1)*coeres_v)/work1(:,1,1)
664 WHERE(prevp(:,1).lt.0.)
665 prevp(:,1) = max(prevp(:,1),-qrs(:,k,1)/dtcld)
666 prevp(:,1) = max(prevp(:,1),satdt_v/2)
667 ELSEWHERE
668 prevp(:,1) = min(prevp(:,1),satdt_v/2)
669 ENDWHERE
670 ENDWHERE
671 ENDWHERE
672
673 !jm fused-K enddo
674 !
675 !===============================================================
676 !
677 ! cold rain processes
678 !
679 ! - follows the revised ice microphysics processes in HDC
680 ! - the processes same as in RH83 and RH84 and LFO behave
681 ! following ice crystal hapits defined in HDC, inclduing
682 ! intercept parameter for snow (n0s), ice crystal number
683 ! concentration (ni), ice nuclei number concentration
684 ! (n0i), ice diameter (d)
685 !
686 !===============================================================
687 !
688 !jm fused-K do k = kts, kte
689 WHERE(lmask)
690 !---------------------------------------------------------------
691 ! n0s: Intercept parameter for snow [m-4] [HDC 6]
692 !---------------------------------------------------------------
693 WHERE(qrs(:,k,2).le.qcrmin)
694 rslope_v(:,2) = rslopesmax
695 rslopeb_v(:,2) = rslopesbmax
696 rslope2_v(:,2) = rslopes2max
697 rslope3_v(:,2) = rslopes3max
698 ELSEWHERE
699 rslope_v(:,2) = 1./LAMDAS(qrs(:,k,2),den(:,k),max(min(exp(alpha*(t0c-t(:,k))),n0smax/n0s),1.))
700 rslopeb_v(:,2) = exp(log(rslope_v(:,2))*(bvts))
701 rslope2_v(:,2) = rslope_v(:,2)*rslope_v(:,2)
702 rslope3_v(:,2) = rslope2_v(:,2)*rslope_v(:,2)
703 ENDWHERE
704 ENDWHERE
705
706 WHERE(lmask)
707 supcol_v = t0c-t(:,k)
708 n0sfac(:,1) = max(min(exp(alpha*supcol_v),n0smax/n0s),1.)
709 supsat_v = max(q(:,k),qmin)-qs(:,k,2)
710 satdt_v = supsat_v/dtcld
711 ifsat_v(:) = 0
712 !-------------------------------------------------------------
713 ! Ni: ice crystal number concentraiton [HDC 5c]
714 !-------------------------------------------------------------
715 temp_v = (den(:,k)*max(qci(:,k,2),qmin))
716 temp_v = sqrt(sqrt(temp_v*temp_v*temp_v))
717 xni(:,1) = min(max(5.38e7*temp_v,1.e3),1.e6)
718 !
719 psaci(:,1) = 0.
720 psacw(:,1) = 0.
721 WHERE( supcol_v .gt. 0 )
722 WHERE(qrs(:,k,2).gt.qcrmin.and.qci(:,k,2).gt.qmin)
723 diameter_v = min(dicon * sqrt(den(:,k)*qci(:,k,2)/xni(:,1)),dimax)
724 !-------------------------------------------------------------
725 ! psaci: Accretion of cloud ice by rain [HDC 10]
726 ! (T<T0: I->S)
727 !-------------------------------------------------------------
728 psaci(:,1) = pi*qci(:,k,2)*(exp(0.07*(-supcol_v)))*n0s*n0sfac(:,1) &
729 *abs((pvts*rslopeb_v(:,2)*denfac(:,k)) &
730 -(1.49e4*diameter_v**1.31))*(2.*rslope3_v(:,2)+2. &
731 *diameter_v*rslope2_v(:,2) &
732 +diameter_v*diameter_v*rslope_v(:,2))/4.
733 ENDWHERE
734 ENDWHERE
735 !-------------------------------------------------------------
736 ! psacw: Accretion of cloud water by snow [HL A7] [LFO 24]
737 ! (T<T0: C->S, and T>=T0: C->R)
738 !-------------------------------------------------------------
739 WHERE(qrs(:,k,2).gt.qcrmin.and.qci(:,k,1).gt.qmin)
740 psacw(:,1) = min(pacrc*n0sfac(:,1)*rslope3_v(:,2) &
741 *rslopeb_v(:,2)*qci(:,k,1)*denfac(:,k) &
742 ,qci(:,k,1)*rdtcld)
743 ENDWHERE
744 pidep(:,1) = 0.
745 psdep(:,1) = 0.
746 WHERE(supcol_v .gt. 0)
747 !-------------------------------------------------------------
748 ! pidep: Deposition/Sublimation rate of ice [HDC 9]
749 ! (T<T0: V->I or I->V)
750 !-------------------------------------------------------------
751 WHERE(qci(:,k,2).gt.0.and.ifsat_v(:).ne.1)
752 diameter_v = dicon * sqrt(den(:,k)*qci(:,k,2)/xni(:,1))
753 pidep(:,1) = 4.*diameter_v*xni(:,1)*(rh(:,k,2)-1.)/work1(:,1,2)
754 supice_v = satdt_v-prevp(:,1)
755 WHERE(pidep(:,1).lt.0.)
756 pidep(:,1) = max(max(pidep(:,1),satdt_v*.5),supice_v)
757 pidep(:,1) = max(pidep(:,1),-qci(:,k,2)*rdtcld)
758 ELSEWHERE
759 pidep(:,1) = min(min(pidep(:,1),satdt_v*.5),supice_v)
760 ENDWHERE
761 WHERE(abs(prevp(:,1)+pidep(:,1)).ge.abs(satdt_v)) ifsat_v(:) = 1
762 ENDWHERE
763 !-------------------------------------------------------------
764 ! psdep: deposition/sublimation rate of snow [HDC 14]
765 ! (V->S or S->V)
766 !-------------------------------------------------------------
767 WHERE(qrs(:,k,2).gt.0..and.ifsat_v(:).ne.1)
768 psdep(:,1) = (rh(:,k,2)-1.)*n0sfac(:,1) &
769 *(precs1*rslope2_v(:,2)+precs2 &
770 *work2(:,1)*(rslope2_v(:,2)*sqrt(rslope_v(:,2)*rslopeb_v(:,2))))/work1(:,1,2)
771 supice_v = satdt_v-prevp(:,1)-pidep(:,1)
772 WHERE(psdep(:,1).lt.0.)
773 psdep(:,1) = max(psdep(:,1),-qrs(:,k,2)*rdtcld)
774 psdep(:,1) = max(max(psdep(:,1),satdt_v*.5),supice_v)
775 ELSEWHERE
776 psdep(:,1) = min(min(psdep(:,1),satdt_v*.5),supice_v)
777 ENDWHERE
778 WHERE(abs(prevp(:,1)+pidep(:,1)+psdep(:,1)).ge.abs(satdt_v)) &
779 ifsat_v(:) = 1
780 ENDWHERE
781 !-------------------------------------------------------------
782 ! pigen: generation(nucleation) of ice from vapor [HL A50] [HDC 7-8]
783 ! (T<T0: V->I)
784 !-------------------------------------------------------------
785 pigen(:,1) = 0.
786 WHERE(supsat_v.gt.0.and.ifsat_v(:).ne.1)
787 supice_v = satdt_v-prevp(:,1)-pidep(:,1)-psdep(:,1)
788 pigen(:,1) = max(0.,((4.92e-11*exp(log(1.e3*exp(0.1*supcol_v)) &
789 *(1.33)))/den(:,k)-max(qci(:,k,2),0.)) &
790 *rdtcld)
791 pigen(:,1) = min(min(pigen(:,1),satdt_v),supice_v)
792 ENDWHERE
793 !
794 !-------------------------------------------------------------
795 ! psaut: conversion(aggregation) of ice to snow [HDC 12]
796 ! (T<T0: I->S)
797 !-------------------------------------------------------------
798 psaut(:,1) = 0.
799 WHERE(qci(:,k,2).gt.0.)
800 psaut(:,1) = max(0.,(qci(:,k,2)-(roqimax/den(:,k)))*rdtcld)
801 ENDWHERE
802 ENDWHERE
803 !-------------------------------------------------------------
804 ! psevp: Evaporation of melting snow [HL A35] [RH83 A27]
805 ! (T>T0: S->V)
806 !-------------------------------------------------------------
807 psevp(:,1) = 0.
808 WHERE(supcol_v.lt.0.)
809 WHERE(qrs(:,k,2).gt.0..and.rh(:,k,1).lt.1.)
810 psevp(:,1) = psdep(:,1)*work1(:,1,2)/work1(:,1,1)
811 ENDWHERE
812 psevp(:,1) = min(max(psevp(:,1),-qrs(:,k,2)*rdtcld),0.)
813 ENDWHERE
814 ENDWHERE
815
816 !JM fuse-K enddo
817 !
818 !
819 !----------------------------------------------------------------
820 ! check mass conservation of generation terms and feedback to the
821 ! large scale
822 !
823
824 !JM fuse-K do k = kts, kte
825 do i = its, min(irestrict,ite)
826 if(t(i,k).le.t0c) then
827 !
828 ! cloud water
829 !
830 value = max(qmin,qci(i,k,1))
831 source = (praut(i,1)+pracw(i,1)+psacw(i,1))*dtcld
832 if (source.gt.value) then
833 factor = value/source
834 praut(i,1) = praut(i,1)*factor
835 pracw(i,1) = pracw(i,1)*factor
836 psacw(i,1) = psacw(i,1)*factor
837 endif
838 !
839 ! cloud ice
840 !
841 value = max(qmin,qci(i,k,2))
842 source = (psaut(i,1)+psaci(i,1)-pigen(i,1)-pidep(i,1))*dtcld
843 if (source.gt.value) then
844 factor = value/source
845 psaut(i,1) = psaut(i,1)*factor
846 psaci(i,1) = psaci(i,1)*factor
847 pigen(i,1) = pigen(i,1)*factor
848 pidep(i,1) = pidep(i,1)*factor
849 endif
850 !
851 ! rain
852 !
853 !
854 value = max(qmin,qrs(i,k,1))
855 source = (-praut(i,1)-pracw(i,1)-prevp(i,1))*dtcld
856 if (source.gt.value) then
857 factor = value/source
858 praut(i,1) = praut(i,1)*factor
859 pracw(i,1) = pracw(i,1)*factor
860 prevp(i,1) = prevp(i,1)*factor
861 endif
862 !
863 ! snow
864 !
865 value = max(qmin,qrs(i,k,2))
866 source = (-psdep(i,1)-psaut(i,1)-psaci(i,1)-psacw(i,1))*dtcld
867 if (source.gt.value) then
868 factor = value/source
869 psdep(i,1) = psdep(i,1)*factor
870 psaut(i,1) = psaut(i,1)*factor
871 psaci(i,1) = psaci(i,1)*factor
872 psacw(i,1) = psacw(i,1)*factor
873 endif
874 !
875 work2(i,1)=-(prevp(i,1)+psdep(i,1)+pigen(i,1)+pidep(i,1))
876 ! update
877 q(i,k) = q(i,k)+work2(i,1)*dtcld
878 qci(i,k,1) = max(qci(i,k,1)-(praut(i,1)+pracw(i,1) &
879 +psacw(i,1))*dtcld,0.)
880 qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,1)+pracw(i,1) &
881 +prevp(i,1))*dtcld,0.)
882 qci(i,k,2) = max(qci(i,k,2)-(psaut(i,1)+psaci(i,1) &
883 -pigen(i,1)-pidep(i,1))*dtcld,0.)
884 qrs(i,k,2) = max(qrs(i,k,2)+(psdep(i,1)+psaut(i,1) &
885 +psaci(i,1)+psacw(i,1))*dtcld,0.)
886 xlf = xls-xl(i,k)
887 xlwork2 = -xls*(psdep(i,1)+pidep(i,1)+pigen(i,1)) &
888 -xl(i,k)*prevp(i,1)-xlf*psacw(i,1)
889 t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
890 else
891 !
892 ! cloud water
893 !
894 value = max(qmin,qci(i,k,1))
895 source=(praut(i,1)+pracw(i,1)+psacw(i,1))*dtcld
896 if (source.gt.value) then
897 factor = value/source
898 praut(i,1) = praut(i,1)*factor
899 pracw(i,1) = pracw(i,1)*factor
900 psacw(i,1) = psacw(i,1)*factor
901 endif
902 !
903 ! rain
904 !
905 value = max(qmin,qrs(i,k,1))
906 source = (-praut(i,1)-pracw(i,1)-prevp(i,1)-psacw(i,1))*dtcld
907 if (source.gt.value) then
908 factor = value/source
909 praut(i,1) = praut(i,1)*factor
910 pracw(i,1) = pracw(i,1)*factor
911 prevp(i,1) = prevp(i,1)*factor
912 psacw(i,1) = psacw(i,1)*factor
913 endif
914 !
915 ! snow
916 !
917 value = max(qcrmin,qrs(i,k,2))
918 source=(-psevp(i,1))*dtcld
919 if (source.gt.value) then
920 factor = value/source
921 psevp(i,1) = psevp(i,1)*factor
922 endif
923 work2(i,1)=-(prevp(i,1)+psevp(i,1))
924 ! update
925 q(i,k) = q(i,k)+work2(i,1)*dtcld
926 qci(i,k,1) = max(qci(i,k,1)-(praut(i,1)+pracw(i,1) &
927 +psacw(i,1))*dtcld,0.)
928 qrs(i,k,1) = max(qrs(i,k,1)+(praut(i,1)+pracw(i,1) &
929 +prevp(i,1) +psacw(i,1))*dtcld,0.)
930 qrs(i,k,2) = max(qrs(i,k,2)+psevp(i,1)*dtcld,0.)
931 xlf = xls-xl(i,k)
932 xlwork2 = -xl(i,k)*(prevp(i,1)+psevp(i,1))
933 t(i,k) = t(i,k)-xlwork2/cpm(i,k)*dtcld
934
935 endif
936 enddo
937 enddo ! fuse K
938 !
939 ! Inline expansion for fpvs
940 ! qs(i,k,1) = fpvs(t(i,k),0,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
941 ! qs(i,k,2) = fpvs(t(i,k),1,rd,rv,cpv,cliq,cice,xlv0,xls,psat,t0c)
942 hsub = xls
943 hvap = xlv0
944 cvap = cpv
945 ttp=t0c+0.01
946 dldt=cvap-cliq
947 xa=-dldt/rv
948 xb=xa+hvap/(rv*ttp)
949 dldti=cvap-cice
950 xai=-dldti/rv
951 xbi=xai+hsub/(rv*ttp)
952 do k = kts, kte
953 WHERE(lmask)
954 tr_v=ttp/t(:,k)
955 logtr_v = log(tr_v)
956 qs(:,k,1)=psat*exp(logtr_v*(xa)+xb*(1.-tr_v))
957 qs(:,k,1) = min(qs(:,k,1),0.99*p(:,k))
958 qs(:,k,1) = ep2 * qs(:,k,1) / (p(:,k) - qs(:,k,1))
959 qs(:,k,1) = max(qs(:,k,1),qmin)
960 ENDWHERE
961 enddo
962 !
963 !----------------------------------------------------------------
964 ! pcond: condensational/evaporational rate of cloud water [HL A46] [RH83 A6]
965 ! if there exists additional water vapor condensated/if
966 ! evaporation of cloud water is not enough to remove subsaturation
967 !
968 do k = kts, kte
969 do i = its, min(irestrict,ite)
970 work1(i,1,1) = ((max(q(i,k),qmin)-(qs(i,k,1)))/(1.+(xl(i,k)) &
971 *(xl(i,k))/(rv*(cpm(i,k)))*(qs(i,k,1)) &
972 /((t(i,k))*(t(i,k)))))
973 pcond(i,1) = min(max(work1(i,1,1)/dtcld,0.),max(q(i,k),0.)/dtcld)
974 if(qci(i,k,1).gt.0..and.work1(i,1,1).lt.0.) &
975 pcond(i,1) = max(work1(i,1,1),-qci(i,k,1))/dtcld
976 q(i,k) = q(i,k)-pcond(i,1)*dtcld
977 qci(i,k,1) = max(qci(i,k,1)+pcond(i,1)*dtcld,0.)
978 t(i,k) = t(i,k)+pcond(i,1)*xl(i,k)/cpm(i,k)*dtcld
979 enddo
980 enddo
981 !
982 !
983 !----------------------------------------------------------------
984 ! padding for small values
985 !
986 do k = kts, kte
987 do i = its, min(irestrict,ite)
988 if(qci(i,k,1).le.qmin) qci(i,k,1) = 0.0
989 if(qci(i,k,2).le.qmin) qci(i,k,2) = 0.0
990 enddo
991 enddo
992 enddo ! big loops
993
994 END SUBROUTINE wsm52d
995 !------------------------------------------------------------------------------
996 #if 0
997 SUBROUTINE slope_wsm5(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3, &
998 vt,irestrict,kts,kte,lmask)
999 IMPLICIT NONE
1000 INTEGER :: irestrict,kts,kte
1001 REAL, DIMENSION( its:ite , kts:kte,2) :: &
1002 qrs, &
1003 rslope, &
1004 rslopeb, &
1005 rslope2, &
1006 rslope3, &
1007 vt
1008 REAL, DIMENSION( its:ite , kts:kte) :: &
1009 den, &
1010 denfac, &
1011 t
1012 REAL, PARAMETER :: t0c = 273.15
1013 REAL, DIMENSION( its:ite , kts:kte ) :: &
1014 n0sfac
1015 REAL :: lamdar, lamdas, x, y, z, supcol
1016 REAL :: supcol_v(CHUNK)
1017 integer :: i, j, k
1018 LOGICAL*4 lmask(CHUNK)
1019 !----------------------------------------------------------------
1020 ! size distributions: (x=mixing ratio, y=air density):
1021 ! valid for mixing ratio > 1.e-9 kg/kg.
1022 lamdar(x,y)= sqrt(sqrt(pidn0r/(x*y))) ! (pidn0r/(x*y))**.25
1023 lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y))) ! (pidn0s*z/(x*y))**.25
1024 #define LAMDAR(x,y) sqrt(sqrt(pidn0r/((x)*(y))))
1025 #define LAMDAS(x,y,z) sqrt(sqrt(pidn0s*(z)/((x)*(y))))
1026 !
1027 do k = kts, kte
1028 WHERE(lmask)
1029 supcol_v = t0c-t(:,k)
1030 !---------------------------------------------------------------
1031 ! n0s: Intercept parameter for snow [m-4] [HDC 6]
1032 !---------------------------------------------------------------
1033 n0sfac(:,k) = max(min(exp(alpha*supcol_v),n0smax/n0s),1.)
1034 WHERE(qrs(:,k,1).le.qcrmin)
1035 rslope(:,k,1) = rslopermax
1036 rslopeb(:,k,1) = rsloperbmax
1037 rslope2(:,k,1) = rsloper2max
1038 rslope3(:,k,1) = rsloper3max
1039 ELSEWHERE
1040 rslope(:,k,1) = 1./LAMDAR(qrs(:,k,1),den(:,k))
1041 rslopeb(:,k,1) = exp(log(rslope(:,k,1))*(bvtr))
1042 rslope2(:,k,1) = rslope(:,k,1)*rslope(:,k,1)
1043 rslope3(:,k,1) = rslope2(:,k,1)*rslope(:,k,1)
1044 ENDWHERE
1045 WHERE(qrs(:,k,2).le.qcrmin)
1046 rslope(:,k,2) = rslopesmax
1047 rslopeb(:,k,2) = rslopesbmax
1048 rslope2(:,k,2) = rslopes2max
1049 rslope3(:,k,2) = rslopes3max
1050 ELSEWHERE
1051 rslope(:,k,2) = 1./LAMDAS(qrs(:,k,2),den(:,k),n0sfac(:,k))
1052 rslopeb(:,k,2) = exp(log(rslope(:,k,2))*(bvts))
1053 rslope2(:,k,2) = rslope(:,k,2)*rslope(:,k,2)
1054 rslope3(:,k,2) = rslope2(:,k,2)*rslope(:,k,2)
1055 ENDWHERE
1056 WHERE (qrs(:,k,1).le.0.0)
1057 vt(:,k,1) = 0.0
1058 ELSEWHERE
1059 vt(:,k,1) = pvtr*rslopeb(:,k,1)*denfac(:,k)
1060 ENDWHERE
1061 WHERE (qrs(:,k,2).le.0.0)
1062 vt(:,k,2) = 0.0
1063 ELSEWHERE
1064 vt(:,k,2) = pvts*rslopeb(:,k,2)*denfac(:,k)
1065 ENDWHERE
1066 ENDWHERE
1067 enddo
1068 #undef LAMDAR
1069 #undef LAMDAS
1070 END SUBROUTINE slope_wsm5
1071 #endif
1072
1073 #if 0
1074 #undef nx
1075 #undef nk
1076 #undef its
1077 #undef ite
1078 #undef ims
1079 #undef ime
1080 #undef kms
1081 #undef kme
1082 #undef im
1083 #undef km
1084 subroutine slope_rain(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3, &
1085 vt,its,ite,kts,kte)
1086 IMPLICIT NONE
1087 INTEGER :: its,ite, jts,jte, kts,kte
1088 REAL, DIMENSION( its:ite , kts:kte) :: &
1089 qrs, &
1090 rslope, &
1091 rslopeb, &
1092 rslope2, &
1093 rslope3, &
1094 vt, &
1095 den, &
1096 denfac, &
1097 t
1098 REAL, PARAMETER :: t0c = 273.15
1099 REAL, DIMENSION( its:ite , kts:kte ) :: &
1100 n0sfac
1101 REAL :: lamdar, x, y, z, supcol
1102 integer :: i, j, k
1103 !----------------------------------------------------------------
1104 ! size distributions: (x=mixing ratio, y=air density):
1105 ! valid for mixing ratio > 1.e-9 kg/kg.
1106 lamdar(x,y)= sqrt(sqrt(pidn0r/(x*y))) ! (pidn0r/(x*y))**.25
1107 !
1108 do k = kts, kte
1109 do i = its, ite
1110 if(qrs(i,k).le.qcrmin)then
1111 rslope(i,k) = rslopermax
1112 rslopeb(i,k) = rsloperbmax
1113 rslope2(i,k) = rsloper2max
1114 rslope3(i,k) = rsloper3max
1115 else
1116 rslope(i,k) = 1./lamdar(qrs(i,k),den(i,k))
1117 rslopeb(i,k) = rslope(i,k)**bvtr
1118 rslope2(i,k) = rslope(i,k)*rslope(i,k)
1119 rslope3(i,k) = rslope2(i,k)*rslope(i,k)
1120 endif
1121 vt(i,k) = pvtr*rslopeb(i,k)*denfac(i,k)
1122 if(qrs(i,k).le.0.0) vt(i,k) = 0.0
1123 enddo
1124 enddo
1125 END subroutine slope_rain
1126 !------------------------------------------------------------------------------
1127 subroutine slope_snow(qrs,den,denfac,t,rslope,rslopeb,rslope2,rslope3, &
1128 vt,its,ite,kts,kte)
1129 IMPLICIT NONE
1130 INTEGER :: its,ite, jts,jte, kts,kte
1131 REAL, DIMENSION( its:ite , kts:kte) :: &
1132 qrs, &
1133 rslope, &
1134 rslopeb, &
1135 rslope2, &
1136 rslope3, &
1137 vt, &
1138 den, &
1139 denfac, &
1140 t
1141 REAL, PARAMETER :: t0c = 273.15
1142 REAL, DIMENSION( its:ite , kts:kte ) :: &
1143 n0sfac
1144 REAL :: lamdas, x, y, z, supcol
1145 integer :: i, j, k
1146 !----------------------------------------------------------------
1147 ! size distributions: (x=mixing ratio, y=air density):
1148 ! valid for mixing ratio > 1.e-9 kg/kg.
1149 lamdas(x,y,z)= sqrt(sqrt(pidn0s*z/(x*y))) ! (pidn0s*z/(x*y))**.25
1150 !
1151 do k = kts, kte
1152 do i = its, ite
1153 supcol = t0c-t(i,k)
1154 !---------------------------------------------------------------
1155 ! n0s: Intercept parameter for snow [m-4] [HDC 6]
1156 !---------------------------------------------------------------
1157 n0sfac(i,k) = max(min(exp(alpha*supcol),n0smax/n0s),1.)
1158 if(qrs(i,k).le.qcrmin)then
1159 rslope(i,k) = rslopesmax
1160 rslopeb(i,k) = rslopesbmax
1161 rslope2(i,k) = rslopes2max
1162 rslope3(i,k) = rslopes3max
1163 else
1164 rslope(i,k) = 1./lamdas(qrs(i,k),den(i,k),n0sfac(i,k))
1165 rslopeb(i,k) = rslope(i,k)**bvts
1166 rslope2(i,k) = rslope(i,k)*rslope(i,k)
1167 rslope3(i,k) = rslope2(i,k)*rslope(i,k)
1168 endif
1169 vt(i,k) = pvts*rslopeb(i,k)*denfac(i,k)
1170 if(qrs(i,k).le.0.0) vt(i,k) = 0.0
1171 enddo
1172 enddo
1173 END subroutine slope_snow
1174 ! SUBROUTINE nislfv_rain_plm(im0,km,den, denfac, tk, dz, ww0,qq0,precip0,dt,id,iter,irestrict,lon,lat,doit,call)
1175 SUBROUTINE nislfv_rain_plm(im, km,denl,denfacl,tkl,dzl,wwl,rql,precip ,dt,id,iter, &
1176 irestrict, lon, lat, doit, call )
1177 !-------------------------------------------------------------------
1178 !
1179 ! for non-iteration semi-Lagrangain forward advection for cloud
1180 ! with mass conservation and positive definite advection
1181 ! 2nd order interpolation with monotonic piecewise linear method
1182 ! this routine is under assumption of decfl < 1 for semi_Lagrangian
1183 !
1184 ! dzl depth of model layer in meter
1185 ! wwl terminal velocity at model layer m/s
1186 ! rql cloud density*mixing ration
1187 ! precip precipitation
1188 ! dt time step
1189 ! id kind of precip: 0 test case; 1 raindrop 2: snow
1190 ! iter how many time to guess mean terminal velocity: 0 pure forward.
1191 ! 0 : use departure wind for advection
1192 ! 1 : use mean wind for advection
1193 ! > 1 : use mean wind after iter-1 iterations
1194 !
1195 ! author: hann-ming henry juang <henry.juang@noaa.gov>
1196 ! implemented by song-you hong
1197 !
1198 implicit none
1199 integer irestrict,lon,lat,call
1200 logical doit
1201 integer im,km,id
1202 real dt
1203 real dzl(im,km),wwl(im,km),rql(im,km),precip(im)
1204 real denl(im,km),denfacl(im,km),tkl(im,km)
1205 !
1206 integer i,k,n,m,kk,kb,kt,iter
1207 real tl,tl2,qql,dql,qqd
1208 real th,th2,qqh,dqh
1209 real zsum,qsum,dim,dip,c1,con1,fa1,fa2
1210 real allold, allnew, zz, dzamin, cflmax, decfl
1211 real dz(km), ww(km), qq(km), wd(km), wa(km), was(km)
1212 real den(km), denfac(km), tk(km)
1213 real wi(km+1), zi(km+1), za(km+1)
1214 real qn(km), qr(km),tmp(km),tmp1(km),tmp2(km),tmp3(km)
1215 real dza(km+1), qa(km+1), qmi(km+1), qpi(km+1)
1216 !
1217 precip(:) = 0.0
1218 !
1219 i_loop : do i=1,irestrict !im
1220 ! -----------------------------------
1221 dz(:) = dzl(i,:)
1222 qq(:) = rql(i,:)
1223 ww(:) = wwl(i,:)
1224 den(:) = denl(i,:)
1225 denfac(:) = denfacl(i,:)
1226 tk(:) = tkl(i,:)
1227 ! skip for no precipitation for all layers
1228 allold = 0.0
1229 do k=1,km
1230 allold = allold + qq(k)
1231 enddo
1232 if(allold.le.0.0) then
1233 cycle i_loop
1234 endif
1235 !
1236 ! compute interface values
1237 zi(1)=0.0
1238 do k=1,km
1239 zi(k+1) = zi(k)+dz(k)
1240 enddo
1241 !
1242 ! save departure wind
1243 wd(:) = ww(:)
1244 n=1
1245 100 continue
1246 ! plm is 2nd order, we can use 2nd order wi or 3rd order wi
1247 ! 2nd order interpolation to get wi
1248 wi(1) = ww(1)
1249 wi(km+1) = ww(km)
1250 do k=2,km
1251 wi(k) = (ww(k)*dz(k-1)+ww(k-1)*dz(k))/(dz(k-1)+dz(k))
1252 enddo
1253 ! 3rd order interpolation to get wi
1254 fa1 = 9./16.
1255 fa2 = 1./16.
1256 wi(1) = ww(1)
1257 wi(2) = 0.5*(ww(2)+ww(1))
1258 do k=3,km-1
1259 wi(k) = fa1*(ww(k)+ww(k-1))-fa2*(ww(k+1)+ww(k-2))
1260 enddo
1261 wi(km) = 0.5*(ww(km)+ww(km-1))
1262 wi(km+1) = ww(km)
1263 !
1264 ! terminate of top of raingroup
1265 do k=2,km
1266 if( ww(k).eq.0.0 ) wi(k)=ww(k-1)
1267 enddo
1268 !
1269 ! diffusivity of wi
1270 con1 = 0.05
1271 do k=km,1,-1
1272 decfl = (wi(k+1)-wi(k))*dt/dz(k)
1273 if( decfl .gt. con1 ) then
1274 wi(k) = wi(k+1) - con1*dz(k)/dt
1275 endif
1276 enddo
1277 ! compute arrival point
1278 do k=1,km+1
1279 za(k) = zi(k) - wi(k)*dt
1280 enddo
1281 !
1282 do k=1,km
1283 dza(k) = za(k+1)-za(k)
1284 enddo
1285 dza(km+1) = zi(km+1) - za(km+1)
1286 !
1287 ! computer deformation at arrival point
1288 do k=1,km
1289 qa(k) = qq(k)*dz(k)/dza(k)
1290 qr(k) = qa(k)/den(k)
1291 enddo
1292 qa(km+1) = 0.0
1293 ! call maxmin(km,1,qa,' arrival points ')
1294 !
1295 ! compute arrival terminal velocity, and estimate mean terminal velocity
1296 ! then back to use mean terminal velocity
1297 if( n.le.iter ) then
1298 if (id.eq.1) then
1299 call slope_rain(qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,1,1,1,km)
1300 else
1301 call slope_snow(qr,den,denfac,tk,tmp,tmp1,tmp2,tmp3,wa,1,1,1,km)
1302 endif
1303 if( n.ge.2 ) wa(1:km)=0.5*(wa(1:km)+was(1:km))
1304 do k=1,km
1305 !#ifdef DEBUG
1306 ! print*,' slope_wsm3 ',qr(k)*1000.,den(k),denfac(k),tk(k),tmp(k),tmp1(k),tmp2(k),ww(k),wa(k)
1307 !#endif
1308 ! mean wind is average of departure and new arrival winds
1309 ww(k) = 0.5* ( wd(k)+wa(k) )
1310 enddo
1311 was(:) = wa(:)
1312 n=n+1
1313 go to 100
1314 endif
1315 !
1316 ! estimate values at arrival cell interface with monotone
1317 do k=2,km
1318 dip=(qa(k+1)-qa(k))/(dza(k+1)+dza(k))
1319 dim=(qa(k)-qa(k-1))/(dza(k-1)+dza(k))
1320 if( dip*dim.le.0.0 ) then
1321 qmi(k)=qa(k)
1322 qpi(k)=qa(k)
1323 else
1324 qpi(k)=qa(k)+0.5*(dip+dim)*dza(k)
1325 qmi(k)=2.0*qa(k)-qpi(k)
1326 if( qpi(k).lt.0.0 .or. qmi(k).lt.0.0 ) then
1327 qpi(k) = qa(k)
1328 qmi(k) = qa(k)
1329 endif
1330 endif
1331 enddo
1332 qpi(1)=qa(1)
1333 qmi(1)=qa(1)
1334 qmi(km+1)=qa(km+1)
1335 qpi(km+1)=qa(km+1)
1336 !
1337 ! interpolation to regular point
1338 qn = 0.0
1339 kb=1
1340 kt=1
1341 intp : do k=1,km
1342 kb=max(kb-1,1)
1343 kt=max(kt-1,1)
1344 ! find kb and kt
1345 if( zi(k).ge.za(km+1) ) then
1346 exit intp
1347 else
1348 find_kb : do kk=kb,km
1349 if( zi(k).le.za(kk+1) ) then
1350 kb = kk
1351 exit find_kb
1352 else
1353 cycle find_kb
1354 endif
1355 enddo find_kb
1356 find_kt : do kk=kt,km
1357 if( zi(k+1).le.za(kk) ) then
1358 kt = kk
1359 exit find_kt
1360 else
1361 cycle find_kt
1362 endif
1363 enddo find_kt
1364 kt = kt - 1
1365 ! compute q with piecewise constant method
1366 if( kt.eq.kb ) then
1367 tl=(zi(k)-za(kb))/dza(kb)
1368 th=(zi(k+1)-za(kb))/dza(kb)
1369 tl2=tl*tl
1370 th2=th*th
1371 qqd=0.5*(qpi(kb)-qmi(kb))
1372 qqh=qqd*th2+qmi(kb)*th
1373 qql=qqd*tl2+qmi(kb)*tl
1374 qn(k) = (qqh-qql)/(th-tl)
1375 else if( kt.gt.kb ) then
1376 tl=(zi(k)-za(kb))/dza(kb)
1377 tl2=tl*tl
1378 qqd=0.5*(qpi(kb)-qmi(kb))
1379 qql=qqd*tl2+qmi(kb)*tl
1380 dql = qa(kb)-qql
1381 zsum = (1.-tl)*dza(kb)
1382 qsum = dql*dza(kb)
1383 if( kt-kb.gt.1 ) then
1384 do m=kb+1,kt-1
1385 zsum = zsum + dza(m)
1386 qsum = qsum + qa(m) * dza(m)
1387 enddo
1388 endif
1389 th=(zi(k+1)-za(kt))/dza(kt)
1390 th2=th*th
1391 qqd=0.5*(qpi(kt)-qmi(kt))
1392 dqh=qqd*th2+qmi(kt)*th
1393 zsum = zsum + th*dza(kt)
1394 qsum = qsum + dqh*dza(kt)
1395 qn(k) = qsum/zsum
1396 endif
1397 cycle intp
1398 endif
1399 !
1400 enddo intp
1401 !
1402 ! rain out
1403 sum_precip: do k=1,km
1404 if( za(k).lt.0.0 .and. za(k+1).lt.0.0 ) then
1405 precip(i) = precip(i) + qa(k)*dza(k)
1406 cycle sum_precip
1407 else if ( za(k).lt.0.0 .and. za(k+1).ge.0.0 ) then
1408 precip(i) = precip(i) + qa(k)*(0.0-za(k))
1409 exit sum_precip
1410 endif
1411 exit sum_precip
1412 enddo sum_precip
1413 !
1414 ! replace the new values
1415 rql(i,:) = qn(:)
1416
1417 !
1418 ! ----------------------------------
1419 enddo i_loop
1420 !
1421 END SUBROUTINE nislfv_rain_plm
1422 #else
1423 SUBROUTINE slope_rain(qrs,den,denfac,t,rslope,rslopeb, &
1424 vt,irestrict,kts,kte,lmask)
1425 IMPLICIT NONE
1426 INTEGER :: irestrict,kts,kte
1427 REAL, DIMENSION( its:ite , kts:kte) :: &
1428 qrs, &
1429 vt, &
1430 den, &
1431 denfac, &
1432 t
1433 REAL, DIMENSION( its:ite ) :: &
1434 rslope, &
1435 rslopeb
1436 REAL, PARAMETER :: t0c = 273.15
1437 REAL :: lamdar, x, y, z
1438 LOGICAL*4 :: lmask(its:ite)
1439 integer :: i, j, k
1440 !----------------------------------------------------------------
1441 ! size distributions: (x=mixing ratio, y=air density):
1442 ! valid for mixing ratio > 1.e-9 kg/kg.
1443 lamdar(x,y)= sqrt(sqrt(pidn0r/(x*y))) ! (pidn0r/(x*y))**.25
1444 #define LAMDAR(x,y) sqrt(sqrt(pidn0r/((x)*(y))))
1445 !
1446 do k = kts, kte
1447 WHERE(lmask)
1448 WHERE(qrs(:,k).le.qcrmin)
1449 rslope(:) = rslopermax
1450 rslopeb(:) = rsloperbmax
1451 ELSEWHERE
1452 rslope(:) = 1./LAMDAR(qrs(:,k),den(:,k))
1453 rslopeb(:) = rslope(:)**bvtr
1454 ! rslopeb(:) = exp(log(rslope(:))*(bvtr))
1455 ENDWHERE
1456 WHERE(qrs(:,k).le.0.0)
1457 vt(:,k) = 0.0
1458 ELSEWHERE
1459 vt(:,k) = pvtr*rslopeb(:)*denfac(:,k)
1460 ENDWHERE
1461 ENDWHERE
1462 enddo
1463 #undef LAMDAR
1464 END SUBROUTINE slope_rain
1465 !------------------------------------------------------------------------------
1466 SUBROUTINE slope_snow(qrs,den,denfac,t,rslope,rslopeb, &
1467 vt,irestrict,kts,kte,lmask)
1468 IMPLICIT NONE
1469 INTEGER :: irestrict,kts,kte
1470 REAL, DIMENSION( its:ite , kts:kte) :: &
1471 qrs, &
1472 vt, &
1473 den, &
1474 denfac, &
1475 t
1476 REAL, DIMENSION( its:ite ) :: &
1477 rslope, &
1478 rslopeb
1479 REAL, PARAMETER :: t0c = 273.15
1480 integer :: i, j, k
1481 LOGICAL*4 lmask(CHUNK)
1482 !----------------------------------------------------------------
1483 ! size distributions: (x=mixing ratio, y=air density):
1484 ! valid for mixing ratio > 1.e-9 kg/kg.
1485 #define LAMDAS(x,y,z) sqrt(sqrt(pidn0s*(z)/((x)*(y))))
1486 !
1487 do k = kts, kte
1488 WHERE(lmask)
1489 !---------------------------------------------------------------
1490 ! n0s: Intercept parameter for snow [m-4] [HDC 6]
1491 !---------------------------------------------------------------
1492 WHERE(qrs(:,k).le.qcrmin)
1493 rslope(:) = rslopesmax
1494 rslopeb(:) = rslopesbmax
1495 ELSEWHERE
1496 rslope(:) = 1./LAMDAS(qrs(:,k),den(:,k),(max(min(exp(alpha*(t0c-t(:,k))),n0smax/n0s),1.)))
1497 rslopeb(:) = rslope(:)**bvts
1498 ! rslopeb(:) = exp(log(rslope(:))*(bvts))
1499 ENDWHERE
1500 WHERE(qrs(:,k).le.0.0)
1501 vt(:,k) = 0.0
1502 ELSEWHERE
1503 vt(:,k) = pvts*rslopeb(:)*denfac(:,k)
1504 ENDWHERE
1505 ENDWHERE
1506 enddo
1507 #undef LAMDAS
1508 END SUBROUTINE slope_snow
1509 !-------------------------------------------------------------------
1510 ! SUBROUTINE nislfv_rain_plm(im,km,denl,denfacl,tkl,dzl,wwl,rql,precip,dt,id,iter)
1511 SUBROUTINE nislfv_rain_plm(im0,km,den,denfac,tk,dz,ww0,qq0,precip0,dt,id,iter,irestrict,lon,lat,doit,call)
1512 !-------------------------------------------------------------------
1513 !
1514 ! for non-iteration semi-Lagrangain forward advection for cloud
1515 ! with mass conservation and positive definite advection
1516 ! 2nd order interpolation with monotonic piecewise linear method
1517 ! this routine is under assumption of decfl < 1 for semi_Lagrangian
1518 !
1519 ! dzl depth of model layer in meter
1520 ! wwl terminal velocity at model layer m/s
1521 ! rql cloud density*mixing ration
1522 ! precip precipitation
1523 ! dt time step
1524 ! id kind of precip: 0 test case; 1 raindrop 2: snow
1525 ! iter how many time to guess mean terminal velocity: 0 pure forward.
1526 ! 0 : use departure wind for advection
1527 ! 1 : use mean wind for advection
1528 ! > 1 : use mean wind after iter-1 iterations
1529 !
1530 ! author: hann-ming henry juang <henry.juang@noaa.gov>
1531 ! implemented by song-you hong
1532 !
1533 #define im nx
1534 !#define km nk
1535 implicit none
1536 integer,intent(in ) :: im0,km,id,irestrict
1537 real, intent(in ) :: den(im,km)
1538 real, intent(in ) :: denfac(im,km)
1539 real, intent(in ) :: tk(im,km)
1540 real, intent(in ) :: dz(im,km)
1541 real, intent(in ) :: ww0(im,km)
1542 real, intent( out) :: precip0(im)
1543 real, intent(inout) :: qq0(im,km)
1544 real, intent(in ) :: dt
1545 logical, intent(in) :: doit
1546 integer :: call
1547 !
1548 integer i,k,m,kk,iter
1549 integer n
1550 real dim(im),dip(im),c1,con1,fa1,fa2
1551 real allold(im), allnew(im), zz(im), dzamin(im), cflmax(im), decfl(im)
1552
1553 real qq(im,km), ww(im,km),precip(im)
1554 real wd(im,km), wa(im,km), was(im,km)
1555 real wi(im,km+1), zi(im,km+1), za(im,km+1)
1556 real qn(im,km),qr(im,km),tmp(im,km),tmp1(im,km),tmp2(im,km),tmp3(im,km)
1557 real dza(im,km+1), qa(im,km+1), qmi(im,km+1), qpi(im,km+1)
1558 logical*4 lmask(im)
1559
1560 INTEGER minkb, minkt
1561 LOGICAL, DIMENSION(CHUNK) :: intp_mask, tmask
1562 INTEGER, DIMENSION(CHUNK) :: kb, kt
1563 REAL, DIMENSION(CHUNK) :: tl,tl2,th,th2,qqd,qqh,qql,zsum,qsum,dql,dqh
1564 REAL, DIMENSION(CHUNK) :: za_gath_t,za_gath_b
1565 REAL, DIMENSION(CHUNK) :: qa_gath_b
1566 REAL, DIMENSION(CHUNK) :: dza_gath_t,dza_gath_b
1567 REAL, DIMENSION(CHUNK) :: qpi_gath_t,qpi_gath_b
1568 REAL, DIMENSION(CHUNK) :: qmi_gath_t,qmi_gath_b
1569
1570 integer, intent(in) :: lat,lon
1571 !
1572 precip = 0.0
1573
1574 ! -----------------------------------
1575 qq = qq0
1576 ww = ww0
1577 ! skip for no precipitation for all layers
1578 lmask = .false.
1579 do i = 1,min( irestrict, im )
1580 lmask(i) = .true.
1581 enddo
1582 allold = 0.0
1583 do k=1,km
1584 where(lmask)allold = allold + qq(:,k)
1585 enddo
1586 lmask = lmask .and. ( allold .gt. 0.0 )
1587 IF ( .NOT. ANY(lmask) ) THEN
1588 precip0 = precip
1589 RETURN
1590 ENDIF
1591
1592 !
1593 ! compute interface values
1594 zi(:,1)=0.0
1595 do k=1,km
1596 where(lmask) zi(:,k+1) = zi(:,k)+dz(:,k)
1597 enddo
1598 !
1599 ! save departure wind
1600 wd = ww
1601 DO n = 0, iter
1602 where(lmask)
1603 ! plm is 2nd order, we can use 2nd order wi or 3rd order wi
1604 ! 2nd order interpolation to get wi
1605 wi(:,1) = ww(:,1)
1606 wi(:,km+1) = ww(:,km)
1607 endwhere
1608 do k=2,km
1609 where(lmask)wi(:,k) = (ww(:,k)*dz(:,k-1)+ww(:,k-1)*dz(:,k))/(dz(:,k-1)+dz(:,k))
1610 enddo
1611
1612 ! 3rd order interpolation to get wi
1613 fa1 = 9./16.
1614 fa2 = 1./16.
1615 where(lmask)
1616 wi(:,1) = ww(:,1)
1617 wi(:,2) = 0.5*(ww(:,2)+ww(:,1))
1618 endwhere
1619 do k=3,km-1
1620 where(lmask)wi(:,k) = fa1*(ww(:,k)+ww(:,k-1))-fa2*(ww(:,k+1)+ww(:,k-2))
1621 enddo
1622 where(lmask)
1623 wi(:,km) = 0.5*(ww(:,km)+ww(:,km-1))
1624 wi(:,km+1) = ww(:,km)
1625 endwhere
1626 !
1627 ! terminate of top of raingroup
1628 do k=2,km
1629 where(lmask .and. ww(:,k).eq.0.0 ) wi(:,k)=ww(:,k-1)
1630 enddo
1631 !
1632 ! diffusivity of wi
1633 con1 = 0.05
1634 do k=km,1,-1
1635 where (lmask)
1636 decfl = (wi(:,k+1)-wi(:,k))*dt/dz(:,k)
1637 elsewhere
1638 decfl = 0.
1639 endwhere
1640 where (lmask )
1641 where (decfl .gt. con1 )
1642 wi(:,k) = wi(:,k+1) - con1*dz(:,k)/dt
1643 endwhere
1644 endwhere
1645 enddo
1646
1647 ! compute arrival point
1648 do k=1,km+1
1649 where (lmask) za(:,k) = zi(:,k) - wi(:,k)*dt
1650 enddo
1651 !
1652 do k=1,km
1653 where (lmask) dza(:,k) = za(:,k+1)-za(:,k)
1654 enddo
1655 where (lmask) dza(:,km+1) = zi(:,km+1) - za(:,km+1)
1656
1657 !
1658 ! compute deformation at arrival point
1659 do k=1,km
1660 where (lmask) qa(:,k) = qq(:,k)*dz(:,k)/dza(:,k)
1661 where (lmask) qr(:,k) = qa(:,k)/den(:,k)
1662 enddo
1663 where (lmask) qa(:,km+1) = 0.0
1664 ! compute arrival terminal velocity, and estimate mean terminal velocity
1665 ! then back to use mean terminal velocity
1666 if( n.le.iter-1 ) then
1667 if (id.eq.1) then
1668 call slope_rain(qr,den,denfac,tk,tmp,tmp1,wa,irestrict,1,km,lmask)
1669 else
1670 call slope_snow(qr,den,denfac,tk,tmp,tmp1,wa,irestrict,1,km,lmask)
1671 endif
1672 do k=1,km
1673 if( n.ge.1 ) where (lmask) wa(:,k)=0.5*(wa(:,k)+was(:,k))
1674 where (lmask ) ww(:,k) = 0.5* ( wd(:,k)+wa(:,k) )
1675 enddo
1676 was = wa
1677 endif
1678 ENDDO ! n loop
1679 !
1680 ! estimate values at arrival cell interface with monotone
1681 do k=2,km
1682 where (lmask )
1683 dip=(qa(:,k+1)-qa(:,k))/(dza(:,k+1)+dza(:,k))
1684 dim=(qa(:,k)-qa(:,k-1))/(dza(:,k-1)+dza(:,k))
1685 where( dip*dim.le.0.0 )
1686 qmi(:,k)=qa(:,k)
1687 qpi(:,k)=qa(:,k)
1688 elsewhere
1689 qpi(:,k)=qa(:,k)+0.5*(dip+dim)*dza(:,k)
1690 qmi(:,k)=2.0*qa(:,k)-qpi(:,k)
1691 where( qpi(:,k).lt.0.0 .or. qmi(:,k).lt.0.0 )
1692 qpi(:,k) = qa(:,k)
1693 qmi(:,k) = qa(:,k)
1694 endwhere
1695 endwhere
1696 endwhere
1697 enddo
1698 where (lmask )
1699 qpi(:,1)=qa(:,1)
1700 qmi(:,1)=qa(:,1)
1701 qmi(:,km+1)=qa(:,km+1)
1702 qpi(:,km+1)=qa(:,km+1)
1703 endwhere
1704 !
1705 ! interpolation to regular point
1706 qn = 0.0
1707 kb=1 ! kb is a vector
1708 kt=1 ! kt is a vector
1709 INTP : do k=1,km
1710 kb=max(kb-1,1)
1711 kt=max(kt-1,1)
1712 ! find kb and kt
1713 intp_mask = ( zi(:,k).lt.za(:,km+1) .AND. lmask )
1714 tmask = intp_mask
1715 minkb = 999
1716 minkt = 999
1717 DO i=1,CHUNK
1718 IF ( tmask(i) .AND. kb(i) .lt. minkb ) minkb = kb(i)
1719 IF ( tmask(i) .AND. kt(i) .lt. minkt ) minkt = kt(i)
1720 ENDDO
1721 find_kb : do kk=minkb,km
1722 WHERE ( tmask .AND. zi(:,k).le.za(:,kk+1) )
1723 kb = kk
1724 tmask = .FALSE.
1725 END WHERE
1726 enddo find_kb
1727
1728 tmask = intp_mask
1729 find_kt : do kk=minkt,km
1730 WHERE ( tmask .AND. zi(:,k+1).le.za(:,kk) )
1731 kt = kk
1732 tmask = .FALSE.
1733 END WHERE
1734 enddo find_kt
1735 kt = max(kt - 1,1)
1736
1737 !#define RANGE_CHECKING
1738 #ifndef RANGE_CHECKING
1739 # define DX1 (i+(kb(i)-1)*im),1
1740 # define DX2 (i+(kt(i)-1)*im),1
1741 #else
1742 # define DX1 i,kb(i)
1743 # define DX2 i,kt(i)
1744 #endif
1745 !DEC$ SIMD
1746 DO i = 1, CHUNK
1747 qa_gath_b(i) = qa(DX1)
1748 za_gath_b(i) = za(DX1)
1749 dza_gath_b(i) = dza(DX1)
1750 qpi_gath_b(i) = qpi(DX1)
1751 qmi_gath_b(i) = qmi(DX1)
1752 ENDDO
1753 !DEC$ SIMD
1754 DO i = 1, CHUNK
1755 za_gath_t(i) = za(DX2)
1756 dza_gath_t(i) = dza(DX2)
1757 qpi_gath_t(i) = qpi(DX2)
1758 qmi_gath_t(i) = qmi(DX2)
1759 ENDDO
1760
1761 WHERE ( kt .eq. kb .AND. intp_mask )
1762 tl=(zi(:,k)-za_gath_b)/dza_gath_b
1763 th=(zi(:,k+1)-za_gath_b)/dza_gath_b
1764 tl2 = tl*tl
1765 th2 = th*th
1766 qqd=0.5*(qpi_gath_b-qmi_gath_b)
1767 qqh=qqd*th2+qmi_gath_b*th
1768 qql=qqd*tl2+qmi_gath_b*tl
1769 qn(:,k) = (qqh-qql)/(th-tl)
1770 ELSE WHERE ( kt .gt. kb .AND. intp_mask )
1771 tl=(zi(:,k)-za_gath_b)/dza_gath_b
1772 tl2=tl*tl
1773 qqd=0.5*(qpi_gath_b-qmi_gath_b)
1774 qql=qqd*tl2+qmi_gath_b*tl
1775 dql = qa_gath_b-qql
1776 zsum = (1.-tl)*dza_gath_b
1777 qsum = dql*dza_gath_b
1778 END WHERE
1779 DO i = 1, CHUNK
1780 if( kt(i)-kb(i).gt.1 .AND. intp_mask(i) ) then
1781 do m=kb(i)+1,kt(i)-1
1782 zsum(i) = zsum(i) + dza(i,m)
1783 qsum(i) = qsum(i) + qa(i,m) * dza(i,m)
1784 enddo
1785 endif
1786 enddo
1787 WHERE ( kt .gt. kb .AND. intp_mask )
1788 th=(zi(:,k+1)-za_gath_t)/dza_gath_t
1789 th2 = th*th
1790 qqd=0.5*(qpi_gath_t-qmi_gath_t)
1791 dqh=qqd*th2+qmi_gath_t*th
1792 zsum = zsum + th*dza_gath_t
1793 qsum = qsum + dqh*dza_gath_t
1794 qn(:,k) = qsum/zsum
1795 END WHERE
1796 ENDDO intp
1797 !
1798 ! rain out
1799 intp_mask = lmask
1800 sum_precip: do k=1,km
1801 WHERE ( za(:,k).lt.0.0 .and. za(:,k+1).lt.0.0 .AND. intp_mask )
1802 precip = precip + qa(:,k)*dza(:,k)
1803 ELSE WHERE ( za(:,k).lt.0.0 .and. za(:,k+1).ge.0.0 .AND. intp_mask)
1804 precip = precip + qa(:,k)*(0.0-za(:,k))
1805 intp_mask = .FALSE.
1806 END WHERE
1807 enddo sum_precip
1808
1809 ! ENDDO ! i loop
1810
1811
1812 !
1813 ! replace the new values
1814 do k = 1,km
1815 where(lmask) qq0(:,k) = qn(:,k)
1816 enddo
1817 precip0 = 0.
1818 where (lmask) precip0 = precip
1819
1820 !
1821 ! ----------------------------------
1822 ! enddo i_loop
1823 !
1824
1825 END SUBROUTINE nislfv_rain_plm
1826
1827 #undef nx
1828 #undef nk
1829 #undef its
1830 #undef ite
1831 #undef ims
1832 #undef ime
1833 #undef kms
1834 #undef kme
1835 #undef im
1836 #undef km
1837
1838 #endif
1839
Mirror of the official repository at https://github.com/wrf-model/WRF