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