| blob | 46d18ffee3db8a1fe049d2fe399ead5485a8188c |
1 !WRF:MODEL_LAYER:PHYSICS
2 !
4 MODULE module_mp_lin
6 USE module_wrf_error
7 !
8 REAL , PARAMETER, PRIVATE :: RH = 1.0
9 ! REAL , PARAMETER, PRIVATE :: episp0 = 0.622*611.21
10 REAL , PARAMETER, PRIVATE :: xnor = 8.0e6
11 REAL , PARAMETER, PRIVATE :: xnos = 3.0e6
13 ! Lin
14 ! REAL , PARAMETER, PRIVATE :: xnog = 4.0e4
15 ! REAL , PARAMETER, PRIVATE :: rhograul = 917.
17 ! Hobbs
18 REAL , PARAMETER, PRIVATE :: xnog = 4.0e6
19 REAL , PARAMETER, PRIVATE :: rhograul = 400.
21 !
22 REAL , PARAMETER, PRIVATE :: &
23 qi0 = 1.0e-3, ql0 = 7.0e-4, qs0 = 6.0E-4, &
24 xmi50 = 4.8e-10, xmi40 = 2.46e-10, &
25 constb = 0.8, constd = 0.25, &
26 o6 = 1./6., cdrag = 0.6, &
27 avisc = 1.49628e-6, adiffwv = 8.7602e-5, &
28 axka = 1.4132e3, di50 = 1.0e-4, xmi = 4.19e-13, &
29 cw = 4.187e3, vf1s = 0.78, vf2s = 0.31, &
30 xni0 = 1.0e-2, xmnin = 1.05e-18, bni = 0.5, &
31 ci = 2.093e3
32 CONTAINS
34 !-------------------------------------------------------------------
35 ! Lin et al., 1983, JAM, 1065-1092, and
36 ! Rutledge and Hobbs, 1984, JAS, 2949-2972
37 !-------------------------------------------------------------------
38 SUBROUTINE lin_et_al(th &
39 ,qv, ql, qr &
40 ,qi, qs &
41 ,rho, pii, p &
42 ,dt_in &
43 ,z,ht, dz8w &
44 ,grav, cp, Rair, rvapor &
45 ,XLS, XLV, XLF, rhowater, rhosnow &
46 ,EP2,SVP1,SVP2,SVP3,SVPT0 &
47 , RAINNC, RAINNCV &
48 ,ids,ide, jds,jde, kds,kde &
49 ,ims,ime, jms,jme, kms,kme &
50 ,its,ite, jts,jte, kts,kte &
51 ! Optional
52 ,qlsink, precr, preci, precs, precg &
53 , F_QG,F_QNDROP &
54 , qg, qndrop &
55 )
56 !-------------------------------------------------------------------
57 IMPLICIT NONE
58 !-------------------------------------------------------------------
59 !
60 ! Shuhua 12/17/99
61 !
62 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
63 ims,ime, jms,jme, kms,kme , &
64 its,ite, jts,jte, kts,kte
66 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
67 INTENT(INOUT) :: &
68 th, &
69 qv, &
70 ql, &
71 qr
73 !
74 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
75 INTENT(IN ) :: &
76 rho, &
77 pii, &
78 p, &
79 dz8w
81 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
82 INTENT(IN ) :: z
86 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: ht
88 REAL, INTENT(IN ) :: dt_in, &
89 grav, &
90 Rair, &
91 rvapor, &
92 cp, &
93 XLS, &
94 XLV, &
95 XLF, &
96 rhowater, &
97 rhosnow
99 REAL, INTENT(IN ) :: EP2,SVP1,SVP2,SVP3,SVPT0
101 REAL, DIMENSION( ims:ime , jms:jme ), &
102 INTENT(INOUT) :: RAINNC, &
103 RAINNCV
105 ! Optional
107 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
108 OPTIONAL, &
109 INTENT(INOUT) :: &
110 qi, &
111 qs, &
112 qg, &
113 qndrop
115 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
116 OPTIONAL, INTENT(OUT ) :: &
117 qlsink, & ! cloud water conversion to rain (/s)
118 precr, & ! rain precipitation rate at all levels (kg/m2/s)
119 preci, & ! ice precipitation rate at all levels (kg/m2/s)
120 precs, & ! snow precipitation rate at all levels (kg/m2/s)
121 precg ! graupel precipitation rate at all levels (kg/m2/s)
123 LOGICAL, INTENT(IN), OPTIONAL :: F_QG, F_QNDROP
125 ! LOCAL VAR
127 INTEGER :: min_q, max_q
129 REAL, DIMENSION( its:ite , jts:jte ) &
130 :: rain, snow, graupel,ice
132 REAL, DIMENSION( kts:kte ) :: qvz, qlz, qrz, &
133 qiz, qsz, qgz, &
134 thz, &
135 tothz, rhoz, &
136 orhoz, sqrhoz, &
137 prez, zz, &
138 precrz, preciz, precsz, precgz, &
139 qndropz, &
140 dzw, preclw
142 LOGICAL :: flag_qg, flag_qndrop
143 !
144 REAL :: dt, pptrain, pptsnow, pptgraul, rhoe_s, &
145 gindex, pptice
146 real :: qndropconst
148 INTEGER :: i,j,k
149 !
150 flag_qg = .false.
151 flag_qndrop = .false.
152 IF ( PRESENT ( f_qg ) ) flag_qg = f_qg
153 IF ( PRESENT ( f_qndrop ) ) flag_qndrop = f_qndrop
154 !
155 dt=dt_in
156 rhoe_s=1.29
157 qndropconst=100.e6 !sg
158 gindex=1.0
160 IF (.not.flag_qg) gindex=0.
162 j_loop: DO j = jts, jte
163 i_loop: DO i = its, ite
164 !
165 !- write data from 3-D to 1-D
166 !
167 DO k = kts, kte
168 qvz(k)=qv(i,k,j)
169 qlz(k)=ql(i,k,j)
170 qrz(k)=qr(i,k,j)
171 thz(k)=th(i,k,j)
172 !
173 rhoz(k)=rho(i,k,j)
174 orhoz(k)=1./rhoz(k)
175 prez(k)=p(i,k,j)
176 sqrhoz(k)=sqrt(rhoe_s*orhoz(k))
177 tothz(k)=pii(i,k,j)
178 zz(k)=z(i,k,j)
179 dzw(k)=dz8w(i,k,j)
180 END DO
182 IF (flag_qndrop .AND. PRESENT( qndrop )) THEN
183 DO k = kts, kte
184 qndropz(k)=qndrop(i,k,j)
185 ENDDO
186 ELSE
187 DO k = kts, kte
188 qndropz(k)=qndropconst
189 ENDDO
190 ENDIF
192 DO k = kts, kte
193 qiz(k)=qi(i,k,j)
194 qsz(k)=qs(i,k,j)
195 ENDDO
197 IF ( flag_qg .AND. PRESENT( qg ) ) THEN
198 DO k = kts, kte
199 qgz(k)=qg(i,k,j)
200 ENDDO
201 ELSE
202 DO k = kts, kte
203 qgz(k)=0.
204 ENDDO
205 ENDIF
206 !
207 pptrain=0.
208 pptsnow=0.
209 pptgraul=0.
210 pptice=0.
211 CALL clphy1d( dt, qvz, qlz, qrz, qiz, qsz, qgz, &
212 qndropz,flag_qndrop, &
213 thz, tothz, rhoz, orhoz, sqrhoz, &
214 prez, zz, dzw, ht(I,J), preclw, &
215 precrz, preciz, precsz, precgz, &
216 pptrain, pptsnow, pptgraul, pptice, &
217 grav, cp, Rair, rvapor, gindex, &
218 XLS, XLV, XLF, rhowater, rhosnow, &
219 EP2,SVP1,SVP2,SVP3,SVPT0, &
220 kts, kte, i, j )
222 !
223 ! Precipitation from cloud microphysics -- only for one time step
224 !
225 ! unit is transferred from m to mm
227 !
228 rain(i,j)=pptrain
229 snow(i,j)=pptsnow
230 graupel(i,j)=pptgraul
231 ice(i,j)=pptice
232 !
233 RAINNCV(i,j)= pptrain + pptsnow + pptgraul + pptice
234 RAINNC(i,j)=RAINNC(i,j) + pptrain + pptsnow + pptgraul + pptice
236 !
237 !- update data from 1-D back to 3-D
238 !
239 !
240 IF ( present(qlsink) .and. present(precr) ) THEN !sg beg
241 DO k = kts, kte
242 if(ql(i,k,j)>1.e-20) then
243 qlsink(i,k,j)=-preclw(k)/ql(i,k,j)
244 else
245 qlsink(i,k,j)=0.
246 endif
247 precr(i,k,j)=precrz(k)
248 END DO
249 END IF !sg end
251 DO k = kts, kte
252 qv(i,k,j)=qvz(k)
253 ql(i,k,j)=qlz(k)
254 qr(i,k,j)=qrz(k)
255 th(i,k,j)=thz(k)
256 END DO
257 !
258 IF ( flag_qndrop .AND. PRESENT( qndrop ) ) THEN !sg beg
259 DO k = kts, kte
260 qndrop(i,k,j)=qndropz(k)
261 ENDDO
262 END IF !sg end
264 DO k = kts, kte
265 qi(i,k,j)=qiz(k)
266 qs(i,k,j)=qsz(k)
267 ENDDO
269 IF ( present(preci) ) THEN !sg beg
270 DO k = kts, kte
271 preci(i,k,j)=preciz(k)
272 ENDDO
273 END IF
275 IF ( present(precs) ) THEN
276 DO k = kts, kte
277 precs(i,k,j)=precsz(k)
278 ENDDO
279 END IF !sg end
281 IF ( flag_qg .AND. PRESENT( qg ) ) THEN
282 DO k = kts, kte
283 qg(i,k,j)=qgz(k)
284 ENDDO
286 IF ( present(precg) ) THEN !sg beg
287 DO k = kts, kte
288 precg(i,k,j)=precgz(k)
289 ENDDO !sg end
290 END IF
291 ELSE !sg beg
292 IF ( present(precg) ) precg(i,:,j)=0. !sg end
293 ENDIF
294 !
295 ENDDO i_loop
296 ENDDO j_loop
298 END SUBROUTINE lin_et_al
301 !-----------------------------------------------------------------------
302 SUBROUTINE clphy1d(dt, qvz, qlz, qrz, qiz, qsz, qgz, &
303 qndropz,flag_qndrop, &
304 thz, tothz, rho, orho, sqrho, &
305 prez, zz, dzw, zsfc, preclw, &
306 precrz, preciz, precsz, precgz, &
307 pptrain, pptsnow, pptgraul, &
308 pptice, grav, cp, Rair, rvapor, gindex, &
309 XLS, XLV, XLF, rhowater, rhosnow, &
310 EP2,SVP1,SVP2,SVP3,SVPT0, &
311 kts, kte, i, j )
312 !-----------------------------------------------------------------------
313 IMPLICIT NONE
314 !-----------------------------------------------------------------------
315 ! This program handles the vertical 1-D cloud micphysics
316 !-----------------------------------------------------------------------
317 ! avisc: constant in empirical formular for dynamic viscosity of air
318 ! =1.49628e-6 [kg/m/s] = 1.49628e-5 [g/cm/s]
319 ! adiffwv: constant in empirical formular for diffusivity of water
320 ! vapor in air
321 ! = 8.7602e-5 [kgm/s3] = 8.7602 [gcm/s3]
322 ! axka: constant in empirical formular for thermal conductivity of air
323 ! = 1.4132e3 [m2/s2/K] = 1.4132e7 [cm2/s2/K]
324 ! qi0: mixing ratio threshold for cloud ice aggregation [kg/kg]
325 ! xmi50: mass of a 50 micron ice crystal
326 ! = 4.8e-10 [kg] =4.8e-7 [g]
327 ! xmi40: mass of a 40 micron ice crystal
328 ! = 2.46e-10 [kg] = 2.46e-7 [g]
329 ! di50: diameter of a 50 micro (radius) ice crystal
330 ! =1.0e-4 [m]
331 ! xmi: mass of one cloud ice crystal
332 ! =4.19e-13 [kg] = 4.19e-10 [g]
333 ! oxmi=1.0/xmi
334 !
335 ! xni0=1.0e-2 [m-3] The value given in Lin et al. is wrong.(see
336 ! Hsie et al.(1980) and Rutledge and Hobbs(1983) )
337 ! bni=0.5 [K-1]
338 ! xmnin: mass of a natural ice nucleus
339 ! = 1.05e-18 [kg] = 1.05e-15 [g] This values is suggested by
340 ! Hsie et al. (1980)
341 ! = 1.0e-12 [kg] suggested by Rutlegde and Hobbs (1983)
342 ! rhowater: density of water=1.0 g/cm3=1000.0 kg/m3
343 ! consta: constant in empirical formular for terminal
344 ! velocity of raindrop
345 ! =2115.0 [cm**(1-b)/s] = 2115.0*0.01**(1-b) [m**(1-b)/s]
346 ! constb: constant in empirical formular for terminal
347 ! velocity of raindrop
348 ! =0.8
349 ! xnor: intercept parameter of the raindrop size distribution
350 ! = 0.08 cm-4 = 8.0e6 m-4
351 ! araut: time sacle for autoconversion of cloud water to raindrops
352 ! =1.0e-3 [s-1]
353 ! ql0: mixing ratio threshold for cloud watercoalescence [kg/kg]
354 ! vf1r: ventilation factors for rain =0.78
355 ! vf2r: ventilation factors for rain =0.31
356 ! rhosnow: density of snow=0.1 g/cm3=100.0 kg/m3
357 ! constc: constant in empirical formular for terminal
358 ! velocity of snow
359 ! =152.93 [cm**(1-d)/s] = 152.93*0.01**(1-d) [m**(1-d)/s]
360 ! constd: constant in empirical formular for terminal
361 ! velocity of snow
362 ! =0.25
363 ! xnos: intercept parameter of the snowflake size distribution
364 ! vf1s: ventilation factors for snow =0.78
365 ! vf2s: ventilation factors for snow =0.31
366 !
367 !----------------------------------------------------------------------
369 INTEGER, INTENT(IN ) :: kts, kte, i, j
371 REAL, DIMENSION( kts:kte ), &
372 INTENT(INOUT) :: qvz, qlz, qrz, qiz, qsz, &
373 qndropz, &
374 qgz, thz
376 REAL, DIMENSION( kts:kte ), &
377 INTENT(IN ) :: tothz, rho, orho, sqrho, &
378 prez, zz, dzw
380 REAL, INTENT(IN ) :: dt, grav, cp, Rair, rvapor, &
381 XLS, XLV, XLF, rhowater, &
382 rhosnow,EP2,SVP1,SVP2,SVP3,SVPT0
384 REAL, DIMENSION( kts:kte ), INTENT(OUT) :: preclw, &
385 precrz, preciz, precsz, precgz
387 REAL, INTENT(INOUT) :: pptrain, pptsnow, pptgraul, pptice
389 REAL, INTENT(IN ) :: zsfc
390 logical, intent(in) :: flag_qndrop !sg
392 ! local vars
394 REAL :: obp4, bp3, bp5, bp6, odp4, &
395 dp3, dp5, dp5o2
398 ! temperary vars
400 REAL :: tmp, tmp0, tmp1, tmp2,tmp3, &
401 tmp4,delta2,delta3, delta4, &
402 tmpa,tmpb,tmpc,tmpd,alpha1, &
403 qic, abi,abr, abg, odtberg, &
404 vti50,eiw,eri,esi,esr, esw, &
405 erw,delrs,term0,term1,araut, &
406 constg2, vf1r, vf2r,alpha2, &
407 Ap, Bp, egw, egi, egs, egr, &
408 constg, gdelta4, g1sdelt4, &
409 factor, tmp_r, tmp_s,tmp_g, &
410 qlpqi, rsat, a1, a2, xnin
412 INTEGER :: k
413 !
414 REAL, DIMENSION( kts:kte ) :: oprez, tem, temcc, theiz, qswz, &
415 qsiz, qvoqswz, qvoqsiz, qvzodt, &
416 qlzodt, qizodt, qszodt, qrzodt, &
417 qgzodt
419 REAL, DIMENSION( kts:kte ) :: psnow, psaut, psfw, psfi, praci, &
420 piacr, psaci, psacw, psdep, pssub, &
421 pracs, psacr, psmlt, psmltevp, &
422 prain, praut, pracw, prevp, pvapor, &
423 pclw, pladj, pcli, pimlt, pihom, &
424 pidw, piadj, pgraupel, pgaut, &
425 pgfr, pgacw, pgaci, pgacr, pgacs, &
426 pgacip,pgacrp,pgacsp,pgwet, pdry, &
427 pgsub, pgdep, pgmlt, pgmltevp, &
428 qschg, qgchg
429 !
431 REAL, DIMENSION( kts:kte ) :: qvsbar, rs0, viscmu, visc, diffwv, &
432 schmidt, xka
434 REAL, DIMENSION( kts:kte ) :: vtr, vts, vtg, &
435 vtrold, vtsold, vtgold, vtiold, &
436 xlambdar, xlambdas, xlambdag, &
437 olambdar, olambdas, olambdag
439 REAL :: episp0k, dtb, odtb, pi, pio4, &
440 pio6, oxLf, xLvocp, xLfocp, consta, &
441 constc, ocdrag, gambp4, gamdp4, &
442 gam4pt5, Cpor, oxmi, gambp3, gamdp3,&
443 gambp6, gam3pt5, gam2pt75, gambp5o2,&
444 gamdp5o2, cwoxlf, ocp, xni50, es
445 !
446 REAL :: qvmin=1.e-20
447 REAL :: gindex
448 REAL :: temc1,save1,save2,xni50mx
450 ! for terminal velocity flux
452 INTEGER :: min_q, max_q
453 REAL :: t_del_tv, del_tv, flux, fluxin, fluxout ,tmpqrz
454 LOGICAL :: notlast
455 !
457 !sg: begin
458 ! liqconc = liquid water content in gcm^-3
459 ! capn = droplet number concentration cm^-3
460 ! dis = relative dispersion (dimensionless) between 0.2 and 1.
461 ! Written by Yangang Liu based on Liu et al., GRL 32, 2005.
462 ! Autoconversion rate P = P0*T
463 ! p0 = rate function
464 ! kappa = constant in Long kernel
465 ! beta = Condensation rate constant
466 ! xc = Normalized critical mass
467 ! ***********************************************************
468 real liqconc, dis, beta, kappa, p0, xc, capn,rhocgs
469 if(flag_qndrop)then
470 dis = 0.5 ! droplet dispersion, set to 0.5 per SG 8-Nov-2006
471 ! Give empirical constants
472 kappa=1.1d10
473 ! Calculate Condensation rate constant
474 beta = (1.0d0+3.0d0*dis**2)*(1.0d0+4.0d0*dis**2)* &
475 (1.0d0+5.0d0*dis**2)/((1.0d0+dis**2)*(1.0d0+2.0d0*dis**2))
476 endif
477 !sg: end
479 dtb=dt
480 odtb=1./dtb
481 pi=acos(-1.)
482 pio4=acos(-1.)/4.
483 pio6=acos(-1.)/6.
484 ocp=1./cp
485 oxLf=1./xLf
486 xLvocp=xLv/cp
487 xLfocp=xLf/cp
488 consta=2115.0*0.01**(1-constb)
489 constc=152.93*0.01**(1-constd)
490 ocdrag=1./Cdrag
491 ! episp0k=RH*episp0
492 episp0k=RH*ep2*1000.*svp1
493 !
494 gambp4=ggamma(constb+4.)
495 gamdp4=ggamma(constd+4.)
496 gam4pt5=ggamma(4.5)
497 Cpor=cp/Rair
498 oxmi=1.0/xmi
499 gambp3=ggamma(constb+3.)
500 gamdp3=ggamma(constd+3.)
501 gambp6=ggamma(constb+6)
502 gam3pt5=ggamma(3.5)
503 gam2pt75=ggamma(2.75)
504 gambp5o2=ggamma((constb+5.)/2.)
505 gamdp5o2=ggamma((constd+5.)/2.)
506 cwoxlf=cw/xlf
507 delta2=0.
508 delta3=0.
509 delta4=0.
510 !
511 !-----------------------------------------------------------------------
512 ! oprez 1./prez ( prez : pressure)
513 ! qsw saturated mixing ratio on water surface
514 ! qsi saturated mixing ratio on ice surface
515 ! episp0k RH*e*saturated pressure at 273.15 K
516 ! qvoqsw qv/qsw
517 ! qvoqsi qv/qsi
518 ! qvzodt qv/dt
519 ! qlzodt ql/dt
520 ! qizodt qi/dt
521 ! qszodt qs/dt
522 ! qrzodt qr/dt
523 ! qgzodt qg/dt
524 !
525 ! temcc temperature in dregee C
526 !
528 obp4=1.0/(constb+4.0)
529 bp3=constb+3.0
530 bp5=constb+5.0
531 bp6=constb+6.0
532 odp4=1.0/(constd+4.0)
533 dp3=constd+3.0
534 dp5=constd+5.0
535 dp5o2=0.5*(constd+5.0)
536 !
537 do k=kts,kte
538 oprez(k)=1./prez(k)
539 enddo
541 do k=kts,kte
542 qlz(k)=amax1( 0.0,qlz(k) )
543 qiz(k)=amax1( 0.0,qiz(k) )
544 qvz(k)=amax1( qvmin,qvz(k) )
545 qsz(k)=amax1( 0.0,qsz(k) )
546 qrz(k)=amax1( 0.0,qrz(k) )
547 qgz(k)=amax1( 0.0,qgz(k) )
548 qndropz(k)=amax1( 0.0,qndropz(k) ) !sg
549 !
550 tem(k)=thz(k)*tothz(k)
551 temcc(k)=tem(k)-273.15
552 !
553 ! qswz(k)=episp0k*oprez(k)* &
554 ! exp( svp2*temcc(k)/(tem(k)-svp3) )
555 es=1000.*svp1*exp( svp2*temcc(k)/(tem(k)-svp3) )
556 qswz(k)=ep2*es/(prez(k)-es)
557 if (tem(k) .lt. 273.15 ) then
558 ! qsiz(k)=episp0k*oprez(k)* &
559 ! exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
560 es=1000.*svp1*exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
561 qsiz(k)=ep2*es/(prez(k)-es)
562 if (temcc(k) .lt. -40.0) qswz(k)=qsiz(k)
563 else
564 qsiz(k)=qswz(k)
565 endif
566 !
567 qvoqswz(k)=qvz(k)/qswz(k)
568 qvoqsiz(k)=qvz(k)/qsiz(k)
569 qvzodt(k)=amax1( 0.0,odtb*qvz(k) )
570 qlzodt(k)=amax1( 0.0,odtb*qlz(k) )
571 qizodt(k)=amax1( 0.0,odtb*qiz(k) )
572 qszodt(k)=amax1( 0.0,odtb*qsz(k) )
573 qrzodt(k)=amax1( 0.0,odtb*qrz(k) )
574 qgzodt(k)=amax1( 0.0,odtb*qgz(k) )
576 theiz(k)=thz(k)+(xlvocp*qvz(k)-xlfocp*qiz(k))/tothz(k)
577 enddo
580 !
581 !
582 !-----------------------------------------------------------------------
583 ! In this simple stable cloud parameterization scheme, only five
584 ! forms of water substance (water vapor, cloud water, cloud ice,
585 ! rain and snow are considered. The prognostic variables are total
586 ! water (qp),cloud water (ql), and cloud ice (qi). Rain and snow are
587 ! diagnosed following Nagata and Ogura, 1991, MWR, 1309-1337. Eq (A7).
588 ! the micro physics are based on (1) Hsie et al.,1980, JAM, 950-977 ;
589 ! (2) Lin et al., 1983, JAM, 1065-1092 ; (3) Rutledge and Hobbs, 1983,
590 ! JAS, 1185-1206 ; (4) Rutledge and Hobbs, 1984, JAS, 2949-2972.
591 !-----------------------------------------------------------------------
592 !
593 ! rhowater: density of water=1.0 g/cm3=1000.0 kg/m3
594 ! rhosnow: density of snow=0.1 g/cm3=100.0 kg/m3
595 ! xnor: intercept parameter of the raindrop size distribution
596 ! = 0.08 cm-4 = 8.0e6 m-4
597 ! xnos: intercept parameter of the snowflake size distribution
598 ! = 0.03 cm-4 = 3.0e6 m-4
599 ! xnog: intercept parameter of the graupel size distribution
600 ! = 4.0e-4 cm-4 = 4.0e4 m-4
601 ! consta: constant in empirical formular for terminal
602 ! velocity of raindrop
603 ! =2115.0 [cm**(1-b)/s] = 2115.0*0.01**(1-b) [m**(1-b)/s]
604 ! constb: constant in empirical formular for terminal
605 ! velocity of raindrop
606 ! =0.8
607 ! constc: constant in empirical formular for terminal
608 ! velocity of snow
609 ! =152.93 [cm**(1-d)/s] = 152.93*0.01**(1-d) [m**(1-d)/s]
610 ! constd: constant in empirical formular for terminal
611 ! velocity of snow
612 ! =0.25
613 ! avisc: constant in empirical formular for dynamic viscosity of air
614 ! =1.49628e-6 [kg/m/s] = 1.49628e-5 [g/cm/s]
615 ! adiffwv: constant in empirical formular for diffusivity of water
616 ! vapor in air
617 ! = 8.7602e-5 [kgm/s3] = 8.7602 [gcm/s3]
618 ! axka: constant in empirical formular for thermal conductivity of air
619 ! = 1.4132e3 [m2/s2/K] = 1.4132e7 [cm2/s2/K]
620 ! qi0: mixing ratio threshold for cloud ice aggregation [kg/kg]
621 ! = 1.0e-3 g/g = 1.0e-3 kg/gk
622 ! ql0: mixing ratio threshold for cloud watercoalescence [kg/kg]
623 ! = 2.0e-3 g/g = 2.0e-3 kg/gk
624 ! qs0: mixing ratio threshold for snow aggregation
625 ! = 6.0e-4 g/g = 6.0e-4 kg/gk
626 ! xmi50: mass of a 50 micron ice crystal
627 ! = 4.8e-10 [kg] =4.8e-7 [g]
628 ! xmi40: mass of a 40 micron ice crystal
629 ! = 2.46e-10 [kg] = 2.46e-7 [g]
630 ! di50: diameter of a 50 micro (radius) ice crystal
631 ! =1.0e-4 [m]
632 ! xmi: mass of one cloud ice crystal
633 ! =4.19e-13 [kg] = 4.19e-10 [g]
634 ! oxmi=1.0/xmi
635 !
638 ! if gindex=1.0 include graupel
639 ! if gindex=0. no graupel
640 !
641 !
642 do k=kts,kte
643 psnow(k)=0.0
644 psaut(k)=0.0
645 psfw(k)=0.0
646 psfi(k)=0.0
647 praci(k)=0.0
648 piacr(k)=0.0
649 psaci(k)=0.0
650 psacw(k)=0.0
651 psdep(k)=0.0
652 pssub(k)=0.0
653 pracs(k)=0.0
654 psacr(k)=0.0
655 psmlt(k)=0.0
656 psmltevp(k)=0.0
657 !
658 prain(k)=0.0
659 praut(k)=0.0
660 pracw(k)=0.0
661 prevp(k)=0.0
662 !
663 pvapor(k)=0.0
664 !
665 pclw(k)=0.0
666 preclw(k)=0.0 !sg
667 pladj(k)=0.0
668 !
669 pcli(k)=0.0
670 pimlt(k)=0.0
671 pihom(k)=0.0
672 pidw(k)=0.0
673 piadj(k)=0.0
674 enddo
676 !
677 !!! graupel
678 !
679 do k=kts,kte
680 pgraupel(k)=0.0
681 pgaut(k)=0.0
682 pgfr(k)=0.0
683 pgacw(k)=0.0
684 pgaci(k)=0.0
685 pgacr(k)=0.0
686 pgacs(k)=0.0
687 pgacip(k)=0.0
688 pgacrP(k)=0.0
689 pgacsp(k)=0.0
690 pgwet(k)=0.0
691 pdry(k)=0.0
692 pgsub(k)=0.0
693 pgdep(k)=0.0
694 pgmlt(k)=0.0
695 pgmltevp(k)=0.0
696 qschg(k)=0.
697 qgchg(k)=0.
698 enddo
699 !
700 !
701 ! Set rs0=episp0*oprez(k)
702 ! episp0=e*saturated pressure at 273.15 K
703 ! e = 0.622
704 !
705 DO k=kts,kte
706 rs0(k)=ep2*1000.*svp1/(prez(k)-1000.*svp1)
707 END DO
708 !
709 !***********************************************************************
710 ! Calculate precipitation fluxes due to terminal velocities.
711 !***********************************************************************
712 !
713 !- Calculate termianl velocity (vt?) of precipitation q?z
714 !- Find maximum vt? to determine the small delta t
715 !
716 !-- rain
717 !
718 t_del_tv=0.
719 del_tv=dtb
720 notlast=.true.
721 DO while (notlast)
722 !
723 min_q=kte
724 max_q=kts-1
725 !
726 do k=kts,kte-1
727 if (qrz(k) .gt. 1.0e-8) then
728 min_q=min0(min_q,k)
729 max_q=max0(max_q,k)
730 tmp1=sqrt(pi*rhowater*xnor/rho(k)/qrz(k))
731 tmp1=sqrt(tmp1)
732 vtrold(k)=o6*consta*gambp4*sqrho(k)/tmp1**constb
733 if (k .eq. 1) then
734 del_tv=amin1(del_tv,0.9*(zz(k)-zsfc)/vtrold(k))
735 else
736 del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtrold(k))
737 endif
738 else
739 vtrold(k)=0.
740 endif
741 enddo
743 if (max_q .ge. min_q) then
744 !
745 !- Check if the summation of the small delta t >= big delta t
746 ! (t_del_tv) (del_tv) (dtb)
748 t_del_tv=t_del_tv+del_tv
749 !
750 if ( t_del_tv .ge. dtb ) then
751 notlast=.false.
752 del_tv=dtb+del_tv-t_del_tv
753 endif
754 !
755 fluxin=0.
756 do k=max_q,min_q,-1
757 fluxout=rho(k)*vtrold(k)*qrz(k)
758 flux=(fluxin-fluxout)/rho(k)/dzw(k)
759 tmpqrz=qrz(k)
760 qrz(k)=qrz(k)+del_tv*flux
761 fluxin=fluxout
762 enddo
763 if (min_q .eq. 1) then
764 pptrain=pptrain+fluxin*del_tv
765 else
766 qrz(min_q-1)=qrz(min_q-1)+del_tv* &
767 fluxin/rho(min_q-1)/dzw(min_q-1)
768 endif
769 !
770 else
771 notlast=.false.
772 endif
773 ENDDO
775 !
776 !-- snow
777 !
778 t_del_tv=0.
779 del_tv=dtb
780 notlast=.true.
782 DO while (notlast)
783 !
784 min_q=kte
785 max_q=kts-1
786 !
787 do k=kts,kte-1
788 if (qsz(k) .gt. 1.0e-8) then
789 min_q=min0(min_q,k)
790 max_q=max0(max_q,k)
791 tmp1=sqrt(pi*rhosnow*xnos/rho(k)/qsz(k))
792 tmp1=sqrt(tmp1)
793 vtsold(k)=o6*constc*gamdp4*sqrho(k)/tmp1**constd
794 if (k .eq. 1) then
795 del_tv=amin1(del_tv,0.9*(zz(k)-zsfc)/vtsold(k))
796 else
797 del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtsold(k))
798 endif
799 else
800 vtsold(k)=0.
801 endif
802 enddo
804 if (max_q .ge. min_q) then
805 !
806 !
807 !- Check if the summation of the small delta t >= big delta t
808 ! (t_del_tv) (del_tv) (dtb)
810 t_del_tv=t_del_tv+del_tv
812 if ( t_del_tv .ge. dtb ) then
813 notlast=.false.
814 del_tv=dtb+del_tv-t_del_tv
815 endif
816 !
817 fluxin=0.
818 do k=max_q,min_q,-1
819 fluxout=rho(k)*vtsold(k)*qsz(k)
820 flux=(fluxin-fluxout)/rho(k)/dzw(k)
821 qsz(k)=qsz(k)+del_tv*flux
822 qsz(k)=amax1(0.,qsz(k))
823 fluxin=fluxout
824 enddo
825 if (min_q .eq. 1) then
826 pptsnow=pptsnow+fluxin*del_tv
827 else
828 qsz(min_q-1)=qsz(min_q-1)+del_tv* &
829 fluxin/rho(min_q-1)/dzw(min_q-1)
830 endif
831 !
832 else
833 notlast=.false.
834 endif
836 ENDDO
837 !
838 !-- grauupel
839 !
840 t_del_tv=0.
841 del_tv=dtb
842 notlast=.true.
843 !
844 DO while (notlast)
845 !
846 min_q=kte
847 max_q=kts-1
848 !
849 do k=kts,kte-1
850 if (qgz(k) .gt. 1.0e-8) then
851 min_q=min0(min_q,k)
852 max_q=max0(max_q,k)
853 tmp1=sqrt(pi*rhograul*xnog/rho(k)/qgz(k))
854 tmp1=sqrt(tmp1)
855 term0=sqrt(4.*grav*rhograul*0.33334/rho(k)/cdrag)
856 vtgold(k)=o6*gam4pt5*term0*sqrt(1./tmp1)
857 if (k .eq. 1) then
858 del_tv=amin1(del_tv,0.9*(zz(k)-zsfc)/vtgold(k))
859 else
860 del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtgold(k))
861 endif
862 else
863 vtgold(k)=0.
864 endif
865 enddo
867 if (max_q .ge. min_q) then
868 !
869 !
870 !- Check if the summation of the small delta t >= big delta t
871 ! (t_del_tv) (del_tv) (dtb)
873 t_del_tv=t_del_tv+del_tv
875 if ( t_del_tv .ge. dtb ) then
876 notlast=.false.
877 del_tv=dtb+del_tv-t_del_tv
878 endif
880 !
881 fluxin=0.
882 do k=max_q,min_q,-1
883 fluxout=rho(k)*vtgold(k)*qgz(k)
884 flux=(fluxin-fluxout)/rho(k)/dzw(k)
885 qgz(k)=qgz(k)+del_tv*flux
886 qgz(k)=amax1(0.,qgz(k))
887 fluxin=fluxout
888 enddo
889 if (min_q .eq. 1) then
890 pptgraul=pptgraul+fluxin*del_tv
891 else
892 qgz(min_q-1)=qgz(min_q-1)+del_tv* &
893 fluxin/rho(min_q-1)/dzw(min_q-1)
894 endif
895 !
896 else
897 notlast=.false.
898 endif
899 !
900 ENDDO
902 !
903 !-- cloud ice (03/21/02) follow Vaughan T.J. Phillips at GFDL
904 !
905 t_del_tv=0.
906 del_tv=dtb
907 notlast=.true.
908 !
909 DO while (notlast)
910 !
911 min_q=kte
912 max_q=kts-1
913 !
914 do k=kts,kte-1
915 if (qiz(k) .gt. 1.0e-8) then
916 min_q=min0(min_q,k)
917 max_q=max0(max_q,k)
918 vtiold(k)= 3.29 * (rho(k)* qiz(k))** 0.16 ! Heymsfield and Donner
919 if (k .eq. 1) then
920 del_tv=amin1(del_tv,0.9*(zz(k)-zsfc)/vtiold(k))
921 else
922 del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtiold(k))
923 endif
924 else
925 vtiold(k)=0.
926 endif
927 enddo
929 if (max_q .ge. min_q) then
930 !
931 !
932 !- Check if the summation of the small delta t >= big delta t
933 ! (t_del_tv) (del_tv) (dtb)
935 t_del_tv=t_del_tv+del_tv
937 if ( t_del_tv .ge. dtb ) then
938 notlast=.false.
939 del_tv=dtb+del_tv-t_del_tv
940 endif
942 fluxin=0.
943 do k=max_q,min_q,-1
944 fluxout=rho(k)*vtiold(k)*qiz(k)
945 flux=(fluxin-fluxout)/rho(k)/dzw(k)
946 qiz(k)=qiz(k)+del_tv*flux
947 qiz(k)=amax1(0.,qiz(k))
948 fluxin=fluxout
949 enddo
950 if (min_q .eq. 1) then
951 pptice=pptice+fluxin*del_tv
952 else
953 qiz(min_q-1)=qiz(min_q-1)+del_tv* &
954 fluxin/rho(min_q-1)/dzw(min_q-1)
955 endif
956 !
957 else
958 notlast=.false.
959 endif
960 !
961 ENDDO
962 do k=kts,kte-1 !sg beg
963 precrz(k)=rho(k)*vtrold(k)*qrz(k)
964 preciz(k)=rho(k)*vtiold(k)*qiz(k)
965 precsz(k)=rho(k)*vtsold(k)*qsz(k)
966 precgz(k)=rho(k)*vtgold(k)*qgz(k)
967 enddo !sg end
968 precrz(kte)=0. !wig - top level never set for vtXold vars
969 preciz(kte)=0. !wig
970 precsz(kte)=0. !wig
971 precgz(kte)=0. !wig
974 ! Microphpysics processes
975 !
976 DO 2000 k=kts,kte
977 !
978 !***********************************************************************
979 !***** diagnose mixing ratios (qrz,qsz), terminal *****
980 !***** velocities (vtr,vts), and slope parameters in size *****
981 !***** distribution(xlambdar,xlambdas) of rain and snow *****
982 !***** follows Nagata and Ogura, 1991, MWR, 1309-1337. Eq (A7) *****
983 !***********************************************************************
984 !
985 !**** assuming no cloud water can exist in the top two levels due to
986 !**** radiation consideration
987 !
988 !! if
989 !! unsaturated,
990 !! no cloud water, rain, ice, snow and graupel
991 !! then
992 !! skip these processes and jump to line 2000
993 !
994 !
995 tmp=qiz(k)+qlz(k)+qsz(k)+qrz(k)+qgz(k)*gindex
996 if( qvz(k)+qlz(k)+qiz(k) .lt. qsiz(k) &
997 .and. tmp .eq. 0.0 ) go to 2000
999 !! calculate terminal velocity of rain
1000 !
1001 if (qrz(k) .gt. 1.0e-8) then
1002 tmp1=sqrt(pi*rhowater*xnor*orho(k)/qrz(k))
1003 xlambdar(k)=sqrt(tmp1)
1004 olambdar(k)=1.0/xlambdar(k)
1005 vtrold(k)=o6*consta*gambp4*sqrho(k)*olambdar(k)**constb
1006 else
1007 vtrold(k)=0.
1008 olambdar(k)=0.
1009 endif
1010 !
1011 ! if (qrz(k) .gt. 1.0e-12) then
1012 if (qrz(k) .gt. 1.0e-8) then
1013 tmp1=sqrt(pi*rhowater*xnor*orho(k)/qrz(k))
1014 xlambdar(k)=sqrt(tmp1)
1015 olambdar(k)=1.0/xlambdar(k)
1016 vtr(k)=o6*consta*gambp4*sqrho(k)*olambdar(k)**constb
1017 else
1018 vtr(k)=0.
1019 olambdar(k)=0.
1020 endif
1021 !
1022 !! calculate terminal velocity of snow
1023 !
1024 if (qsz(k) .gt. 1.0e-8) then
1025 tmp1=sqrt(pi*rhosnow*xnos*orho(k)/qsz(k))
1026 xlambdas(k)=sqrt(tmp1)
1027 olambdas(k)=1.0/xlambdas(k)
1028 vtsold(k)=o6*constc*gamdp4*sqrho(k)*olambdas(k)**constd
1029 else
1030 vtsold(k)=0.
1031 olambdas(k)=0.
1032 endif
1033 !
1034 ! if (qsz(k) .gt. 1.0e-12) then
1035 if (qsz(k) .gt. 1.0e-8) then
1036 tmp1=sqrt(pi*rhosnow*xnos*orho(k)/qsz(k))
1037 xlambdas(k)=sqrt(tmp1)
1038 olambdas(k)=1.0/xlambdas(k)
1039 vts(k)=o6*constc*gamdp4*sqrho(k)*olambdas(k)**constd
1040 else
1041 vts(k)=0.
1042 olambdas(k)=0.
1043 endif
1044 !
1045 !! calculate terminal velocity of graupel
1046 !
1047 if (qgz(k) .gt. 1.0e-8) then
1048 tmp1=sqrt( pi*rhograul*xnog*orho(k)/qgz(k))
1049 xlambdag(k)=sqrt(tmp1)
1050 olambdag(k)=1.0/xlambdag(k)
1051 term0=sqrt(4.*grav*rhograul*0.33334*orho(k)*ocdrag)
1052 vtgold(k)=o6*gam4pt5*term0*sqrt(olambdag(k))
1053 else
1054 vtgold(k)=0.
1055 olambdag(k)=0.
1056 endif
1057 !
1058 ! if (qgz(k) .gt. 1.0e-12) then
1059 if (qgz(k) .gt. 1.0e-8) then
1060 tmp1=sqrt( pi*rhograul*xnog*orho(k)/qgz(k))
1061 xlambdag(k)=sqrt(tmp1)
1062 olambdag(k)=1.0/xlambdag(k)
1063 term0=sqrt(4.*grav*rhograul*0.33334*orho(k)*ocdrag)
1064 vtg(k)=o6*gam4pt5*term0*sqrt(olambdag(k))
1065 else
1066 vtg(k)=0.
1067 olambdag(k)=0.
1068 endif
1069 !
1070 !***********************************************************************
1071 !***** compute viscosity,difusivity,thermal conductivity, and ******
1072 !***** Schmidt number ******
1073 !***********************************************************************
1074 !c------------------------------------------------------------------
1075 !c viscmu: dynamic viscosity of air kg/m/s
1076 !c visc: kinematic viscosity of air = viscmu/rho (m2/s)
1077 !c avisc=1.49628e-6 kg/m/s=1.49628e-5 g/cm/s
1078 !c viscmu=1.718e-5 kg/m/s in RH
1079 !c diffwv: Diffusivity of water vapor in air
1080 !c adiffwv = 8.7602e-5 (8.794e-5 in MM5) kgm/s3
1081 !c = 8.7602 (8.794 in MM5) gcm/s3
1082 !c diffwv(k)=2.26e-5 m2/s
1083 !c schmidt: Schmidt number=visc/diffwv
1084 !c xka: thermal conductivity of air J/m/s/K (Kgm/s3/K)
1085 !c xka(k)=2.43e-2 J/m/s/K in RH
1086 !c axka=1.4132e3 (1.414e3 in MM5) m2/s2/k = 1.4132e7 cm2/s2/k
1087 !c------------------------------------------------------------------
1089 viscmu(k)=avisc*tem(k)**1.5/(tem(k)+120.0)
1090 visc(k)=viscmu(k)*orho(k)
1091 diffwv(k)=adiffwv*tem(k)**1.81*oprez(k)
1092 schmidt(k)=visc(k)/diffwv(k)
1093 xka(k)=axka*viscmu(k)
1095 if (tem(k) .lt. 273.15) then
1097 !
1098 !***********************************************************************
1099 !********* snow production processes for T < 0 C **********
1100 !***********************************************************************
1101 !c
1102 !c (1) ICE CRYSTAL AGGREGATION TO SNOW (Psaut): Lin (21)
1103 !c! psaut=alpha1*(qi-qi0)
1104 !c! alpha1=1.0e-3*exp(0.025*(T-T0))
1105 !c
1106 ! alpha1=1.0e-3*exp( 0.025*temcc(k) )
1108 alpha1=1.0e-3*exp( 0.025*temcc(k) )
1109 !
1110 if(temcc(k) .lt. -20.0) then
1111 tmp1=-7.6+4.0*exp( -0.2443e-3*(abs(temcc(k))-20)**2.455 )
1112 qic=1.0e-3*exp(tmp1)*orho(k)
1113 else
1114 qic=qi0
1115 end if
1116 !testing
1117 ! tmp1=amax1( 0.0,alpha1*(qiz(k)-qic) )
1118 ! psaut(k)=amin1( tmp1,qizodt(k) )
1120 tmp1=odtb*(qiz(k)-qic)*(1.0-exp(-alpha1*dtb))
1121 psaut(k)=amax1( 0.0,tmp1 )
1123 !c
1124 !c (2) BERGERON PROCESS TRANSFER OF CLOUD WATER TO SNOW (Psfw)
1125 !c this process only considered when -31 C < T < 0 C
1126 !c Lin (33) and Hsie (17)
1127 !c
1128 !c!
1129 !c! parama1 and parama2 functions must be user supplied
1130 !c!
1132 ! testing
1133 if( qlz(k) .gt. 1.0e-10 ) then
1134 temc1=amax1(-30.99,temcc(k))
1135 ! print*,'temc1',temc1,qlz(k)
1136 a1=parama1( temc1 )
1137 a2=parama2( temc1 )
1138 tmp1=1.0-a2
1139 !! change unit from cgs to mks
1140 a1=a1*0.001**tmp1
1141 !c! dtberg is the time needed for a crystal to grow from 40 to 50 um
1142 !c ! odtberg=1.0/dtberg
1143 odtberg=(a1*tmp1)/(xmi50**tmp1-xmi40**tmp1)
1144 !
1145 !c! compute terminal velocity of a 50 micron ice cystal
1146 !
1147 vti50=constc*di50**constd*sqrho(k)
1148 !
1149 eiw=1.0
1150 save1=a1*xmi50**a2
1151 save2=0.25*pi*eiw*rho(k)*di50*di50*vti50
1152 !
1153 tmp2=( save1 + save2*qlz(k) )
1154 !
1155 !! maximum number of 50 micron crystals limited by the amount
1156 !! of supercool water
1157 !
1158 xni50mx=qlzodt(k)/tmp2
1159 !
1160 !! number of 50 micron crystals produced
1161 !
1162 !
1163 xni50=qiz(k)*( 1.0-exp(-dtb*odtberg) )/xmi50
1164 xni50=amin1(xni50,xni50mx)
1165 !
1166 tmp3=odtb*tmp2/save2*( 1.0-exp(-save2*xni50*dtb) )
1167 psfw(k)=amin1( tmp3,qlzodt(k) )
1168 !testing
1169 ! psfw(k)=0.
1171 !0915 if( temcc(k).gt.-30.99 ) then
1172 !0915 a1=parama1( temcc(k) )
1173 !0915 a2=parama2( temcc(k) )
1174 !0915 tmp1=1.0-a2
1175 !! change unit from cgs to mks
1176 !0915 a1=a1*0.001**tmp1
1178 !c! dtberg is the time needed for a crystal to grow from 40 to 50 um
1179 !c! odtberg=1.0/dtberg
1180 !0915 odtberg=(a1*tmp1)/(xmi50**tmp1-xmi40**tmp1)
1182 !c! number of 50 micron crystals produced
1183 !0915 xni50=qiz(k)*dtb*odtberg/xmi50
1185 !c! need to calculate the terminal velocity of a 50 micron
1186 !c! ice cystal
1187 !0915 vti50=constc*di50**constd*sqrho(k)
1188 !0915 eiw=1.0
1189 !0915 tmp2=xni50*( a1*xmi50**a2 + &
1190 !0915 0.25*qlz(k)*pi*eiw*rho(k)*di50*di50*vti50 )
1191 !0915 psfw(k)=amin1( tmp2,qlzodt(k) )
1192 !0915 psfw(k)=0.
1193 !c
1194 !c (3) REDUCTION OF CLOUD ICE BY BERGERON PROCESS (Psfi): Lin (34)
1195 !c this process only considered when -31 C < T < 0 C
1196 !c
1197 tmp1=xni50*xmi50-psfw(k)
1198 psfi(k)=amin1(tmp1,qizodt(k))
1199 ! testing
1200 ! psfi(k)=0.
1201 end if
1202 !
1204 !0915 tmp1=qiz(k)*odtberg
1205 !0915 psfi(k)=amin1(tmp1,qizodt(k))
1206 ! testing
1207 !0915 psfi(k)=0.
1208 !0915 end if
1209 !
1210 if(qrz(k) .le. 0.0) go to 1000
1211 !
1212 ! Processes (4) and (5) only need when qrz > 0.0
1213 !
1214 !c
1215 !c (4) CLOUD ICE ACCRETION BY RAIN (Praci): Lin (25)
1216 !c may produce snow or graupel
1217 !c
1218 eri=1.0
1219 !0915 tmp1=qiz(k)*pio4*eri*xnor*consta*sqrho(k)
1220 !0915 tmp2=tmp1*gambp3*olambdar(k)**bp3
1221 !0915 praci(k)=amin1( tmp2,qizodt(k) )
1223 save1=pio4*eri*xnor*consta*sqrho(k)
1224 tmp1=save1*gambp3*olambdar(k)**bp3
1225 praci(k)=qizodt(k)*( 1.0-exp(-tmp1*dtb) )
1227 !c
1228 !c (5) RAIN ACCRETION BY CLOUD ICE (Piacr): Lin (26)
1229 !c
1230 !0915 tmp2=tmp1*rho(k)*pio6*rhowater*gambp6*oxmi* &
1231 !0915 olambdar(k)**bp6
1232 !0915 piacr(k)=amin1( tmp2,qrzodt(k) )
1234 tmp2=qiz(k)*save1*rho(k)*pio6*rhowater*gambp6*oxmi* &
1235 olambdar(k)**bp6
1236 piacr(k)=amin1( tmp2,qrzodt(k) )
1238 !
1239 1000 continue
1240 !
1241 if(qsz(k) .le. 0.0) go to 1200
1242 !
1243 ! Compute the following processes only when qsz > 0.0
1244 !
1245 !c
1246 !c (6) ICE CRYSTAL ACCRETION BY SNOW (Psaci): Lin (22)
1247 !c
1248 esi=exp( 0.025*temcc(k) )
1249 save1=pio4*xnos*constc*gamdp3*sqrho(k)* &
1250 olambdas(k)**dp3
1251 tmp1=esi*save1
1252 psaci(k)=qizodt(k)*( 1.0-exp(-tmp1*dtb) )
1254 !0915 tmp1=pio4*xnos*constc*gamdp3*sqrho(k)* &
1255 !0915 olambdas(k)**dp3
1256 !0915 tmp2=qiz(k)*esi*tmp1
1257 !0915 psaci(k)=amin1( tmp2,qizodt(k) )
1258 !c
1259 !c (7) CLOUD WATER ACCRETION BY SNOW (Psacw): Lin (24)
1260 !c
1261 esw=1.0
1262 tmp1=esw*save1
1263 psacw(k)=qlzodt(K)*( 1.0-exp(-tmp1*dtb) )
1265 !0915 tmp2=qlz(k)*esw*tmp1
1266 !0915 psacw(k)=amin1( tmp2,qlzodt(k) )
1267 !c
1268 !c (8) DEPOSITION/SUBLIMATION OF SNOW (Psdep/Pssub): Lin (31)
1269 !c includes consideration of ventilation effect
1270 !c
1271 !c abi=2*pi*(Si-1)/rho/(A"+B")
1272 !c
1273 tmpa=rvapor*xka(k)*tem(k)*tem(k)
1274 tmpb=xls*xls*rho(k)*qsiz(k)*diffwv(k)
1275 tmpc=tmpa*qsiz(k)*diffwv(k)
1276 abi=2.0*pi*(qvoqsiz(k)-1.0)*tmpc/(tmpa+tmpb)
1277 !
1278 !c vf1s,vf2s=ventilation factors for snow
1279 !c vf1s=0.78,vf2s=0.31 in LIN
1280 !
1281 tmp1=constc*sqrho(k)*olambdas(k)**dp5/visc(k)
1282 tmp2=abi*xnos*( vf1s*olambdas(k)*olambdas(k)+ &
1283 vf2s*schmidt(k)**0.33334*gamdp5o2*sqrt(tmp1) )
1284 tmp3=odtb*( qvz(k)-qsiz(k) )
1285 !
1286 if( tmp2 .le. 0.0) then
1287 tmp2=amax1( tmp2,tmp3)
1288 pssub(k)=amax1( tmp2,-qszodt(k) )
1289 psdep(k)=0.0
1290 else
1291 psdep(k)=amin1( tmp2,tmp3 )
1292 pssub(k)=0.0
1293 end if
1295 !0915 psdep(k)=amax1(0.0,tmp2)
1296 !0915 pssub(k)=amin1(0.0,tmp2)
1297 !0915 pssub(k)=amax1( pssub(k),-qszodt(k) )
1298 !
1299 if(qrz(k) .le. 0.0) go to 1200
1300 !
1301 ! Compute processes (9) and (10) only when qsz > 0.0 and qrz > 0.0
1302 !
1303 !c
1304 !c (9) ACCRETION OF SNOW BY RAIN (Pracs): Lin (27)
1305 !c
1306 esr=1.0
1307 tmpa=olambdar(k)*olambdar(k)
1308 tmpb=olambdas(k)*olambdas(k)
1309 tmpc=olambdar(k)*olambdas(k)
1310 tmp1=pi*pi*esr*xnor*xnos*abs( vtr(k)-vts(k) )*orho(k)
1311 tmp2=tmpb*tmpb*olambdar(k)*(5.0*tmpb+2.0*tmpc+0.5*tmpa)
1312 tmp3=tmp1*rhosnow*tmp2
1313 pracs(k)=amin1( tmp3,qszodt(k) )
1314 !c
1315 !c (10) ACCRETION OF RAIN BY SNOW (Psacr): Lin (28)
1316 !c
1317 tmp3=tmpa*tmpa*olambdas(k)*(5.0*tmpa+2.0*tmpc+0.5*tmpb)
1318 tmp4=tmp1*rhowater*tmp3
1319 psacr(k)=amin1( tmp4,qrzodt(k) )
1320 !
1321 1200 continue
1322 !
1323 else
1324 !
1325 !***********************************************************************
1326 !********* snow production processes for T > 0 C **********
1327 !***********************************************************************
1328 !
1329 if (qsz(k) .le. 0.0) go to 1400
1330 !c
1331 !c (1) CLOUD WATER ACCRETION BY SNOW (Psacw): Lin (24)
1332 !c
1333 esw=1.0
1335 tmp1=esw*pio4*xnos*constc*gamdp3*sqrho(k)* &
1336 olambdas(k)**dp3
1337 psacw(k)=qlzodt(k)*( 1.0-exp(-tmp1*dtb) )
1339 !0915 tmp1=pio4*xnos*constc*gamdp3*sqrho(k)* &
1340 !0915 olambdas(k)**dp3
1341 !0915 tmp2=qlz(k)*esw*tmp1
1342 !0915 psacw(k)=amin1( tmp2,qlzodt(k) )
1343 !c
1344 !c (2) ACCRETION OF RAIN BY SNOW (Psacr): Lin (28)
1345 !c
1346 esr=1.0
1347 tmpa=olambdar(k)*olambdar(k)
1348 tmpb=olambdas(k)*olambdas(k)
1349 tmpc=olambdar(k)*olambdas(k)
1350 tmp1=pi*pi*esr*xnor*xnos*abs( vtr(k)-vts(k) )*orho(k)
1351 tmp2=tmpa*tmpa*olambdas(k)*(5.0*tmpa+2.0*tmpc+0.5*tmpb)
1352 tmp3=tmp1*rhowater*tmp2
1353 psacr(k)=amin1( tmp3,qrzodt(k) )
1354 !c
1355 !c (3) MELTING OF SNOW (Psmlt): Lin (32)
1356 !c Psmlt is negative value
1357 !
1358 delrs=rs0(k)-qvz(k)
1359 term1=2.0*pi*orho(k)*( xlv*diffwv(k)*rho(k)*delrs- &
1360 xka(k)*temcc(k) )
1361 tmp1=constc*sqrho(k)*olambdas(k)**dp5/visc(k)
1362 tmp2=xnos*( vf1s*olambdas(k)*olambdas(k)+ &
1363 vf2s*schmidt(k)**0.33334*gamdp5o2*sqrt(tmp1) )
1364 tmp3=term1*oxlf*tmp2-cwoxlf*temcc(k)*( psacw(k)+psacr(k) )
1365 tmp4=amin1(0.0,tmp3)
1366 psmlt(k)=amax1( tmp4,-qszodt(k) )
1367 !c
1368 !c (4) EVAPORATION OF MELTING SNOW (Psmltevp): HR (A27)
1369 !c but use Lin et al. coefficience
1370 !c Psmltevp is a negative value
1371 !c
1372 tmpa=rvapor*xka(k)*tem(k)*tem(k)
1373 tmpb=xlv*xlv*rho(k)*qswz(k)*diffwv(k)
1374 tmpc=tmpa*qswz(k)*diffwv(k)
1375 tmpd=amin1( 0.0,(qvoqswz(k)-0.90)*qswz(k)*odtb )
1377 ! abr=2.0*pi*(qvoqswz(k)-1.0)*tmpc/(tmpa+tmpb)
1379 abr=2.0*pi*(qvoqswz(k)-0.90)*tmpc/(tmpa+tmpb)
1380 !
1381 !**** allow evaporation to occur when RH less than 90%
1382 !**** here not using 100% because the evaporation cooling
1383 !**** of temperature is not taking into account yet; hence,
1384 !**** the qsw value is a little bit larger. This will avoid
1385 !**** evaporation can generate cloud.
1386 !
1387 !c vf1s,vf2s=ventilation factors for snow
1388 !c vf1s=0.78,vf2s=0.31 in LIN
1389 !
1390 tmp1=constc*sqrho(k)*olambdas(k)**dp5/visc(k)
1391 tmp2=abr*xnos*( vf1s*olambdas(k)*olambdas(k)+ &
1392 vf2s*schmidt(k)**0.33334*gamdp5o2*sqrt(tmp1) )
1393 tmp3=amin1(0.0,tmp2)
1394 tmp3=amax1( tmp3,tmpd )
1395 psmltevp(k)=amax1( tmp3,-qszodt(k) )
1396 1400 continue
1397 !
1398 end if
1400 !***********************************************************************
1401 !********* rain production processes **********
1402 !***********************************************************************
1403 !
1404 !c
1405 !c (1) AUTOCONVERSION OF RAIN (Praut): RH
1406 !sg: begin
1407 if(flag_qndrop)then
1408 if( qndropz(k) >= 1. ) then
1409 ! Liu et al. autoconversion scheme
1410 rhocgs=rho(k)*1.e-3
1411 liqconc=rhocgs*qlz(k)
1412 capn=rhocgs*qndropz(k)
1413 ! rate function
1414 if(liqconc.gt.1.e-10)then
1415 p0=kappa*beta/capn*(liqconc*liqconc*liqconc)
1416 xc=9.7d-17*capn*sqrt(capn)/(liqconc*liqconc)
1417 ! Calculate autoconversion rate (g/g/s)
1418 if(xc.lt.10.)then
1419 praut(k)=p0/rhocgs*0.5d0*(xc*xc+2*xc+2.0d0)* &
1420 (1.0d0+xc)*dexp(-2.0d0*xc)
1421 endif
1422 endif
1423 endif
1424 else
1425 !sg: end
1426 !c araut=afa*rho
1427 !c afa=0.001 Rate coefficient for autoconvergence
1428 !c
1429 !c araut=1.0e-3
1430 !c
1431 araut=0.001
1432 !testing
1433 ! tmp1=amax1( 0.0,araut*(qlz(k)-ql0) )
1434 ! praut(k)=amin1( tmp1,qlzodt(k) )
1435 tmp1=odtb*(qlz(k)-ql0)*( 1.0-exp(-araut*dtb) )
1436 praut(k)=amax1( 0.0,tmp1 )
1437 endif !sg
1439 !c
1440 !c (2) ACCRETION OF CLOUD WATER BY RAIN (Pracw): Lin (51)
1441 !c
1442 erw=1.0
1443 ! tmp1=qlz(k)*pio4*erw*xnor*consta*sqrho(k)
1444 ! tmp2=tmp1*gambp3*olambdar(k)**bp3
1445 ! pracw(k)=amin1( tmp2,qlzodt(k) )
1447 tmp1=pio4*erw*xnor*consta*sqrho(k)* &
1448 gambp3*olambdar(k)**bp3
1449 pracw(k)=qlzodt(k)*( 1.0-exp(-tmp1*dtb) )
1451 !c
1452 !c (3) EVAPORATION OF RAIN (Prevp): Lin (52)
1453 !c Prevp is negative value
1454 !c
1455 !c Sw=qvoqsw : saturation ratio
1456 !c
1457 tmpa=rvapor*xka(k)*tem(k)*tem(k)
1458 tmpb=xlv*xlv*rho(k)*qswz(k)*diffwv(k)
1459 tmpc=tmpa*qswz(k)*diffwv(k)
1460 tmpd=amin1(0.0,(qvoqswz(k)-0.90)*qswz(k)*odtb)
1461 !
1462 ! abr=2.0*pi*(qvoqswz(k)-1.0)*tmpc/(tmpa+tmpb)
1464 abr=2.0*pi*(qvoqswz(k)-0.90)*tmpc/(tmpa+tmpb)
1465 !
1466 !c vf1r,vf2r=ventilation factors for rain
1467 !c vf1r=0.78,vf2r=0.31 in RH, LIN and MM5
1468 !
1469 vf1r=0.78
1470 vf2r=0.31
1471 tmp1=consta*sqrho(k)*olambdar(k)**bp5/visc(k)
1472 tmp2=abr*xnor*( vf1r*olambdar(k)*olambdar(k)+ &
1473 vf2r*schmidt(k)**0.33334*gambp5o2*sqrt(tmp1) )
1474 tmp3=amin1( 0.0,tmp2 )
1475 tmp3=amax1( tmp3,tmpd )
1476 prevp(k)=amax1( tmp3,-qrzodt(k) )
1478 !
1479 ! if(iout .gt. 0) write(20,*)'tmp1,tmp2,tmp3=',tmp1,tmp2,tmp3
1480 ! if(iout .gt. 0) write(20,*)'qlz,qiz,qrz=',qlz(k),qiz(k),qrz(k)
1481 ! if(iout .gt. 0) write(20,*)'tem,qsz,qvz=',tem(k),qsz(k),qvz(k)
1485 ! if (gindex .eq. 0.) goto 900
1486 !
1487 if (tem(k) .lt. 273.15) then
1488 !
1489 !
1490 !-- graupel
1491 !***********************************************************************
1492 !********* graupel production processes for T < 0 C **********
1493 !***********************************************************************
1494 !c
1495 !c (1) AUTOCONVERSION OF SNOW TO FORM GRAUPEL (Pgaut): Lin (37)
1496 !c pgaut=alpha2*(qsz-qs0)
1497 !c qs0=6.0E-4
1498 !c alpha2=1.0e-3*exp(0.09*temcc(k)) Lin (38)
1499 !
1500 alpha2=1.0e-3*exp(0.09*temcc(k))
1501 !
1503 ! testing
1504 ! tmp1=alpha2*(qsz(k)-qs0)
1505 ! tmp1=amax1(0.0,tmp1)
1506 ! pgaut(k)=amin1( tmp1,qszodt(k) )
1508 tmp1=odtb*(qsz(k)-qs0)*(1.0-exp(-alpha2*dtb))
1509 pgaut(k)=amax1( 0.0,tmp1 )
1511 !c
1512 !c (2) FREEZING OF RAIN TO FORM GRAUPEL (Pgfr): Lin (45)
1513 !c positive value
1514 !c Constant in Bigg freezing Aplume=Ap=0.66 /k
1515 !c Constant in raindrop freezing equ. Bplume=Bp=100./m/m/m/s
1516 !
1518 if (qrz(k) .gt. 1.e-8 ) then
1519 Bp=100.
1520 Ap=0.66
1521 tmp1=olambdar(k)*olambdar(k)*olambdar(k)
1522 tmp2=20.*pi*pi*Bp*xnor*rhowater*orho(k)* &
1523 (exp(-Ap*temcc(k))-1.0)*tmp1*tmp1*olambdar(k)
1524 Pgfr(k)=amin1( tmp2,qrzodt(k) )
1525 else
1526 Pgfr(k)=0
1527 endif
1529 !c
1530 !c if (qgz(k) = 0.0) skip the other step below about graupel
1531 !c
1532 if (qgz(k) .eq. 0.0) goto 4000
1534 !c
1535 !c Comparing Pgwet(wet process) and Pdry(dry process),
1536 !c we will pick up the small one.
1537 !c
1539 !c ---------------
1540 !c | dry processes |
1541 !c ---------------
1542 !c
1543 !c (3) ACCRETION OF CLOUD WATER BY GRAUPEL (Pgacw): Lin (40)
1544 !c egw=1.0
1545 !c Cdrag=0.6 drag coefficients for hairstone
1546 !c constg=sqrt(4.*grav*rhograul*0.33334*orho(k)/Cdrag)
1547 !c
1548 egw=1.0
1549 constg=sqrt(4.*grav*rhograul*0.33334*orho(k)*oCdrag)
1550 tmp1=pio4*xnog*gam3pt5*constg*olambdag(k)**3.5
1551 tmp2=qlz(k)*egw*tmp1
1552 Pgacw(k)=amin1( tmp2,qlzodt(k) )
1553 !c
1554 !c (4) ACCRETION OF ICE CRYSTAL BY GRAUPEL (Pgaci): Lin (41)
1555 !c egi=1. for wet growth
1556 !c egi=0.1 for dry growth
1557 !c
1558 egi=0.1
1559 tmp2=qiz(k)*egi*tmp1
1560 pgaci(k)=amin1( tmp2,qizodt(k) )
1563 !c
1564 !c (5) ACCRETION OF SNOW BY GRAUPEL (Pgacs) : Lin (29)
1565 !c Compute processes (6) only when qsz > 0.0 and qgz > 0.0
1566 !c
1567 egs=exp(0.09*temcc(k))
1568 tmpa=olambdas(k)*olambdas(k)
1569 tmpb=olambdag(k)*olambdag(k)
1570 tmpc=olambdas(k)*olambdag(k)
1571 tmp1=pi*pi*xnos*xnog*abs( vts(k)-vtg(k) )*orho(k)
1572 tmp2=tmpa*tmpa*olambdag(k)*(5.0*tmpa+2.0*tmpc+0.5*tmpb)
1573 tmp3=tmp1*egs*rhosnow*tmp2
1574 Pgacs(k)=amin1( tmp3,qszodt(k) )
1577 !c
1578 !c (6) ACCRETION OF RAIN BY GRAUPEL (Pgacr): Lin (42)
1579 !c Compute processes (6) only when qrz > 0.0 and qgz > 0.0
1580 !c egr=1.
1581 !c
1582 egr=1.
1583 tmpa=olambdar(k)*olambdar(k)
1584 tmpb=olambdag(k)*olambdag(k)
1585 tmpc=olambdar(k)*olambdag(k)
1586 tmp1=pi*pi*xnor*xnog*abs( vtr(k)-vtg(k) )*orho(k)
1587 tmp2=tmpa*tmpa*olambdag(k)*(5.0*tmpa+2.0*tmpc+0.5*tmpb)
1588 tmp3=tmp1*egr*rhowater*tmp2
1589 pgacr(k)=amin1( tmp3,qrzodt(k) )
1591 !c
1592 !c (7) Calculate total dry process effect Pdry(k)
1593 !c
1594 Pdry(k)=Pgacw(k)+pgaci(k)+Pgacs(k)+pgacr(k)
1596 !c ---------------
1597 !c | wet processes |
1598 !c ---------------
1599 !c
1600 !c (3) ACCRETION OF ICE CRYSTAL BY GRAUPEL (Pgacip): Lin (41)
1601 !c egi=1. for wet growth
1602 !c egi=0.1 for dry growth
1603 !c
1604 tmp2=10.*pgaci(k)
1605 pgacip(k)=amin1( tmp2,qizodt(k) )
1607 !c
1608 !c (4) ACCRETION OF SNOW BY GRAUPEL ((Pgacsp) : Lin (29)
1609 !c Compute processes (6) only when qsz > 0.0 and qgz > 0.0
1610 !c egs=exp(0.09*(tem(k)-273.15)) when T < 273.15 k
1611 !c
1612 tmp3=Pgacs(k)*1.0/egs
1613 Pgacsp(k)=amin1( tmp3,qszodt(k) )
1615 !c
1616 !c (5) WET GROWTH OF GRAUPEL (Pgwet) : Lin (43)
1617 !c may involve Pgacs or Pgaci and
1618 !c must include PPgacw or Pgacr, or both.
1619 !c ( The amount of Pgacw which is not able
1620 !c to freeze is shed to rain. )
1621 IF(temcc(k).gt.-40.)THEN
1623 term0=constg*olambdag(k)**5.5/visc(k)
1625 !c
1626 !c vf1s,vf2s=ventilation factors for graupel
1627 !c vf1s=0.78,vf2s=0.31 in LIN
1628 !c Cdrag=0.6 drag coefficient for hairstone
1629 !c constg2=vf1s*olambdag(k)*olambdag(k)+
1630 !c vf2s*schmidt(k)**0.33334*gam2pt75*sqrt(term0)
1632 delrs=rs0(k)-qvz(k)
1633 tmp0=1./(xlf+cw*temcc(k))
1634 tmp1=2.*pi*xnog*(rho(k)*xlv*diffwv(k)*delrs-xka(k)* &
1635 temcc(k))*orho(k)*tmp0
1636 constg2=vf1s*olambdag(k)*olambdag(k)+ &
1637 vf2s*schmidt(k)**0.33334*gam2pt75*sqrt(term0)
1638 tmp3=tmp1*constg2+(Pgacip(k)+Pgacsp(k))* &
1639 (1-Ci*temcc(k)*tmp0)
1640 tmp3=amax1(0.0,tmp3)
1641 Pgwet(k)=amax1(tmp3,qlzodt(k)+qszodt(k)+qizodt(k) )
1643 !c
1644 !c Comparing Pgwet(wet process) and Pdry(dry process),
1645 !c we will apply the small one.
1646 !c if dry processes then delta4=1.0
1647 !c if wet processes then delta4=0.0
1648 !
1649 if ( Pdry(k) .lt. Pgwet(k) ) then
1650 delta4=1.0
1651 else
1652 delta4=0.0
1653 endif
1654 ELSE
1655 delta4=1.0
1656 ENDIF
1658 !c
1659 !c
1660 !c (6) Pgacrp(k)=Pgwet(k)-Pgacw(k)-Pgacip(k)-Pgacsp(k)
1661 !c if Pgacrp(k) > 0. then some of the rain is frozen to hail
1662 !c if Pgacrp(k) < 0. then some of the cloud water collected
1663 !c by the hail is unable to freeze and is
1664 !c shed as rain.
1665 !c
1666 Pgacrp(k)=Pgwet(k)-Pgacw(k)-Pgacip(k)-Pgacsp(k)
1668 !c
1669 !c (8) DEPOSITION/SUBLIMATION OF GRAUPEL (Pgdep/Pgsub): Lin (46)
1670 !c includes ventilation effect
1671 !c constg=sqrt(4.*grav*rhograul*0.33334*orho(k)/Cdrag)
1672 !c constg2=vf1s*olambdag(k)*olambdag(k)+
1673 !c vf2s*schmidt(k)**0.33334*gam2pt75*constg
1674 !c
1675 !c abg=2*pi*(Si-1)/rho/(A"+B")
1676 !c
1677 tmpa=rvapor*xka(k)*tem(k)*tem(k)
1678 tmpb=xls*xls*rho(k)*qsiz(k)*diffwv(k)
1679 tmpc=tmpa*qsiz(k)*diffwv(k)
1680 abg=2.0*pi*(qvoqsiz(k)-1.0)*tmpc/(tmpa+tmpb)
1681 !c
1682 !c vf1s,vf2s=ventilation factors for graupel
1683 !c vf1s=0.78,vf2s=0.31 in LIN
1684 !c Cdrag=0.6 drag coefficient for hairstone
1685 !c
1686 term0=constg*olambdag(k)**5.5/visc(k)
1687 constg2=vf1s*olambdag(k)*olambdag(k)+ &
1688 vf2s*schmidt(k)**0.33334*gam2pt75*sqrt(term0)
1689 tmp2=abg*xnog*constg2
1690 pgdep(k)=amax1(0.0,tmp2)
1691 pgsub(k)=amin1(0.0,tmp2)
1692 pgsub(k)=amax1( pgsub(k),-qgzodt(k) )
1694 4000 continue
1695 else
1696 !
1697 !***********************************************************************
1698 !********* graupel production processes for T > 0 C **********
1699 !***********************************************************************
1700 !
1701 !c
1702 !c (1) ACCRETION OF CLOUD WATER BY GRAUPEL (Pgacw): Lin (40)
1703 !c egw=1.0
1704 !c Cdrag=0.6 drag coefficients for hairstone
1705 !c constg=sqrt(4.*grav*rhograul*0.33334*orho(k)/Cdrag)
1707 egw=1.0
1708 constg=sqrt(4.*grav*rhograul*0.33334*orho(k)*oCdrag)
1709 tmp1=pio4*xnog*gam3pt5*constg*olambdag(k)**3.5
1710 tmp2=qlz(k)*egw*tmp1
1711 Pgacw(k)=amin1( tmp2,qlzodt(k) )
1713 !c
1714 !c (2) ACCRETION OF RAIN BY GRAUPEL (Pgacr): Lin (42)
1715 !c Compute processes (5) only when qrz > 0.0 and qgz > 0.0
1716 !c egr=1.
1717 !c
1718 egr=1.
1719 tmpa=olambdar(k)*olambdar(k)
1720 tmpb=olambdag(k)*olambdag(k)
1721 tmpc=olambdar(k)*olambdag(k)
1722 tmp1=pi*pi*xnor*xnog*abs( vtr(k)-vtg(k) )*orho(k)
1723 tmp2=tmpa*tmpa*olambdag(k)*(5.0*tmpa+2.0*tmpc+0.5*tmpb)
1724 tmp3=tmp1*egr*rhowater*tmp2
1725 pgacr(k)=amin1( tmp3,qrzodt(k) )
1728 !c
1729 !c (3) GRAUPEL MELTING TO FORM RAIN (Pgmlt): Lin (47)
1730 !c Pgmlt is negative value
1731 !c constg=sqrt(4.*grav*rhograul*0.33334*orho(k)/Cdrag)
1732 !c constg2=vf1s*olambdag(k)*olambdag(k)+
1733 !c vf2s*schmidt(k)**0.33334*gam2pt75*constg
1734 !c Cdrag=0.6 drag coefficients for hairstone
1735 !
1736 delrs=rs0(k)-qvz(k)
1737 term1=2.0*pi*orho(k)*( xlv*diffwv(k)*rho(k)*delrs- &
1738 xka(k)*temcc(k) )
1739 term0=sqrt(4.*grav*rhograul*0.33334*orho(k)*ocdrag) &
1740 *olambdag(k)**5.5/visc(k)
1742 constg2=vf1s*olambdag(k)*olambdag(k)+ &
1743 vf2s*schmidt(k)**0.33334*gam2pt75*sqrt(term0)
1744 tmp2=xnog*constg2
1745 tmp3=term1*oxlf*tmp2-cwoxlf*temcc(k)*( pgacw(k)+pgacr(k) )
1746 tmp4=amin1(0.0,tmp3)
1747 pgmlt(k)=amax1( tmp4,-qgzodt(k) )
1750 !c
1751 !c (4) EVAPORATION OF MELTING GRAUPEL (Pgmltevp) : HR (A19)
1752 !c but use Lin et al. coefficience
1753 !c Pgmltevp is a negative value
1754 !c abg=2.0*pi*(qvoqsiz(k)-1.0)*tmpc/(tmpa+tmpb)
1755 !c
1756 tmpa=rvapor*xka(k)*tem(k)*tem(k)
1757 tmpb=xlv*xlv*rho(k)*qswz(k)*diffwv(k)
1758 tmpc=tmpa*qswz(k)*diffwv(k)
1759 tmpd=amin1( 0.0,(qvoqswz(k)-0.90)*qswz(k)*odtb )
1761 !c
1762 !c abg=2*pi*(Si-1)/rho/(A"+B")
1763 !c
1764 abg=2.0*pi*(qvoqswz(k)-0.90)*tmpc/(tmpa+tmpb)
1765 !
1766 !**** allow evaporation to occur when RH less than 90%
1767 !**** here not using 100% because the evaporation cooling
1768 !**** of temperature is not taking into account yet; hence,
1769 !**** the qgw value is a little bit larger. This will avoid
1770 !**** evaporation can generate cloud.
1771 !
1772 !c vf1s,vf2s=ventilation factors for snow
1773 !c vf1s=0.78,vf2s=0.31 in LIN
1774 !c constg=sqrt(4.*grav*rhograul*0.33334*orho(k)/Cdrag)
1775 !c constg2=vf1s*olambdag(k)*olambdag(k)+
1776 !c vf2s*schmidt(k)**0.33334*gam2pt75*constg
1777 !
1778 tmp2=abg*xnog*constg2
1779 tmp3=amin1(0.0,tmp2)
1780 tmp3=amax1( tmp3,tmpd )
1781 pgmltevp(k)=amax1( tmp3,-qgzodt(k) )
1783 !c
1784 !c (5) ACCRETION OF SNOW BY GRAUPEL (Pgacs) : Lin (29)
1785 !c Compute processes (3) only when qsz > 0.0 and qgz > 0.0
1786 !c egs=1.0
1787 !c
1788 egs=1.
1789 tmpa=olambdas(k)*olambdas(k)
1790 tmpb=olambdag(k)*olambdag(k)
1791 tmpc=olambdas(k)*olambdag(k)
1792 tmp1=pi*pi*xnos*xnog*abs( vts(k)-vtg(k) )*orho(k)
1793 tmp2=tmpa*tmpa*olambdag(k)*(5.0*tmpa+2.0*tmpc+0.5*tmpb)
1794 tmp3=tmp1*egs*rhosnow*tmp2
1795 Pgacs(k)=amin1( tmp3,qszodt(k) )
1797 endif
1800 !
1801 900 continue
1803 !cc
1804 !c
1805 !c**********************************************************************
1806 !c***** combine all processes together and avoid negative *****
1807 !c***** water substances
1808 !***********************************************************************
1809 !c
1810 if ( temcc(k) .lt. 0.0) then
1811 !,delta4,1.-delta4
1812 !c
1813 !c gdelta4=gindex*delta4
1814 !c g1sdelt4=gindex*(1.-delta4)
1815 !c
1816 gdelta4=gindex*delta4
1817 g1sdelt4=gindex*(1.-delta4)
1818 !c
1819 !c combined water vapor depletions
1820 !c
1821 !cc graupel
1822 tmp=psdep(k)+pgdep(k)*gindex
1823 if ( tmp .gt. qvzodt(k) ) then
1824 factor=qvzodt(k)/tmp
1825 psdep(k)=psdep(k)*factor
1826 pgdep(k)=pgdep(k)*factor*gindex
1827 end if
1828 !c
1829 !c combined cloud water depletions
1830 !c
1831 tmp=praut(k)+psacw(k)+psfw(k)+pracw(k)+gindex*Pgacw(k)
1832 if ( tmp .gt. qlzodt(k) ) then
1833 factor=qlzodt(k)/tmp
1834 praut(k)=praut(k)*factor
1835 psacw(k)=psacw(k)*factor
1836 psfw(k)=psfw(k)*factor
1837 pracw(k)=pracw(k)*factor
1838 !cc graupel
1839 Pgacw(k)=Pgacw(k)*factor*gindex
1840 end if
1841 !c
1842 !c combined cloud ice depletions
1843 !c
1844 tmp=psaut(k)+psaci(k)+praci(k)+psfi(k)+Pgaci(k)*gdelta4 &
1845 +Pgacip(k)*g1sdelt4
1846 if (tmp .gt. qizodt(k) ) then
1847 factor=qizodt(k)/tmp
1848 psaut(k)=psaut(k)*factor
1849 psaci(k)=psaci(k)*factor
1850 praci(k)=praci(k)*factor
1851 psfi(k)=psfi(k)*factor
1852 !cc graupel
1853 Pgaci(k)=Pgaci(k)*factor*gdelta4
1854 Pgacip(k)=Pgacip(k)*factor*g1sdelt4
1855 endif
1856 !c
1857 !c combined all rain processes
1858 !c
1859 tmp_r=piacr(k)+psacr(k)-prevp(k)-praut(k)-pracw(k) &
1860 +Pgfr(k)*gindex+Pgacr(k)*gdelta4 &
1861 +Pgacrp(k)*g1sdelt4
1862 if (tmp_r .gt. qrzodt(k) ) then
1863 factor=qrzodt(k)/tmp_r
1864 piacr(k)=piacr(k)*factor
1865 psacr(k)=psacr(k)*factor
1866 prevp(k)=prevp(k)*factor
1867 !cc graupel
1868 Pgfr(k)=Pgfr(k)*factor*gindex
1869 Pgacr(k)=Pgacr(k)*factor*gdelta4
1870 Pgacrp(k)=Pgacrp(k)*factor*g1sdelt4
1871 endif
1873 !c
1874 !c if qrz < 1.0E-4 and qsz < 1.0E-4 then delta2=1.
1875 !c (all Pracs and Psacr become to snow)
1876 !c if qrz >= 1.0E-4 or qsz >= 1.0E-4 then delta2=0.
1877 !c (all Pracs and Psacr become to graupel)
1878 !c
1879 if (qrz(k) .lt. 1.0E-4 .and. qsz(k) .lt. 1.0E-4) then
1880 delta2=1.0
1881 else
1882 delta2=0.0
1883 endif
1884 !
1885 !cc graupel
1887 !c
1888 !c if qrz(k) < 1.0e-4 then delta3=1. means praci(k) --> qs
1889 !c piacr(k) --> qs
1890 !c if qrz(k) > 1.0e-4 then delta3=0. means praci(k) --> qg
1891 !c piacr(k) --> qg : Lin (20)
1893 if (qrz(k) .lt. 1.0e-4) then
1894 delta3=1.0
1895 else
1896 delta3=0.0
1897 endif
1898 !
1899 !c
1900 !c if gindex = 0.(no graupel) then delta2=1.0
1901 !c delta3=1.0
1902 !c
1903 if (gindex .eq. 0.) then
1904 delta2=1.0
1905 delta3=1.0
1906 endif
1907 !
1908 !c
1909 !c combined all snow processes
1910 !c
1911 tmp_s=-pssub(k)-(psaut(k)+psaci(k)+psacw(k)+psfw(k)+ &
1912 psfi(k)+praci(k)*delta3+piacr(k)*delta3+ &
1913 psdep(k))+Pgaut(k)*gindex+Pgacs(k)*gdelta4+ &
1914 Pgacsp(k)*g1sdelt4+Pracs(k)*(1.-delta2)- &
1915 Psacr(k)*delta2
1916 if ( tmp_s .gt. qszodt(k) ) then
1917 factor=qszodt(k)/tmp_s
1918 pssub(k)=pssub(k)*factor
1919 Pracs(k)=Pracs(k)*factor
1920 !cc graupel
1921 Pgaut(k)=Pgaut(k)*factor*gindex
1922 Pgacs(k)=Pgacs(k)*factor*gdelta4
1923 Pgacsp(k)=Pgacsp(k)*factor*g1sdelt4
1924 endif
1926 !cc graupel
1927 !
1929 ! if (gindex .eq. 0.) goto 998
1930 !c
1931 !c combined all graupel processes
1932 !c
1933 tmp_g=-pgaut(k)-pgfr(k)-Pgacw(k)*delta4-Pgaci(k)*delta4 &
1934 -Pgacr(k)*delta4-Pgacs(k)*delta4 &
1935 -pgwet(k)*(1.-delta4)-pgsub(k)-pgdep(k) &
1936 -psacr(k)*(1-delta2)-Pracs(k)*(1-delta2) &
1937 -praci(k)*(1-delta3)-piacr(k)*(1-delta3)
1938 if (tmp_g .gt. qgzodt(k)) then
1939 factor=qgzodt(k)/tmp_g
1940 pgsub(k)=pgsub(k)*factor
1941 endif
1943 998 continue
1944 !c
1945 !c calculate new water substances, thetae, tem, and qvsbar
1946 !c
1948 !cc graupel
1949 pvapor(k)=-pssub(k)-psdep(k)-prevp(k)-pgsub(k)*gindex &
1950 -pgdep(k)*gindex
1951 qvz(k)=amax1( qvmin,qvz(k)+dtb*pvapor(k) )
1952 pclw(k)=-praut(k)-pracw(k)-psacw(k)-psfw(k)-pgacw(k)*gindex
1953 if(flag_qndrop)then
1954 if( qlz(k) > 1e-20 ) &
1955 qndropz(k)=amax1( 0.0,qndropz(k)+dtb*pclw(k)*qndropz(k)/qlz(k) ) !sg
1956 endif
1957 qlz(k)=amax1( 0.0,qlz(k)+dtb*pclw(k) )
1958 pcli(k)=-psaut(k)-psfi(k)-psaci(k)-praci(k)-pgaci(k)*gdelta4 &
1959 -Pgacip(k)*g1sdelt4
1960 qiz(k)=amax1( 0.0,qiz(k)+dtb*pcli(k) )
1961 tmp_r=piacr(k)+psacr(k)-prevp(k)-praut(k)-pracw(k) &
1962 +Pgfr(k)*gindex+Pgacr(k)*gdelta4 &
1963 +Pgacrp(k)*g1sdelt4
1964 232 format(i2,1x,6(f9.3,1x))
1965 prain(k)=-tmp_r
1966 qrz(k)=amax1( 0.0,qrz(k)+dtb*prain(k) )
1967 tmp_s=-pssub(k)-(psaut(k)+psaci(k)+psacw(k)+psfw(k)+ &
1968 psfi(k)+praci(k)*delta3+piacr(k)*delta3+ &
1969 psdep(k))+Pgaut(k)*gindex+Pgacs(k)*gdelta4+ &
1970 Pgacsp(k)*g1sdelt4+Pracs(k)*(1.-delta2)- &
1971 Psacr(k)*delta2
1972 psnow(k)=-tmp_s
1973 qsz(k)=amax1( 0.0,qsz(k)+dtb*psnow(k) )
1974 qschg(k)=qschg(k)+psnow(k)
1975 qschg(k)=psnow(k)
1976 !cc graupel
1977 tmp_g=-pgaut(k)-pgfr(k)-Pgacw(k)*delta4-Pgaci(k)*delta4 &
1978 -Pgacr(k)*delta4-Pgacs(k)*delta4 &
1979 -pgwet(k)*(1.-delta4)-pgsub(k)-pgdep(k) &
1980 -psacr(k)*(1-delta2)-Pracs(k)*(1-delta2) &
1981 -praci(k)*(1-delta3)-piacr(k)*(1-delta3)
1982 252 format(i2,1x,6(f12.9,1x))
1983 262 format(i2,1x,7(f12.9,1x))
1984 pgraupel(k)=-tmp_g
1985 pgraupel(k)=pgraupel(k)*gindex
1986 qgz(k)=amax1( 0.0,qgz(k)+dtb*pgraupel(k))
1987 ! qgchg(k)=qgchg(k)+pgraupel(k)
1988 qgchg(k)=pgraupel(k)
1989 qgz(k)=qgz(k)*gindex
1991 tmp=ocp/tothz(k)*xLf*(qschg(k)+qgchg(k))
1992 theiz(k)=theiz(k)+dtb*tmp
1993 thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
1994 tem(k)=thz(k)*tothz(k)
1996 temcc(k)=tem(k)-273.15
1998 if( temcc(k) .lt. -40.0 ) qswz(k)=qsiz(k)
1999 qlpqi=qlz(k)+qiz(k)
2000 if ( qlpqi .eq. 0.0 ) then
2001 qvsbar(k)=qsiz(k)
2002 else
2003 qvsbar(k)=( qiz(k)*qsiz(k)+qlz(k)*qswz(k) )/qlpqi
2004 endif
2006 !
2007 else
2008 !c
2009 !c combined cloud water depletions
2010 !c
2011 tmp=praut(k)+psacw(k)+pracw(k)+pgacw(k)*gindex
2012 if ( tmp .gt. qlzodt(k) ) then
2013 factor=qlzodt(k)/tmp
2014 praut(k)=praut(k)*factor
2015 psacw(k)=psacw(k)*factor
2016 pracw(k)=pracw(k)*factor
2017 !cc graupel
2018 pgacw(k)=pgacw(k)*factor*gindex
2019 end if
2020 !c
2021 !c combined all snow processes
2022 !c
2023 tmp_s=-(psmlt(k)+psmltevp(k))+Pgacs(k)*gindex
2024 if (tmp_s .gt. qszodt(k) ) then
2025 factor=qszodt(k)/tmp_s
2026 psmlt(k)=psmlt(k)*factor
2027 psmltevp(k)=psmltevp(k)*factor
2028 !cc graupel
2029 Pgacs(k)=Pgacs(k)*factor*gindex
2030 endif
2032 !c
2033 !c
2034 !cc graupel
2035 !c
2036 ! if (gindex .eq. 0.) goto 997
2038 !c
2039 !c combined all graupel processes
2040 !c
2041 tmp_g=-pgmlt(k)-pgacs(k)-pgmltevp(k)
2042 if (tmp_g .gt. qgzodt(k)) then
2043 factor=qgzodt(k)/tmp_g
2044 pgmltevp(k)=pgmltevp(k)*factor
2045 pgmlt(k)=pgmlt(k)*factor
2046 endif
2047 !c
2048 997 continue
2050 !c
2051 !c combined all rain processes
2052 !c
2053 tmp_r=-prevp(k)-(praut(k)+pracw(k)+psacw(k)-psmlt(k)) &
2054 +pgmlt(k)*gindex-pgacw(k)*gindex
2055 if (tmp_r .gt. qrzodt(k) ) then
2056 factor=qrzodt(k)/tmp_r
2057 prevp(k)=prevp(k)*factor
2058 endif
2059 !c
2060 !c
2061 !c calculate new water substances and thetae
2062 !c
2065 pvapor(k)=-psmltevp(k)-prevp(k)-pgmltevp(k)*gindex
2066 qvz(k)=amax1( qvmin,qvz(k)+dtb*pvapor(k))
2067 pclw(k)=-praut(k)-pracw(k)-psacw(k)-pgacw(k)*gindex
2068 if(flag_qndrop)then
2069 if( qlz(k) > 1e-20 ) &
2070 qndropz(k)=amax1( 0.0,qndropz(k)+dtb*pclw(k)*qndropz(k)/qlz(k) ) !sg
2071 endif
2072 qlz(k)=amax1( 0.0,qlz(k)+dtb*pclw(k) )
2073 pcli(k)=0.0
2074 qiz(k)=amax1( 0.0,qiz(k)+dtb*pcli(k) )
2075 tmp_r=-prevp(k)-(praut(k)+pracw(k)+psacw(k)-psmlt(k)) &
2076 +pgmlt(k)*gindex-pgacw(k)*gindex
2077 242 format(i2,1x,7(f9.6,1x))
2078 prain(k)=-tmp_r
2079 tmpqrz=qrz(k)
2080 qrz(k)=amax1( 0.0,qrz(k)+dtb*prain(k) )
2081 tmp_s=-(psmlt(k)+psmltevp(k))+Pgacs(k)*gindex
2082 psnow(k)=-tmp_s
2083 qsz(k)=amax1( 0.0,qsz(k)+dtb*psnow(k) )
2084 ! qschg(k)=qschg(k)+psnow(k)
2085 qschg(k)=psnow(k)
2086 !cc graupel
2088 tmp_g=-pgmlt(k)-pgacs(k)-pgmltevp(k)
2089 ! write(*,272)k,pgmlt(k),pgacs(k),pgmltevp(k),
2090 272 format(i2,1x,3(f12.9,1x))
2091 pgraupel(k)=-tmp_g*gindex
2092 qgz(k)=amax1( 0.0,qgz(k)+dtb*pgraupel(k))
2093 ! qgchg(k)=qgchg(k)+pgraupel(k)
2094 qgchg(k)=pgraupel(k)
2095 qgz(k)=qgz(k)*gindex
2096 !
2097 tmp=ocp/tothz(k)*xLf*(qschg(k)+qgchg(k))
2098 theiz(k)=theiz(k)+dtb*tmp
2099 thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
2101 tem(k)=thz(k)*tothz(k)
2102 temcc(k)=tem(k)-273.15
2103 ! qswz(k)=episp0k*oprez(k)* &
2104 ! exp( svp2*temcc(k)/(tem(k)-svp3) )
2105 es=1000.*svp1*exp( svp2*temcc(k)/(tem(k)-svp3) )
2106 qswz(k)=ep2*es/(prez(k)-es)
2107 qsiz(k)=qswz(k)
2108 qvsbar(k)=qswz(k)
2109 !
2110 end if
2111 preclw(k)=pclw(k) !sg
2113 !
2114 !***********************************************************************
2115 !********** saturation adjustment **********
2116 !***********************************************************************
2117 !
2118 ! allow supersaturation exits linearly from 0% at 500 mb to 50%
2119 ! above 300 mb
2120 ! 5.0e-5=1.0/(500mb-300mb)
2121 !
2122 rsat=1.0+0.5*(50000.0-prez(k))*5.0e-5
2123 rsat=amax1(1.0,rsat)
2124 rsat=amin1(1.5,rsat)
2125 rsat=1.0
2126 if( qvz(k)+qlz(k)+qiz(k) .lt. rsat*qvsbar(k) ) then
2128 !c
2129 !c unsaturated
2130 !c
2131 qvz(k)=qvz(k)+qlz(k)+qiz(k)
2132 qlz(k)=0.0
2133 qiz(k)=0.0
2135 thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
2136 tem(k)=thz(k)*tothz(k)
2137 temcc(k)=tem(k)-273.15
2139 go to 1800
2140 !
2141 else
2142 !c
2143 !c saturated
2144 !c
2145 !
2146 pladj(k)=qlz(k)
2147 piadj(k)=qiz(k)
2148 !
2150 CALL satadj(qvz, qlz, qiz, prez, theiz, thz, tothz, kts, kte, &
2151 k, xLvocp, xLfocp, episp0k, EP2,SVP1,SVP2,SVP3,SVPT0 )
2153 !
2154 pladj(k)=odtb*(qlz(k)-pladj(k))
2155 piadj(k)=odtb*(qiz(k)-piadj(k))
2156 !
2157 pclw(k)=pclw(k)+pladj(k)
2158 pcli(k)=pcli(k)+piadj(k)
2159 pvapor(k)=pvapor(k)-( pladj(k)+piadj(k) )
2160 !
2161 thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
2162 tem(k)=thz(k)*tothz(k)
2164 temcc(k)=tem(k)-273.15
2166 ! qswz(k)=episp0k*oprez(k)* &
2167 ! exp( svp2*temcc(k)/(tem(k)-svp3) )
2168 es=1000.*svp1*exp( svp2*temcc(k)/(tem(k)-svp3) )
2169 qswz(k)=ep2*es/(prez(k)-es)
2170 if (tem(k) .lt. 273.15 ) then
2171 ! qsiz(k)=episp0k*oprez(k)* &
2172 ! exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
2173 es=1000.*svp1*exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
2174 qsiz(k)=ep2*es/(prez(k)-es)
2175 if (temcc(k) .lt. -40.0) qswz(k)=qsiz(k)
2176 else
2177 qsiz(k)=qswz(k)
2178 endif
2179 qlpqi=qlz(k)+qiz(k)
2180 if ( qlpqi .eq. 0.0 ) then
2181 qvsbar(k)=qsiz(k)
2182 else
2183 qvsbar(k)=( qiz(k)*qsiz(k)+qlz(k)*qswz(k) )/qlpqi
2184 endif
2186 end if
2188 !
2189 !***********************************************************************
2190 !***** melting and freezing of cloud ice and cloud water *****
2191 !***********************************************************************
2192 qlpqi=qlz(k)+qiz(k)
2193 if(qlpqi .le. 0.0) go to 1800
2194 !
2195 !c
2196 !c (1) HOMOGENEOUS NUCLEATION WHEN T< -40 C (Pihom)
2197 !c
2198 if(temcc(k) .lt. -40.0) pihom(k)=qlz(k)*odtb
2199 !c
2200 !c (2) MELTING OF ICE CRYSTAL WHEN T> 0 C (Pimlt)
2201 !c
2202 if(temcc(k) .gt. 0.0) pimlt(k)=qiz(k)*odtb
2203 !c
2204 !c (3) PRODUCTION OF CLOUD ICE BY BERGERON PROCESS (Pidw): Hsie (p957)
2205 !c this process only considered when -31 C < T < 0 C
2206 !c
2207 if(temcc(k) .lt. 0.0 .and. temcc(k) .gt. -31.0) then
2208 !c!
2209 !c! parama1 and parama2 functions must be user supplied
2210 !c!
2211 a1=parama1( temcc(k) )
2212 a2=parama2( temcc(k) )
2213 !! change unit from cgs to mks
2214 a1=a1*0.001**(1.0-a2)
2215 xnin=xni0*exp(-bni*temcc(k))
2216 pidw(k)=xnin*orho(k)*(a1*xmnin**a2)
2217 end if
2218 !
2219 pcli(k)=pcli(k)+pihom(k)-pimlt(k)+pidw(k)
2220 pclw(k)=pclw(k)-pihom(k)+pimlt(k)-pidw(k)
2221 qlz(k)=amax1( 0.0,qlz(k)+dtb*(-pihom(k)+pimlt(k)-pidw(k)) )
2222 qiz(k)=amax1( 0.0,qiz(k)+dtb*(pihom(k)-pimlt(k)+pidw(k)) )
2224 !
2225 CALL satadj(qvz, qlz, qiz, prez, theiz, thz, tothz, kts, kte, &
2226 k, xLvocp, xLfocp, episp0k ,EP2,SVP1,SVP2,SVP3,SVPT0)
2228 thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
2229 tem(k)=thz(k)*tothz(k)
2231 temcc(k)=tem(k)-273.15
2233 ! qswz(k)=episp0k*oprez(k)* &
2234 ! exp( svp2*temcc(k)/(tem(k)-svp3) )
2235 es=1000.*svp1*exp( svp2*temcc(k)/(tem(k)-svp3) )
2236 qswz(k)=ep2*es/(prez(k)-es)
2238 if (tem(k) .lt. 273.15 ) then
2239 ! qsiz(k)=episp0k*oprez(k)* &
2240 ! exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
2241 es=1000.*svp1*exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
2242 qsiz(k)=ep2*es/(prez(k)-es)
2243 if (temcc(k) .lt. -40.0) qswz(k)=qsiz(k)
2244 else
2245 qsiz(k)=qswz(k)
2246 endif
2247 qlpqi=qlz(k)+qiz(k)
2248 if ( qlpqi .eq. 0.0 ) then
2249 qvsbar(k)=qsiz(k)
2250 else
2251 qvsbar(k)=( qiz(k)*qsiz(k)+qlz(k)*qswz(k) )/qlpqi
2252 endif
2254 1800 continue
2255 !
2256 !***********************************************************************
2257 !********** integrate the productions of rain and snow **********
2258 !***********************************************************************
2259 !c
2261 2000 continue
2264 !---------------------------------------------------------------------
2266 !
2267 !***********************************************************************
2268 !****** Write terms in cloud physics to time series dataset *****
2269 !***********************************************************************
2270 !
2271 ! open(unit=24,form='formatted',status='new',
2272 ! & file='cloud.dat')
2274 !9030 format(10e12.6)
2276 ! write(24,*)'tmp'
2277 ! write(24,9030) (tem(k),k=kts+1,kte)
2278 ! write(24,*)'qiz'
2279 ! write(24,9030) (qiz(k),k=kts+1,kte)
2280 ! write(24,*)'qsz'
2281 ! write(24,9030) (qsz(k),k=kts+1,kte)
2282 ! write(24,*)'qrz'
2283 ! write(24,9030) (qrz(k),k=kts+1,kte)
2284 ! write(24,*)'qgz'
2285 ! write(24,9030) (qgz(k),k=kts+1,kte)
2286 ! write(24,*)'qvoqsw'
2287 ! write(24,9030) (qvoqswz(k),k=kts+1,kte)
2288 ! write(24,*)'qvoqsi'
2289 ! write(24,9030) (qvoqsiz(k),k=kts+1,kte)
2290 ! write(24,*)'vtr'
2291 ! write(24,9030) (vtr(k),k=kts+1,kte)
2292 ! write(24,*)'vts'
2293 ! write(24,9030) (vts(k),k=kts+1,kte)
2294 ! write(24,*)'vtg'
2295 ! write(24,9030) (vtg(k),k=kts+1,kte)
2296 ! write(24,*)'pclw'
2297 ! write(24,9030) (pclw(k),k=kts+1,kte)
2298 ! write(24,*)'pvapor'
2299 ! write(24,9030) (pvapor(k),k=kts+1,kte)
2300 ! write(24,*)'pcli'
2301 ! write(24,9030) (pcli(k),k=kts+1,kte)
2302 ! write(24,*)'pimlt'
2303 ! write(24,9030) (pimlt(k),k=kts+1,kte)
2304 ! write(24,*)'pihom'
2305 ! write(24,9030) (pihom(k),k=kts+1,kte)
2306 ! write(24,*)'pidw'
2307 ! write(24,9030) (pidw(k),k=kts+1,kte)
2308 ! write(24,*)'prain'
2309 ! write(24,9030) (prain(k),k=kts+1,kte)
2310 ! write(24,*)'praut'
2311 ! write(24,9030) (praut(k),k=kts+1,kte)
2312 ! write(24,*)'pracw'
2313 ! write(24,9030) (pracw(k),k=kts+1,kte)
2314 ! write(24,*)'prevp'
2315 ! write(24,9030) (prevp(k),k=kts+1,kte)
2316 ! write(24,*)'psnow'
2317 ! write(24,9030) (psnow(k),k=kts+1,kte)
2318 ! write(24,*)'psaut'
2319 ! write(24,9030) (psaut(k),k=kts+1,kte)
2320 ! write(24,*)'psfw'
2321 ! write(24,9030) (psfw(k),k=kts+1,kte)
2322 ! write(24,*)'psfi'
2323 ! write(24,9030) (psfi(k),k=kts+1,kte)
2324 ! write(24,*)'praci'
2325 ! write(24,9030) (praci(k),k=kts+1,kte)
2326 ! write(24,*)'piacr'
2327 ! write(24,9030) (piacr(k),k=kts+1,kte)
2328 ! write(24,*)'psaci'
2329 ! write(24,9030) (psaci(k),k=kts+1,kte)
2330 ! write(24,*)'psacw'
2331 ! write(24,9030) (psacw(k),k=kts+1,kte)
2332 ! write(24,*)'psdep'
2333 ! write(24,9030) (psdep(k),k=kts+1,kte)
2334 ! write(24,*)'pssub'
2335 ! write(24,9030) (pssub(k),k=kts+1,kte)
2336 ! write(24,*)'pracs'
2337 ! write(24,9030) (pracs(k),k=kts+1,kte)
2338 ! write(24,*)'psacr'
2339 ! write(24,9030) (psacr(k),k=kts+1,kte)
2340 ! write(24,*)'psmlt'
2341 ! write(24,9030) (psmlt(k),k=kts+1,kte)
2342 ! write(24,*)'psmltevp'
2343 ! write(24,9030) (psmltevp(k),k=kts+1,kte)
2344 ! write(24,*)'pladj'
2345 ! write(24,9030) (pladj(k),k=kts+1,kte)
2346 ! write(24,*)'piadj'
2347 ! write(24,9030) (piadj(k),k=kts+1,kte)
2348 ! write(24,*)'pgraupel'
2349 ! write(24,9030) (pgraupel(k),k=kts+1,kte)
2350 ! write(24,*)'pgaut'
2351 ! write(24,9030) (pgaut(k),k=kts+1,kte)
2352 ! write(24,*)'pgfr'
2353 ! write(24,9030) (pgfr(k),k=kts+1,kte)
2354 ! write(24,*)'pgacw'
2355 ! write(24,9030) (pgacw(k),k=kts+1,kte)
2356 ! write(24,*)'pgaci'
2357 ! write(24,9030) (pgaci(k),k=kts+1,kte)
2358 ! write(24,*)'pgacr'
2359 ! write(24,9030) (pgacr(k),k=kts+1,kte)
2360 ! write(24,*)'pgacs'
2361 ! write(24,9030) (pgacs(k),k=kts+1,kte)
2362 ! write(24,*)'pgacip'
2363 ! write(24,9030) (pgacip(k),k=kts+1,kte)
2364 ! write(24,*)'pgacrP'
2365 ! write(24,9030) (pgacrP(k),k=kts+1,kte)
2366 ! write(24,*)'pgacsp'
2367 ! write(24,9030) (pgacsp(k),k=kts+1,kte)
2368 ! write(24,*)'pgwet'
2369 ! write(24,9030) (pgwet(k),k=kts+1,kte)
2370 ! write(24,*)'pdry'
2371 ! write(24,9030) (pdry(k),k=kts+1,kte)
2372 ! write(24,*)'pgsub'
2373 ! write(24,9030) (pgsub(k),k=kts+1,kte)
2374 ! write(24,*)'pgdep'
2375 ! write(24,9030) (pgdep(k),k=kts+1,kte)
2376 ! write(24,*)'pgmlt'
2377 ! write(24,9030) (pgmlt(k),k=kts+1,kte)
2378 ! write(24,*)'pgmltevp'
2379 ! write(24,9030) (pgmltevp(k),k=kts+1,kte)
2383 !**** below if qv < qvmin then qv=qvmin, ql=0.0, and qi=0.0
2384 !
2385 do k=kts+1,kte
2386 if ( qvz(k) .lt. qvmin ) then
2387 qlz(k)=0.0
2388 qiz(k)=0.0
2389 qvz(k)=amax1( qvmin,qvz(k)+qlz(k)+qiz(k) )
2390 end if
2391 enddo
2392 !
2393 END SUBROUTINE clphy1d
2396 !---------------------------------------------------------------------
2397 ! SATURATED ADJUSTMENT
2398 !---------------------------------------------------------------------
2399 SUBROUTINE satadj(qvz, qlz, qiz, prez, theiz, thz, tothz, &
2400 kts, kte, k, xLvocp, xLfocp, episp0k, EP2,SVP1,SVP2,SVP3,SVPT0)
2401 !---------------------------------------------------------------------
2402 IMPLICIT NONE
2403 !---------------------------------------------------------------------
2404 ! This program use Newton's method for finding saturated temperature
2405 ! and saturation mixing ratio.
2406 !
2407 ! In this saturation adjustment scheme we assume
2408 ! (1) the saturation mixing ratio is the mass weighted average of
2409 ! saturation values over liquid water (qsw), and ice (qsi)
2410 ! following Lord et al., 1984 and Tao, 1989
2411 !
2412 ! (2) the percentage of cloud liquid and cloud ice will
2413 ! be fixed during the saturation calculation
2414 !---------------------------------------------------------------------
2415 !
2417 INTEGER, INTENT(IN ) :: kts, kte, k
2419 REAL, DIMENSION( kts:kte ), &
2420 INTENT(INOUT) :: qvz, qlz, qiz
2421 !
2422 REAL, DIMENSION( kts:kte ), &
2423 INTENT(IN ) :: prez, theiz, tothz
2425 REAL, INTENT(IN ) :: xLvocp, xLfocp, episp0k
2426 REAL, INTENT(IN ) :: EP2,SVP1,SVP2,SVP3,SVPT0
2428 ! LOCAL VARS
2430 INTEGER :: n
2432 REAL, DIMENSION( kts:kte ) :: thz, tem, temcc, qsiz, &
2433 qswz, qvsbar
2435 REAL :: qsat, qlpqi, ratql, t0, t1, tmp1, ratqi, tsat, absft, &
2436 denom1, denom2, dqvsbar, ftsat, dftsat, qpz, &
2437 gindex, es
2438 !
2439 !---------------------------------------------------------------------
2441 thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
2443 tem(k)=tothz(k)*thz(k)
2444 if (tem(k) .gt. 273.15) then
2445 ! qsat=episp0k/prez(k)* &
2446 ! exp( svp2*(tem(k)-273.15)/(tem(k)-svp3) )
2447 es=1000.*svp1*exp( svp2*(tem(k)-svpt0)/(tem(k)-svp3) )
2448 qsat=ep2*es/(prez(k)-es)
2449 else
2450 qsat=episp0k/prez(k)* &
2451 exp( 21.8745584*(tem(k)-273.15)/(tem(k)-7.66) )
2452 end if
2453 qpz=qvz(k)+qlz(k)+qiz(k)
2454 if (qpz .lt. qsat) then
2455 qvz(k)=qpz
2456 qiz(k)=0.0
2457 qlz(k)=0.0
2458 go to 400
2459 end if
2460 qlpqi=qlz(k)+qiz(k)
2461 if( qlpqi .ge. 1.0e-5) then
2462 ratql=qlz(k)/qlpqi
2463 ratqi=qiz(k)/qlpqi
2464 else
2465 t0=273.15
2466 ! t1=233.15
2467 t1=248.15
2468 tmp1=( t0-tem(k) )/(t0-t1)
2469 tmp1=amin1(1.0,tmp1)
2470 tmp1=amax1(0.0,tmp1)
2471 ratqi=tmp1
2472 ratql=1.0-tmp1
2473 end if
2474 !
2475 !
2476 !-- saturation mixing ratios over water and ice
2477 !-- at the outset we will follow Bolton 1980 MWR for
2478 !-- the water and Murray JAS 1967 for the ice
2479 !
2480 !-- dqvsbar=d(qvsbar)/dT
2481 !-- ftsat=F(Tsat)
2482 !-- dftsat=d(F(T))/dT
2483 !
2484 ! First guess of tsat
2486 tsat=tem(k)
2487 absft=1.0
2488 !
2489 do 200 n=1,20
2490 denom1=1.0/(tsat-svp3)
2491 denom2=1.0/(tsat-7.66)
2492 ! qswz(k)=episp0k/prez(k)* &
2493 ! exp( svp2*denom1*(tsat-273.15) )
2494 es=1000.*svp1*exp( svp2*denom1*(tsat-svpt0) )
2495 qswz(k)=ep2*es/(prez(k)-es)
2496 if (tem(k) .lt. 273.15) then
2497 ! qsiz(k)=episp0k/prez(k)* &
2498 ! exp( 21.8745584*denom2*(tsat-273.15) )
2499 es=1000.*svp1*exp( 21.8745584*denom2*(tsat-273.15) )
2500 qsiz(k)=ep2*es/(prez(k)-es)
2501 if (tem(k) .lt. 233.15) qswz(k)=qsiz(k)
2502 else
2503 qsiz(k)=qswz(k)
2504 endif
2505 qvsbar(k)=ratql*qswz(k)+ratqi*qsiz(k)
2506 !
2507 ! if( absft .lt. 0.01 .and. n .gt. 3 ) go to 300
2508 if( absft .lt. 0.01 ) go to 300
2509 !
2510 dqvsbar=ratql*qswz(k)*svp2*243.5*denom1*denom1+ &
2511 ratqi*qsiz(k)*21.8745584*265.5*denom2*denom2
2512 ftsat=tsat+(xlvocp+ratqi*xlfocp)*qvsbar(k)- &
2513 tothz(k)*theiz(k)-xlfocp*ratqi*(qvz(k)+qlz(k)+qiz(k))
2514 dftsat=1.0+(xlvocp+ratqi*xlfocp)*dqvsbar
2515 tsat=tsat-ftsat/dftsat
2516 absft=abs(ftsat)
2518 200 continue
2519 9020 format(1x,'point can not converge, absft,n=',e12.5,i5)
2520 !
2521 300 continue
2522 if( qpz .gt. qvsbar(k) ) then
2523 qvz(k)=qvsbar(k)
2524 qiz(k)=ratqi*( qpz-qvz(k) )
2525 qlz(k)=ratql*( qpz-qvz(k) )
2526 else
2527 qvz(k)=qpz
2528 qiz(k)=0.0
2529 qlz(k)=0.0
2530 end if
2531 400 continue
2533 END SUBROUTINE satadj
2536 !----------------------------------------------------------------
2537 REAL FUNCTION parama1(temp)
2538 !----------------------------------------------------------------
2539 IMPLICIT NONE
2540 !----------------------------------------------------------------
2541 ! This program calculate the parameter for crystal growth rate
2542 ! in Bergeron process
2543 !----------------------------------------------------------------
2545 REAL, INTENT (IN ) :: temp
2546 REAL, DIMENSION(32) :: a1
2547 INTEGER :: i1, i1p1
2548 REAL :: ratio
2550 data a1/0.100e-10,0.7939e-7,0.7841e-6,0.3369e-5,0.4336e-5, &
2551 0.5285e-5,0.3728e-5,0.1852e-5,0.2991e-6,0.4248e-6, &
2552 0.7434e-6,0.1812e-5,0.4394e-5,0.9145e-5,0.1725e-4, &
2553 0.3348e-4,0.1725e-4,0.9175e-5,0.4412e-5,0.2252e-5, &
2554 0.9115e-6,0.4876e-6,0.3473e-6,0.4758e-6,0.6306e-6, &
2555 0.8573e-6,0.7868e-6,0.7192e-6,0.6513e-6,0.5956e-6, &
2556 0.5333e-6,0.4834e-6/
2558 i1=int(-temp)+1
2559 i1p1=i1+1
2560 ratio=-(temp)-float(i1-1)
2561 parama1=a1(i1)+ratio*( a1(i1p1)-a1(i1) )
2563 END FUNCTION parama1
2565 !----------------------------------------------------------------
2566 REAL FUNCTION parama2(temp)
2567 !----------------------------------------------------------------
2568 IMPLICIT NONE
2569 !----------------------------------------------------------------
2570 ! This program calculate the parameter for crystal growth rate
2571 ! in Bergeron process
2572 !----------------------------------------------------------------
2574 REAL, INTENT (IN ) :: temp
2575 REAL, DIMENSION(32) :: a2
2576 INTEGER :: i1, i1p1
2577 REAL :: ratio
2579 data a2/0.0100,0.4006,0.4831,0.5320,0.5307,0.5319,0.5249, &
2580 0.4888,0.3849,0.4047,0.4318,0.4771,0.5183,0.5463, &
2581 0.5651,0.5813,0.5655,0.5478,0.5203,0.4906,0.4447, &
2582 0.4126,0.3960,0.4149,0.4320,0.4506,0.4483,0.4460, &
2583 0.4433,0.4413,0.4382,0.4361/
2584 i1=int(-temp)+1
2585 i1p1=i1+1
2586 ratio=-(temp)-float(i1-1)
2587 parama2=a2(i1)+ratio*( a2(i1p1)-a2(i1) )
2589 END FUNCTION parama2
2591 !----------------------------------------------------------------
2592 REAL FUNCTION ggamma(X)
2593 !----------------------------------------------------------------
2594 IMPLICIT NONE
2595 !----------------------------------------------------------------
2596 REAL, INTENT(IN ) :: x
2597 REAL, DIMENSION(8) :: B
2598 INTEGER ::j, K1
2599 REAL ::PF, G1TO2 ,TEMP
2601 DATA B/-.577191652,.988205891,-.897056937,.918206857, &
2602 -.756704078,.482199394,-.193527818,.035868343/
2604 PF=1.
2605 TEMP=X
2606 DO 10 J=1,200
2607 IF (TEMP .LE. 2) GO TO 20
2608 TEMP=TEMP-1.
2609 10 PF=PF*TEMP
2610 100 FORMAT(//,5X,'module_mp_lin: INPUT TO GAMMA FUNCTION TOO LARGE, X=',E12.5)
2611 WRITE(wrf_err_message,100)X
2612 CALL wrf_error_fatal(wrf_err_message)
2613 20 G1TO2=1.
2614 TEMP=TEMP - 1.
2615 DO 30 K1=1,8
2616 30 G1TO2=G1TO2 + B(K1)*TEMP**K1
2617 ggamma=PF*G1TO2
2619 END FUNCTION ggamma
2621 !----------------------------------------------------------------
2623 END MODULE module_mp_lin