module_mp_lin.F
References to this file elsewhere.
1 !WRF:MODEL_LAYER:PHYSICS
2 !
3
4 MODULE module_mp_lin
5
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
12
13 ! Lin
14 ! REAL , PARAMETER, PRIVATE :: xnog = 4.0e4
15 ! REAL , PARAMETER, PRIVATE :: rhograul = 917.
16
17 ! Hobbs
18 REAL , PARAMETER, PRIVATE :: xnog = 4.0e6
19 REAL , PARAMETER, PRIVATE :: rhograul = 400.
20
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
33
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
65
66 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
67 INTENT(INOUT) :: &
68 th, &
69 qv, &
70 ql, &
71 qr
72
73 !
74 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
75 INTENT(IN ) :: &
76 rho, &
77 pii, &
78 p, &
79 dz8w
80
81 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
82 INTENT(IN ) :: z
83
84
85
86 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: ht
87
88 REAL, INTENT(IN ) :: dt_in, &
89 grav, &
90 Rair, &
91 rvapor, &
92 cp, &
93 XLS, &
94 XLV, &
95 XLF, &
96 rhowater, &
97 rhosnow
98
99 REAL, INTENT(IN ) :: EP2,SVP1,SVP2,SVP3,SVPT0
100
101 REAL, DIMENSION( ims:ime , jms:jme ), &
102 INTENT(INOUT) :: RAINNC, &
103 RAINNCV
104
105 ! Optional
106
107 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
108 OPTIONAL, &
109 INTENT(INOUT) :: &
110 qi, &
111 qs, &
112 qg, &
113 qndrop
114
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)
122
123 LOGICAL, INTENT(IN), OPTIONAL :: F_QG, F_QNDROP
124
125 ! LOCAL VAR
126
127 INTEGER :: min_q, max_q
128
129 REAL, DIMENSION( its:ite , jts:jte ) &
130 :: rain, snow, graupel,ice
131
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
141
142 LOGICAL :: flag_qg, flag_qndrop
143 !
144 REAL :: dt, pptrain, pptsnow, pptgraul, rhoe_s, &
145 gindex, pptice
146 real :: qndropconst
147
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
159
160 IF (.not.flag_qg) gindex=0.
161
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
181
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
191
192 DO k = kts, kte
193 qiz(k)=qi(i,k,j)
194 qsz(k)=qs(i,k,j)
195 ENDDO
196
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 )
221
222 !
223 ! Precipitation from cloud microphysics -- only for one time step
224 !
225 ! unit is transferred from m to mm
226
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
235
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
250
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
263
264 DO k = kts, kte
265 qi(i,k,j)=qiz(k)
266 qs(i,k,j)=qsz(k)
267 ENDDO
268
269 IF ( present(preci) ) THEN !sg beg
270 DO k = kts, kte
271 preci(i,k,j)=preciz(k)
272 ENDDO
273 END IF
274
275 IF ( present(precs) ) THEN
276 DO k = kts, kte
277 precs(i,k,j)=precsz(k)
278 ENDDO
279 END IF !sg end
280
281 IF ( flag_qg .AND. PRESENT( qg ) ) THEN
282 DO k = kts, kte
283 qg(i,k,j)=qgz(k)
284 ENDDO
285
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
297
298 END SUBROUTINE lin_et_al
299
300
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 !----------------------------------------------------------------------
368
369 INTEGER, INTENT(IN ) :: kts, kte, i, j
370
371 REAL, DIMENSION( kts:kte ), &
372 INTENT(INOUT) :: qvz, qlz, qrz, qiz, qsz, &
373 qndropz, &
374 qgz, thz
375
376 REAL, DIMENSION( kts:kte ), &
377 INTENT(IN ) :: tothz, rho, orho, sqrho, &
378 prez, zz, dzw
379
380 REAL, INTENT(IN ) :: dt, grav, cp, Rair, rvapor, &
381 XLS, XLV, XLF, rhowater, &
382 rhosnow,EP2,SVP1,SVP2,SVP3,SVPT0
383
384 REAL, DIMENSION( kts:kte ), INTENT(OUT) :: preclw, &
385 precrz, preciz, precsz, precgz
386
387 REAL, INTENT(INOUT) :: pptrain, pptsnow, pptgraul, pptice
388
389 REAL, INTENT(IN ) :: zsfc
390 logical, intent(in) :: flag_qndrop !sg
391
392 ! local vars
393
394 REAL :: obp4, bp3, bp5, bp6, odp4, &
395 dp3, dp5, dp5o2
396
397
398 ! temperary vars
399
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
411
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
418
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 !
430
431 REAL, DIMENSION( kts:kte ) :: qvsbar, rs0, viscmu, visc, diffwv, &
432 schmidt, xka
433
434 REAL, DIMENSION( kts:kte ) :: vtr, vts, vtg, &
435 vtrold, vtsold, vtgold, vtiold, &
436 xlambdar, xlambdas, xlambdag, &
437 olambdar, olambdas, olambdag
438
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
449
450 ! for terminal velocity flux
451
452 INTEGER :: min_q, max_q
453 REAL :: t_del_tv, del_tv, flux, fluxin, fluxout ,tmpqrz
454 LOGICAL :: notlast
455 !
456
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
478
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 !
527
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
540
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) )
575
576 theiz(k)=thz(k)+(xlvocp*qvz(k)-xlfocp*qiz(k))/tothz(k)
577 enddo
578
579
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 !
636
637
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
675
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
742
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)
747
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
774
775 !
776 !-- snow
777 !
778 t_del_tv=0.
779 del_tv=dtb
780 notlast=.true.
781
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
803
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)
809
810 t_del_tv=t_del_tv+del_tv
811
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
835
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
866
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)
872
873 t_del_tv=t_del_tv+del_tv
874
875 if ( t_del_tv .ge. dtb ) then
876 notlast=.false.
877 del_tv=dtb+del_tv-t_del_tv
878 endif
879
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
901
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
928
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)
934
935 t_del_tv=t_del_tv+del_tv
936
937 if ( t_del_tv .ge. dtb ) then
938 notlast=.false.
939 del_tv=dtb+del_tv-t_del_tv
940 endif
941
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
972
973
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
998
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------------------------------------------------------------------
1088
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)
1094
1095 if (tem(k) .lt. 273.15) then
1096
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) )
1107
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) )
1119
1120 tmp1=odtb*(qiz(k)-qic)*(1.0-exp(-alpha1*dtb))
1121 psaut(k)=amax1( 0.0,tmp1 )
1122
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!
1131
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.
1170
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
1177
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)
1181
1182 !c! number of 50 micron crystals produced
1183 !0915 xni50=qiz(k)*dtb*odtberg/xmi50
1184
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 !
1203
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) )
1222
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) )
1226
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) )
1233
1234 tmp2=qiz(k)*save1*rho(k)*pio6*rhowater*gambp6*oxmi* &
1235 olambdar(k)**bp6
1236 piacr(k)=amin1( tmp2,qrzodt(k) )
1237
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) )
1253
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) )
1264
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
1294
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
1334
1335 tmp1=esw*pio4*xnos*constc*gamdp3*sqrho(k)* &
1336 olambdas(k)**dp3
1337 psacw(k)=qlzodt(k)*( 1.0-exp(-tmp1*dtb) )
1338
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 )
1376
1377 ! abr=2.0*pi*(qvoqswz(k)-1.0)*tmpc/(tmpa+tmpb)
1378
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
1399
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
1438
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) )
1446
1447 tmp1=pio4*erw*xnor*consta*sqrho(k)* &
1448 gambp3*olambdar(k)**bp3
1449 pracw(k)=qlzodt(k)*( 1.0-exp(-tmp1*dtb) )
1450
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)
1463
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) )
1477
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)
1482
1483
1484
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 !
1502
1503 ! testing
1504 ! tmp1=alpha2*(qsz(k)-qs0)
1505 ! tmp1=amax1(0.0,tmp1)
1506 ! pgaut(k)=amin1( tmp1,qszodt(k) )
1507
1508 tmp1=odtb*(qsz(k)-qs0)*(1.0-exp(-alpha2*dtb))
1509 pgaut(k)=amax1( 0.0,tmp1 )
1510
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 !
1517
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
1528
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
1533
1534 !c
1535 !c Comparing Pgwet(wet process) and Pdry(dry process),
1536 !c we will pick up the small one.
1537 !c
1538
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) )
1561
1562
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) )
1575
1576
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) )
1590
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)
1595
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) )
1606
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) )
1614
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
1622
1623 term0=constg*olambdag(k)**5.5/visc(k)
1624
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)
1631
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) )
1642
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
1657
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)
1667
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 constg2=vf1s*olambdag(k)*olambdag(k)+ &
1687 vf2s*schmidt(k)**0.33334*gam2pt75*sqrt(term0)
1688 tmp2=abg*xnog*constg2
1689 pgdep(k)=amax1(0.0,tmp2)
1690 pgsub(k)=amin1(0.0,tmp2)
1691 pgsub(k)=amax1( pgsub(k),-qgzodt(k) )
1692
1693 4000 continue
1694 else
1695 !
1696 !***********************************************************************
1697 !********* graupel production processes for T > 0 C **********
1698 !***********************************************************************
1699 !
1700 !c
1701 !c (1) ACCRETION OF CLOUD WATER BY GRAUPEL (Pgacw): Lin (40)
1702 !c egw=1.0
1703 !c Cdrag=0.6 drag coefficients for hairstone
1704 !c constg=sqrt(4.*grav*rhograul*0.33334*orho(k)/Cdrag)
1705
1706 egw=1.0
1707 constg=sqrt(4.*grav*rhograul*0.33334*orho(k)*oCdrag)
1708 tmp1=pio4*xnog*gam3pt5*constg*olambdag(k)**3.5
1709 tmp2=qlz(k)*egw*tmp1
1710 Pgacw(k)=amin1( tmp2,qlzodt(k) )
1711
1712 !c
1713 !c (2) ACCRETION OF RAIN BY GRAUPEL (Pgacr): Lin (42)
1714 !c Compute processes (5) only when qrz > 0.0 and qgz > 0.0
1715 !c egr=1.
1716 !c
1717 egr=1.
1718 tmpa=olambdar(k)*olambdar(k)
1719 tmpb=olambdag(k)*olambdag(k)
1720 tmpc=olambdar(k)*olambdag(k)
1721 tmp1=pi*pi*xnor*xnog*abs( vtr(k)-vtg(k) )*orho(k)
1722 tmp2=tmpa*tmpa*olambdag(k)*(5.0*tmpa+2.0*tmpc+0.5*tmpb)
1723 tmp3=tmp1*egr*rhowater*tmp2
1724 pgacr(k)=amin1( tmp3,qrzodt(k) )
1725
1726
1727 !c
1728 !c (3) GRAUPEL MELTING TO FORM RAIN (Pgmlt): Lin (47)
1729 !c Pgmlt is negative value
1730 !c constg=sqrt(4.*grav*rhograul*0.33334*orho(k)/Cdrag)
1731 !c constg2=vf1s*olambdag(k)*olambdag(k)+
1732 !c vf2s*schmidt(k)**0.33334*gam2pt75*constg
1733 !c Cdrag=0.6 drag coefficients for hairstone
1734 !
1735 delrs=rs0(k)-qvz(k)
1736 term1=2.0*pi*orho(k)*( xlv*diffwv(k)*rho(k)*delrs- &
1737 xka(k)*temcc(k) )
1738 term0=sqrt(4.*grav*rhograul*0.33334*orho(k)*ocdrag) &
1739 *olambdag(k)**5.5/visc(k)
1740
1741 constg2=vf1s*olambdag(k)*olambdag(k)+ &
1742 vf2s*schmidt(k)**0.33334*gam2pt75*sqrt(term0)
1743 tmp2=xnog*constg2
1744 tmp3=term1*oxlf*tmp2-cwoxlf*temcc(k)*( pgacw(k)+pgacr(k) )
1745 tmp4=amin1(0.0,tmp3)
1746 pgmlt(k)=amax1( tmp4,-qgzodt(k) )
1747
1748
1749 !c
1750 !c (4) EVAPORATION OF MELTING GRAUPEL (Pgmltevp) : HR (A19)
1751 !c but use Lin et al. coefficience
1752 !c Pgmltevp is a negative value
1753 !c abg=2.0*pi*(qvoqsiz(k)-1.0)*tmpc/(tmpa+tmpb)
1754 !c
1755 tmpa=rvapor*xka(k)*tem(k)*tem(k)
1756 tmpb=xlv*xlv*rho(k)*qswz(k)*diffwv(k)
1757 tmpc=tmpa*qswz(k)*diffwv(k)
1758 tmpd=amin1( 0.0,(qvoqswz(k)-0.90)*qswz(k)*odtb )
1759
1760 !c
1761 !c abg=2*pi*(Si-1)/rho/(A"+B")
1762 !c
1763 abg=2.0*pi*(qvoqswz(k)-0.90)*tmpc/(tmpa+tmpb)
1764 !
1765 !**** allow evaporation to occur when RH less than 90%
1766 !**** here not using 100% because the evaporation cooling
1767 !**** of temperature is not taking into account yet; hence,
1768 !**** the qgw value is a little bit larger. This will avoid
1769 !**** evaporation can generate cloud.
1770 !
1771 !c vf1s,vf2s=ventilation factors for snow
1772 !c vf1s=0.78,vf2s=0.31 in LIN
1773 !c constg=sqrt(4.*grav*rhograul*0.33334*orho(k)/Cdrag)
1774 !c constg2=vf1s*olambdag(k)*olambdag(k)+
1775 !c vf2s*schmidt(k)**0.33334*gam2pt75*constg
1776 !
1777 tmp2=abg*xnog*constg2
1778 tmp3=amin1(0.0,tmp2)
1779 tmp3=amax1( tmp3,tmpd )
1780 pgmltevp(k)=amax1( tmp3,-qgzodt(k) )
1781
1782 !c
1783 !c (5) ACCRETION OF SNOW BY GRAUPEL (Pgacs) : Lin (29)
1784 !c Compute processes (3) only when qsz > 0.0 and qgz > 0.0
1785 !c egs=1.0
1786 !c
1787 egs=1.
1788 tmpa=olambdas(k)*olambdas(k)
1789 tmpb=olambdag(k)*olambdag(k)
1790 tmpc=olambdas(k)*olambdag(k)
1791 tmp1=pi*pi*xnos*xnog*abs( vts(k)-vtg(k) )*orho(k)
1792 tmp2=tmpa*tmpa*olambdag(k)*(5.0*tmpa+2.0*tmpc+0.5*tmpb)
1793 tmp3=tmp1*egs*rhosnow*tmp2
1794 Pgacs(k)=amin1( tmp3,qszodt(k) )
1795
1796 endif
1797
1798
1799 !
1800 900 continue
1801
1802 !cc
1803 !c
1804 !c**********************************************************************
1805 !c***** combine all processes together and avoid negative *****
1806 !c***** water substances
1807 !***********************************************************************
1808 !c
1809 if ( temcc(k) .lt. 0.0) then
1810 !,delta4,1.-delta4
1811 !c
1812 !c gdelta4=gindex*delta4
1813 !c g1sdelt4=gindex*(1.-delta4)
1814 !c
1815 gdelta4=gindex*delta4
1816 g1sdelt4=gindex*(1.-delta4)
1817 !c
1818 !c combined water vapor depletions
1819 !c
1820 !cc graupel
1821 tmp=psdep(k)+pgdep(k)*gindex
1822 if ( tmp .gt. qvzodt(k) ) then
1823 factor=qvzodt(k)/tmp
1824 psdep(k)=psdep(k)*factor
1825 pgdep(k)=pgdep(k)*factor*gindex
1826 end if
1827 !c
1828 !c combined cloud water depletions
1829 !c
1830 tmp=praut(k)+psacw(k)+psfw(k)+pracw(k)+gindex*Pgacw(k)
1831 if ( tmp .gt. qlzodt(k) ) then
1832 factor=qlzodt(k)/tmp
1833 praut(k)=praut(k)*factor
1834 psacw(k)=psacw(k)*factor
1835 psfw(k)=psfw(k)*factor
1836 pracw(k)=pracw(k)*factor
1837 !cc graupel
1838 Pgacw(k)=Pgacw(k)*factor*gindex
1839 end if
1840 !c
1841 !c combined cloud ice depletions
1842 !c
1843 tmp=psaut(k)+psaci(k)+praci(k)+psfi(k)+Pgaci(k)*gdelta4 &
1844 +Pgacip(k)*g1sdelt4
1845 if (tmp .gt. qizodt(k) ) then
1846 factor=qizodt(k)/tmp
1847 psaut(k)=psaut(k)*factor
1848 psaci(k)=psaci(k)*factor
1849 praci(k)=praci(k)*factor
1850 psfi(k)=psfi(k)*factor
1851 !cc graupel
1852 Pgaci(k)=Pgaci(k)*factor*gdelta4
1853 Pgacip(k)=Pgacip(k)*factor*g1sdelt4
1854 endif
1855 !c
1856 !c combined all rain processes
1857 !c
1858 tmp_r=piacr(k)+psacr(k)-prevp(k)-praut(k)-pracw(k) &
1859 +Pgfr(k)*gindex+Pgacr(k)*gdelta4 &
1860 +Pgacrp(k)*g1sdelt4
1861 if (tmp_r .gt. qrzodt(k) ) then
1862 factor=qrzodt(k)/tmp_r
1863 piacr(k)=piacr(k)*factor
1864 psacr(k)=psacr(k)*factor
1865 prevp(k)=prevp(k)*factor
1866 !cc graupel
1867 Pgfr(k)=Pgfr(k)*factor*gindex
1868 Pgacr(k)=Pgacr(k)*factor*gdelta4
1869 Pgacrp(k)=Pgacrp(k)*factor*g1sdelt4
1870 endif
1871
1872 !c
1873 !c if qrz < 1.0E-4 and qsz < 1.0E-4 then delta2=1.
1874 !c (all Pracs and Psacr become to snow)
1875 !c if qrz >= 1.0E-4 or qsz >= 1.0E-4 then delta2=0.
1876 !c (all Pracs and Psacr become to graupel)
1877 !c
1878 if (qrz(k) .lt. 1.0E-4 .and. qsz(k) .lt. 1.0E-4) then
1879 delta2=1.0
1880 else
1881 delta2=0.0
1882 endif
1883 !
1884 !cc graupel
1885
1886 !c
1887 !c if qrz(k) < 1.0e-4 then delta3=1. means praci(k) --> qs
1888 !c piacr(k) --> qs
1889 !c if qrz(k) > 1.0e-4 then delta3=0. means praci(k) --> qg
1890 !c piacr(k) --> qg : Lin (20)
1891
1892 if (qrz(k) .lt. 1.0e-4) then
1893 delta3=1.0
1894 else
1895 delta3=0.0
1896 endif
1897 !
1898 !c
1899 !c if gindex = 0.(no graupel) then delta2=1.0
1900 !c delta3=1.0
1901 !c
1902 if (gindex .eq. 0.) then
1903 delta2=1.0
1904 delta3=1.0
1905 endif
1906 !
1907 !c
1908 !c combined all snow processes
1909 !c
1910 tmp_s=-pssub(k)-(psaut(k)+psaci(k)+psacw(k)+psfw(k)+ &
1911 psfi(k)+praci(k)*delta3+piacr(k)*delta3+ &
1912 psdep(k))+Pgaut(k)*gindex+Pgacs(k)*gdelta4+ &
1913 Pgacsp(k)*g1sdelt4+Pracs(k)*(1.-delta2)- &
1914 Psacr(k)*delta2
1915 if ( tmp_s .gt. qszodt(k) ) then
1916 factor=qszodt(k)/tmp_s
1917 pssub(k)=pssub(k)*factor
1918 Pracs(k)=Pracs(k)*factor
1919 !cc graupel
1920 Pgaut(k)=Pgaut(k)*factor*gindex
1921 Pgacs(k)=Pgacs(k)*factor*gdelta4
1922 Pgacsp(k)=Pgacsp(k)*factor*g1sdelt4
1923 endif
1924
1925 !cc graupel
1926 !
1927
1928 ! if (gindex .eq. 0.) goto 998
1929 !c
1930 !c combined all graupel processes
1931 !c
1932 tmp_g=-pgaut(k)-pgfr(k)-Pgacw(k)*delta4-Pgaci(k)*delta4 &
1933 -Pgacr(k)*delta4-Pgacs(k)*delta4 &
1934 -pgwet(k)*(1.-delta4)-pgsub(k)-pgdep(k) &
1935 -psacr(k)*(1-delta2)-Pracs(k)*(1-delta2) &
1936 -praci(k)*(1-delta3)-piacr(k)*(1-delta3)
1937 if (tmp_g .gt. qgzodt(k)) then
1938 factor=qgzodt(k)/tmp_g
1939 pgsub(k)=pgsub(k)*factor
1940 endif
1941
1942 998 continue
1943 !c
1944 !c calculate new water substances, thetae, tem, and qvsbar
1945 !c
1946
1947 !cc graupel
1948 pvapor(k)=-pssub(k)-psdep(k)-prevp(k)-pgsub(k)*gindex &
1949 -pgdep(k)*gindex
1950 qvz(k)=amax1( qvmin,qvz(k)+dtb*pvapor(k) )
1951 pclw(k)=-praut(k)-pracw(k)-psacw(k)-psfw(k)-pgacw(k)*gindex
1952 if(flag_qndrop)then
1953 if( qlz(k) > 1e-20 ) &
1954 qndropz(k)=amax1( 0.0,qndropz(k)+dtb*pclw(k)*qndropz(k)/qlz(k) ) !sg
1955 endif
1956 qlz(k)=amax1( 0.0,qlz(k)+dtb*pclw(k) )
1957 pcli(k)=-psaut(k)-psfi(k)-psaci(k)-praci(k)-pgaci(k)*gdelta4 &
1958 -Pgacip(k)*g1sdelt4
1959 qiz(k)=amax1( 0.0,qiz(k)+dtb*pcli(k) )
1960 tmp_r=piacr(k)+psacr(k)-prevp(k)-praut(k)-pracw(k) &
1961 +Pgfr(k)*gindex+Pgacr(k)*gdelta4 &
1962 +Pgacrp(k)*g1sdelt4
1963 232 format(i2,1x,6(f9.3,1x))
1964 prain(k)=-tmp_r
1965 qrz(k)=amax1( 0.0,qrz(k)+dtb*prain(k) )
1966 tmp_s=-pssub(k)-(psaut(k)+psaci(k)+psacw(k)+psfw(k)+ &
1967 psfi(k)+praci(k)*delta3+piacr(k)*delta3+ &
1968 psdep(k))+Pgaut(k)*gindex+Pgacs(k)*gdelta4+ &
1969 Pgacsp(k)*g1sdelt4+Pracs(k)*(1.-delta2)- &
1970 Psacr(k)*delta2
1971 psnow(k)=-tmp_s
1972 qsz(k)=amax1( 0.0,qsz(k)+dtb*psnow(k) )
1973 qschg(k)=qschg(k)+psnow(k)
1974 qschg(k)=psnow(k)
1975 !cc graupel
1976 tmp_g=-pgaut(k)-pgfr(k)-Pgacw(k)*delta4-Pgaci(k)*delta4 &
1977 -Pgacr(k)*delta4-Pgacs(k)*delta4 &
1978 -pgwet(k)*(1.-delta4)-pgsub(k)-pgdep(k) &
1979 -psacr(k)*(1-delta2)-Pracs(k)*(1-delta2) &
1980 -praci(k)*(1-delta3)-piacr(k)*(1-delta3)
1981 252 format(i2,1x,6(f12.9,1x))
1982 262 format(i2,1x,7(f12.9,1x))
1983 pgraupel(k)=-tmp_g
1984 pgraupel(k)=pgraupel(k)*gindex
1985 qgz(k)=amax1( 0.0,qgz(k)+dtb*pgraupel(k))
1986 ! qgchg(k)=qgchg(k)+pgraupel(k)
1987 qgchg(k)=pgraupel(k)
1988 qgz(k)=qgz(k)*gindex
1989
1990 tmp=ocp/tothz(k)*xLf*(qschg(k)+qgchg(k))
1991 theiz(k)=theiz(k)+dtb*tmp
1992 thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
1993 tem(k)=thz(k)*tothz(k)
1994
1995 temcc(k)=tem(k)-273.15
1996
1997 if( temcc(k) .lt. -40.0 ) qswz(k)=qsiz(k)
1998 qlpqi=qlz(k)+qiz(k)
1999 if ( qlpqi .eq. 0.0 ) then
2000 qvsbar(k)=qsiz(k)
2001 else
2002 qvsbar(k)=( qiz(k)*qsiz(k)+qlz(k)*qswz(k) )/qlpqi
2003 endif
2004
2005 !
2006 else
2007 !c
2008 !c combined cloud water depletions
2009 !c
2010 tmp=praut(k)+psacw(k)+pracw(k)+pgacw(k)*gindex
2011 if ( tmp .gt. qlzodt(k) ) then
2012 factor=qlzodt(k)/tmp
2013 praut(k)=praut(k)*factor
2014 psacw(k)=psacw(k)*factor
2015 pracw(k)=pracw(k)*factor
2016 !cc graupel
2017 pgacw(k)=pgacw(k)*factor*gindex
2018 end if
2019 !c
2020 !c combined all snow processes
2021 !c
2022 tmp_s=-(psmlt(k)+psmltevp(k))+Pgacs(k)*gindex
2023 if (tmp_s .gt. qszodt(k) ) then
2024 factor=qszodt(k)/tmp_s
2025 psmlt(k)=psmlt(k)*factor
2026 psmltevp(k)=psmltevp(k)*factor
2027 !cc graupel
2028 Pgacs(k)=Pgacs(k)*factor*gindex
2029 endif
2030
2031 !c
2032 !c
2033 !cc graupel
2034 !c
2035 ! if (gindex .eq. 0.) goto 997
2036
2037 !c
2038 !c combined all graupel processes
2039 !c
2040 tmp_g=-pgmlt(k)-pgacs(k)-pgmltevp(k)
2041 if (tmp_g .gt. qgzodt(k)) then
2042 factor=qgzodt(k)/tmp_g
2043 pgmltevp(k)=pgmltevp(k)*factor
2044 pgmlt(k)=pgmlt(k)*factor
2045 endif
2046 !c
2047 997 continue
2048
2049 !c
2050 !c combined all rain processes
2051 !c
2052 tmp_r=-prevp(k)-(praut(k)+pracw(k)+psacw(k)-psmlt(k)) &
2053 +pgmlt(k)*gindex-pgacw(k)*gindex
2054 if (tmp_r .gt. qrzodt(k) ) then
2055 factor=qrzodt(k)/tmp_r
2056 prevp(k)=prevp(k)*factor
2057 endif
2058 !c
2059 !c
2060 !c calculate new water substances and thetae
2061 !c
2062
2063
2064 pvapor(k)=-psmltevp(k)-prevp(k)-pgmltevp(k)*gindex
2065 qvz(k)=amax1( qvmin,qvz(k)+dtb*pvapor(k))
2066 pclw(k)=-praut(k)-pracw(k)-psacw(k)-pgacw(k)*gindex
2067 if(flag_qndrop)then
2068 if( qlz(k) > 1e-20 ) &
2069 qndropz(k)=amax1( 0.0,qndropz(k)+dtb*pclw(k)*qndropz(k)/qlz(k) ) !sg
2070 endif
2071 qlz(k)=amax1( 0.0,qlz(k)+dtb*pclw(k) )
2072 pcli(k)=0.0
2073 qiz(k)=amax1( 0.0,qiz(k)+dtb*pcli(k) )
2074 tmp_r=-prevp(k)-(praut(k)+pracw(k)+psacw(k)-psmlt(k)) &
2075 +pgmlt(k)*gindex-pgacw(k)*gindex
2076 242 format(i2,1x,7(f9.6,1x))
2077 prain(k)=-tmp_r
2078 tmpqrz=qrz(k)
2079 qrz(k)=amax1( 0.0,qrz(k)+dtb*prain(k) )
2080 tmp_s=-(psmlt(k)+psmltevp(k))+Pgacs(k)*gindex
2081 psnow(k)=-tmp_s
2082 qsz(k)=amax1( 0.0,qsz(k)+dtb*psnow(k) )
2083 ! qschg(k)=qschg(k)+psnow(k)
2084 qschg(k)=psnow(k)
2085 !cc graupel
2086
2087 tmp_g=-pgmlt(k)-pgacs(k)-pgmltevp(k)
2088 ! write(*,272)k,pgmlt(k),pgacs(k),pgmltevp(k),
2089 272 format(i2,1x,3(f12.9,1x))
2090 pgraupel(k)=-tmp_g*gindex
2091 qgz(k)=amax1( 0.0,qgz(k)+dtb*pgraupel(k))
2092 ! qgchg(k)=qgchg(k)+pgraupel(k)
2093 qgchg(k)=pgraupel(k)
2094 qgz(k)=qgz(k)*gindex
2095 !
2096 tmp=ocp/tothz(k)*xLf*(qschg(k)+qgchg(k))
2097 theiz(k)=theiz(k)+dtb*tmp
2098 thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
2099
2100 tem(k)=thz(k)*tothz(k)
2101 temcc(k)=tem(k)-273.15
2102 ! qswz(k)=episp0k*oprez(k)* &
2103 ! exp( svp2*temcc(k)/(tem(k)-svp3) )
2104 es=1000.*svp1*exp( svp2*temcc(k)/(tem(k)-svp3) )
2105 qswz(k)=ep2*es/(prez(k)-es)
2106 qsiz(k)=qswz(k)
2107 qvsbar(k)=qswz(k)
2108 !
2109 end if
2110 preclw(k)=pclw(k) !sg
2111
2112 !
2113 !***********************************************************************
2114 !********** saturation adjustment **********
2115 !***********************************************************************
2116 !
2117 ! allow supersaturation exits linearly from 0% at 500 mb to 50%
2118 ! above 300 mb
2119 ! 5.0e-5=1.0/(500mb-300mb)
2120 !
2121 rsat=1.0+0.5*(50000.0-prez(k))*5.0e-5
2122 rsat=amax1(1.0,rsat)
2123 rsat=amin1(1.5,rsat)
2124 rsat=1.0
2125 if( qvz(k)+qlz(k)+qiz(k) .lt. rsat*qvsbar(k) ) then
2126
2127 !c
2128 !c unsaturated
2129 !c
2130 qvz(k)=qvz(k)+qlz(k)+qiz(k)
2131 qlz(k)=0.0
2132 qiz(k)=0.0
2133
2134 thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
2135 tem(k)=thz(k)*tothz(k)
2136 temcc(k)=tem(k)-273.15
2137
2138 go to 1800
2139 !
2140 else
2141 !c
2142 !c saturated
2143 !c
2144 !
2145 pladj(k)=qlz(k)
2146 piadj(k)=qiz(k)
2147 !
2148
2149 CALL satadj(qvz, qlz, qiz, prez, theiz, thz, tothz, kts, kte, &
2150 k, xLvocp, xLfocp, episp0k, EP2,SVP1,SVP2,SVP3,SVPT0 )
2151
2152 !
2153 pladj(k)=odtb*(qlz(k)-pladj(k))
2154 piadj(k)=odtb*(qiz(k)-piadj(k))
2155 !
2156 pclw(k)=pclw(k)+pladj(k)
2157 pcli(k)=pcli(k)+piadj(k)
2158 pvapor(k)=pvapor(k)-( pladj(k)+piadj(k) )
2159 !
2160 thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
2161 tem(k)=thz(k)*tothz(k)
2162
2163 temcc(k)=tem(k)-273.15
2164
2165 ! qswz(k)=episp0k*oprez(k)* &
2166 ! exp( svp2*temcc(k)/(tem(k)-svp3) )
2167 es=1000.*svp1*exp( svp2*temcc(k)/(tem(k)-svp3) )
2168 qswz(k)=ep2*es/(prez(k)-es)
2169 if (tem(k) .lt. 273.15 ) then
2170 ! qsiz(k)=episp0k*oprez(k)* &
2171 ! exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
2172 es=1000.*svp1*exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
2173 qsiz(k)=ep2*es/(prez(k)-es)
2174 if (temcc(k) .lt. -40.0) qswz(k)=qsiz(k)
2175 else
2176 qsiz(k)=qswz(k)
2177 endif
2178 qlpqi=qlz(k)+qiz(k)
2179 if ( qlpqi .eq. 0.0 ) then
2180 qvsbar(k)=qsiz(k)
2181 else
2182 qvsbar(k)=( qiz(k)*qsiz(k)+qlz(k)*qswz(k) )/qlpqi
2183 endif
2184
2185 end if
2186
2187 !
2188 !***********************************************************************
2189 !***** melting and freezing of cloud ice and cloud water *****
2190 !***********************************************************************
2191 qlpqi=qlz(k)+qiz(k)
2192 if(qlpqi .le. 0.0) go to 1800
2193 !
2194 !c
2195 !c (1) HOMOGENEOUS NUCLEATION WHEN T< -40 C (Pihom)
2196 !c
2197 if(temcc(k) .lt. -40.0) pihom(k)=qlz(k)*odtb
2198 !c
2199 !c (2) MELTING OF ICE CRYSTAL WHEN T> 0 C (Pimlt)
2200 !c
2201 if(temcc(k) .gt. 0.0) pimlt(k)=qiz(k)*odtb
2202 !c
2203 !c (3) PRODUCTION OF CLOUD ICE BY BERGERON PROCESS (Pidw): Hsie (p957)
2204 !c this process only considered when -31 C < T < 0 C
2205 !c
2206 if(temcc(k) .lt. 0.0 .and. temcc(k) .gt. -31.0) then
2207 !c!
2208 !c! parama1 and parama2 functions must be user supplied
2209 !c!
2210 a1=parama1( temcc(k) )
2211 a2=parama2( temcc(k) )
2212 !! change unit from cgs to mks
2213 a1=a1*0.001**(1.0-a2)
2214 xnin=xni0*exp(-bni*temcc(k))
2215 pidw(k)=xnin*orho(k)*(a1*xmnin**a2)
2216 end if
2217 !
2218 pcli(k)=pcli(k)+pihom(k)-pimlt(k)+pidw(k)
2219 pclw(k)=pclw(k)-pihom(k)+pimlt(k)-pidw(k)
2220 qlz(k)=amax1( 0.0,qlz(k)+dtb*(-pihom(k)+pimlt(k)-pidw(k)) )
2221 qiz(k)=amax1( 0.0,qiz(k)+dtb*(pihom(k)-pimlt(k)+pidw(k)) )
2222
2223 !
2224 CALL satadj(qvz, qlz, qiz, prez, theiz, thz, tothz, kts, kte, &
2225 k, xLvocp, xLfocp, episp0k ,EP2,SVP1,SVP2,SVP3,SVPT0)
2226
2227 thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
2228 tem(k)=thz(k)*tothz(k)
2229
2230 temcc(k)=tem(k)-273.15
2231
2232 ! qswz(k)=episp0k*oprez(k)* &
2233 ! exp( svp2*temcc(k)/(tem(k)-svp3) )
2234 es=1000.*svp1*exp( svp2*temcc(k)/(tem(k)-svp3) )
2235 qswz(k)=ep2*es/(prez(k)-es)
2236
2237 if (tem(k) .lt. 273.15 ) then
2238 ! qsiz(k)=episp0k*oprez(k)* &
2239 ! exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
2240 es=1000.*svp1*exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
2241 qsiz(k)=ep2*es/(prez(k)-es)
2242 if (temcc(k) .lt. -40.0) qswz(k)=qsiz(k)
2243 else
2244 qsiz(k)=qswz(k)
2245 endif
2246 qlpqi=qlz(k)+qiz(k)
2247 if ( qlpqi .eq. 0.0 ) then
2248 qvsbar(k)=qsiz(k)
2249 else
2250 qvsbar(k)=( qiz(k)*qsiz(k)+qlz(k)*qswz(k) )/qlpqi
2251 endif
2252
2253 1800 continue
2254 !
2255 !***********************************************************************
2256 !********** integrate the productions of rain and snow **********
2257 !***********************************************************************
2258 !c
2259
2260 2000 continue
2261
2262
2263 !---------------------------------------------------------------------
2264
2265 !
2266 !***********************************************************************
2267 !****** Write terms in cloud physics to time series dataset *****
2268 !***********************************************************************
2269 !
2270 ! open(unit=24,form='formatted',status='new',
2271 ! & file='cloud.dat')
2272
2273 !9030 format(10e12.6)
2274
2275 ! write(24,*)'tmp'
2276 ! write(24,9030) (tem(k),k=kts+1,kte)
2277 ! write(24,*)'qiz'
2278 ! write(24,9030) (qiz(k),k=kts+1,kte)
2279 ! write(24,*)'qsz'
2280 ! write(24,9030) (qsz(k),k=kts+1,kte)
2281 ! write(24,*)'qrz'
2282 ! write(24,9030) (qrz(k),k=kts+1,kte)
2283 ! write(24,*)'qgz'
2284 ! write(24,9030) (qgz(k),k=kts+1,kte)
2285 ! write(24,*)'qvoqsw'
2286 ! write(24,9030) (qvoqswz(k),k=kts+1,kte)
2287 ! write(24,*)'qvoqsi'
2288 ! write(24,9030) (qvoqsiz(k),k=kts+1,kte)
2289 ! write(24,*)'vtr'
2290 ! write(24,9030) (vtr(k),k=kts+1,kte)
2291 ! write(24,*)'vts'
2292 ! write(24,9030) (vts(k),k=kts+1,kte)
2293 ! write(24,*)'vtg'
2294 ! write(24,9030) (vtg(k),k=kts+1,kte)
2295 ! write(24,*)'pclw'
2296 ! write(24,9030) (pclw(k),k=kts+1,kte)
2297 ! write(24,*)'pvapor'
2298 ! write(24,9030) (pvapor(k),k=kts+1,kte)
2299 ! write(24,*)'pcli'
2300 ! write(24,9030) (pcli(k),k=kts+1,kte)
2301 ! write(24,*)'pimlt'
2302 ! write(24,9030) (pimlt(k),k=kts+1,kte)
2303 ! write(24,*)'pihom'
2304 ! write(24,9030) (pihom(k),k=kts+1,kte)
2305 ! write(24,*)'pidw'
2306 ! write(24,9030) (pidw(k),k=kts+1,kte)
2307 ! write(24,*)'prain'
2308 ! write(24,9030) (prain(k),k=kts+1,kte)
2309 ! write(24,*)'praut'
2310 ! write(24,9030) (praut(k),k=kts+1,kte)
2311 ! write(24,*)'pracw'
2312 ! write(24,9030) (pracw(k),k=kts+1,kte)
2313 ! write(24,*)'prevp'
2314 ! write(24,9030) (prevp(k),k=kts+1,kte)
2315 ! write(24,*)'psnow'
2316 ! write(24,9030) (psnow(k),k=kts+1,kte)
2317 ! write(24,*)'psaut'
2318 ! write(24,9030) (psaut(k),k=kts+1,kte)
2319 ! write(24,*)'psfw'
2320 ! write(24,9030) (psfw(k),k=kts+1,kte)
2321 ! write(24,*)'psfi'
2322 ! write(24,9030) (psfi(k),k=kts+1,kte)
2323 ! write(24,*)'praci'
2324 ! write(24,9030) (praci(k),k=kts+1,kte)
2325 ! write(24,*)'piacr'
2326 ! write(24,9030) (piacr(k),k=kts+1,kte)
2327 ! write(24,*)'psaci'
2328 ! write(24,9030) (psaci(k),k=kts+1,kte)
2329 ! write(24,*)'psacw'
2330 ! write(24,9030) (psacw(k),k=kts+1,kte)
2331 ! write(24,*)'psdep'
2332 ! write(24,9030) (psdep(k),k=kts+1,kte)
2333 ! write(24,*)'pssub'
2334 ! write(24,9030) (pssub(k),k=kts+1,kte)
2335 ! write(24,*)'pracs'
2336 ! write(24,9030) (pracs(k),k=kts+1,kte)
2337 ! write(24,*)'psacr'
2338 ! write(24,9030) (psacr(k),k=kts+1,kte)
2339 ! write(24,*)'psmlt'
2340 ! write(24,9030) (psmlt(k),k=kts+1,kte)
2341 ! write(24,*)'psmltevp'
2342 ! write(24,9030) (psmltevp(k),k=kts+1,kte)
2343 ! write(24,*)'pladj'
2344 ! write(24,9030) (pladj(k),k=kts+1,kte)
2345 ! write(24,*)'piadj'
2346 ! write(24,9030) (piadj(k),k=kts+1,kte)
2347 ! write(24,*)'pgraupel'
2348 ! write(24,9030) (pgraupel(k),k=kts+1,kte)
2349 ! write(24,*)'pgaut'
2350 ! write(24,9030) (pgaut(k),k=kts+1,kte)
2351 ! write(24,*)'pgfr'
2352 ! write(24,9030) (pgfr(k),k=kts+1,kte)
2353 ! write(24,*)'pgacw'
2354 ! write(24,9030) (pgacw(k),k=kts+1,kte)
2355 ! write(24,*)'pgaci'
2356 ! write(24,9030) (pgaci(k),k=kts+1,kte)
2357 ! write(24,*)'pgacr'
2358 ! write(24,9030) (pgacr(k),k=kts+1,kte)
2359 ! write(24,*)'pgacs'
2360 ! write(24,9030) (pgacs(k),k=kts+1,kte)
2361 ! write(24,*)'pgacip'
2362 ! write(24,9030) (pgacip(k),k=kts+1,kte)
2363 ! write(24,*)'pgacrP'
2364 ! write(24,9030) (pgacrP(k),k=kts+1,kte)
2365 ! write(24,*)'pgacsp'
2366 ! write(24,9030) (pgacsp(k),k=kts+1,kte)
2367 ! write(24,*)'pgwet'
2368 ! write(24,9030) (pgwet(k),k=kts+1,kte)
2369 ! write(24,*)'pdry'
2370 ! write(24,9030) (pdry(k),k=kts+1,kte)
2371 ! write(24,*)'pgsub'
2372 ! write(24,9030) (pgsub(k),k=kts+1,kte)
2373 ! write(24,*)'pgdep'
2374 ! write(24,9030) (pgdep(k),k=kts+1,kte)
2375 ! write(24,*)'pgmlt'
2376 ! write(24,9030) (pgmlt(k),k=kts+1,kte)
2377 ! write(24,*)'pgmltevp'
2378 ! write(24,9030) (pgmltevp(k),k=kts+1,kte)
2379
2380
2381
2382 !**** below if qv < qvmin then qv=qvmin, ql=0.0, and qi=0.0
2383 !
2384 do k=kts+1,kte
2385 if ( qvz(k) .lt. qvmin ) then
2386 qlz(k)=0.0
2387 qiz(k)=0.0
2388 qvz(k)=amax1( qvmin,qvz(k)+qlz(k)+qiz(k) )
2389 end if
2390 enddo
2391 !
2392 END SUBROUTINE clphy1d
2393
2394
2395 !---------------------------------------------------------------------
2396 ! SATURATED ADJUSTMENT
2397 !---------------------------------------------------------------------
2398 SUBROUTINE satadj(qvz, qlz, qiz, prez, theiz, thz, tothz, &
2399 kts, kte, k, xLvocp, xLfocp, episp0k, EP2,SVP1,SVP2,SVP3,SVPT0)
2400 !---------------------------------------------------------------------
2401 IMPLICIT NONE
2402 !---------------------------------------------------------------------
2403 ! This program use Newton's method for finding saturated temperature
2404 ! and saturation mixing ratio.
2405 !
2406 ! In this saturation adjustment scheme we assume
2407 ! (1) the saturation mixing ratio is the mass weighted average of
2408 ! saturation values over liquid water (qsw), and ice (qsi)
2409 ! following Lord et al., 1984 and Tao, 1989
2410 !
2411 ! (2) the percentage of cloud liquid and cloud ice will
2412 ! be fixed during the saturation calculation
2413 !---------------------------------------------------------------------
2414 !
2415
2416 INTEGER, INTENT(IN ) :: kts, kte, k
2417
2418 REAL, DIMENSION( kts:kte ), &
2419 INTENT(INOUT) :: qvz, qlz, qiz
2420 !
2421 REAL, DIMENSION( kts:kte ), &
2422 INTENT(IN ) :: prez, theiz, tothz
2423
2424 REAL, INTENT(IN ) :: xLvocp, xLfocp, episp0k
2425 REAL, INTENT(IN ) :: EP2,SVP1,SVP2,SVP3,SVPT0
2426
2427 ! LOCAL VARS
2428
2429 INTEGER :: n
2430
2431 REAL, DIMENSION( kts:kte ) :: thz, tem, temcc, qsiz, &
2432 qswz, qvsbar
2433
2434 REAL :: qsat, qlpqi, ratql, t0, t1, tmp1, ratqi, tsat, absft, &
2435 denom1, denom2, dqvsbar, ftsat, dftsat, qpz, &
2436 gindex, es
2437 !
2438 !---------------------------------------------------------------------
2439
2440 thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
2441
2442 tem(k)=tothz(k)*thz(k)
2443 if (tem(k) .gt. 273.15) then
2444 ! qsat=episp0k/prez(k)* &
2445 ! exp( svp2*(tem(k)-273.15)/(tem(k)-svp3) )
2446 es=1000.*svp1*exp( svp2*(tem(k)-svpt0)/(tem(k)-svp3) )
2447 qsat=ep2*es/(prez(k)-es)
2448 else
2449 qsat=episp0k/prez(k)* &
2450 exp( 21.8745584*(tem(k)-273.15)/(tem(k)-7.66) )
2451 end if
2452 qpz=qvz(k)+qlz(k)+qiz(k)
2453 if (qpz .lt. qsat) then
2454 qvz(k)=qpz
2455 qiz(k)=0.0
2456 qlz(k)=0.0
2457 go to 400
2458 end if
2459 qlpqi=qlz(k)+qiz(k)
2460 if( qlpqi .ge. 1.0e-5) then
2461 ratql=qlz(k)/qlpqi
2462 ratqi=qiz(k)/qlpqi
2463 else
2464 t0=273.15
2465 ! t1=233.15
2466 t1=248.15
2467 tmp1=( t0-tem(k) )/(t0-t1)
2468 tmp1=amin1(1.0,tmp1)
2469 tmp1=amax1(0.0,tmp1)
2470 ratqi=tmp1
2471 ratql=1.0-tmp1
2472 end if
2473 !
2474 !
2475 !-- saturation mixing ratios over water and ice
2476 !-- at the outset we will follow Bolton 1980 MWR for
2477 !-- the water and Murray JAS 1967 for the ice
2478 !
2479 !-- dqvsbar=d(qvsbar)/dT
2480 !-- ftsat=F(Tsat)
2481 !-- dftsat=d(F(T))/dT
2482 !
2483 ! First guess of tsat
2484
2485 tsat=tem(k)
2486 absft=1.0
2487 !
2488 do 200 n=1,20
2489 denom1=1.0/(tsat-svp3)
2490 denom2=1.0/(tsat-7.66)
2491 ! qswz(k)=episp0k/prez(k)* &
2492 ! exp( svp2*denom1*(tsat-273.15) )
2493 es=1000.*svp1*exp( svp2*denom1*(tsat-svpt0) )
2494 qswz(k)=ep2*es/(prez(k)-es)
2495 if (tem(k) .lt. 273.15) then
2496 ! qsiz(k)=episp0k/prez(k)* &
2497 ! exp( 21.8745584*denom2*(tsat-273.15) )
2498 es=1000.*svp1*exp( 21.8745584*denom2*(tsat-273.15) )
2499 qsiz(k)=ep2*es/(prez(k)-es)
2500 if (tem(k) .lt. 233.15) qswz(k)=qsiz(k)
2501 else
2502 qsiz(k)=qswz(k)
2503 endif
2504 qvsbar(k)=ratql*qswz(k)+ratqi*qsiz(k)
2505 !
2506 ! if( absft .lt. 0.01 .and. n .gt. 3 ) go to 300
2507 if( absft .lt. 0.01 ) go to 300
2508 !
2509 dqvsbar=ratql*qswz(k)*svp2*243.5*denom1*denom1+ &
2510 ratqi*qsiz(k)*21.8745584*265.5*denom2*denom2
2511 ftsat=tsat+(xlvocp+ratqi*xlfocp)*qvsbar(k)- &
2512 tothz(k)*theiz(k)-xlfocp*ratqi*(qvz(k)+qlz(k)+qiz(k))
2513 dftsat=1.0+(xlvocp+ratqi*xlfocp)*dqvsbar
2514 tsat=tsat-ftsat/dftsat
2515 absft=abs(ftsat)
2516
2517 200 continue
2518 9020 format(1x,'point can not converge, absft,n=',e12.5,i5)
2519 !
2520 300 continue
2521 if( qpz .gt. qvsbar(k) ) then
2522 qvz(k)=qvsbar(k)
2523 qiz(k)=ratqi*( qpz-qvz(k) )
2524 qlz(k)=ratql*( qpz-qvz(k) )
2525 else
2526 qvz(k)=qpz
2527 qiz(k)=0.0
2528 qlz(k)=0.0
2529 end if
2530 400 continue
2531
2532 END SUBROUTINE satadj
2533
2534
2535 !----------------------------------------------------------------
2536 REAL FUNCTION parama1(temp)
2537 !----------------------------------------------------------------
2538 IMPLICIT NONE
2539 !----------------------------------------------------------------
2540 ! This program calculate the parameter for crystal growth rate
2541 ! in Bergeron process
2542 !----------------------------------------------------------------
2543
2544 REAL, INTENT (IN ) :: temp
2545 REAL, DIMENSION(32) :: a1
2546 INTEGER :: i1, i1p1
2547 REAL :: ratio
2548
2549 data a1/0.100e-10,0.7939e-7,0.7841e-6,0.3369e-5,0.4336e-5, &
2550 0.5285e-5,0.3728e-5,0.1852e-5,0.2991e-6,0.4248e-6, &
2551 0.7434e-6,0.1812e-5,0.4394e-5,0.9145e-5,0.1725e-4, &
2552 0.3348e-4,0.1725e-4,0.9175e-5,0.4412e-5,0.2252e-5, &
2553 0.9115e-6,0.4876e-6,0.3473e-6,0.4758e-6,0.6306e-6, &
2554 0.8573e-6,0.7868e-6,0.7192e-6,0.6513e-6,0.5956e-6, &
2555 0.5333e-6,0.4834e-6/
2556
2557 i1=int(-temp)+1
2558 i1p1=i1+1
2559 ratio=-(temp)-float(i1-1)
2560 parama1=a1(i1)+ratio*( a1(i1p1)-a1(i1) )
2561
2562 END FUNCTION parama1
2563
2564 !----------------------------------------------------------------
2565 REAL FUNCTION parama2(temp)
2566 !----------------------------------------------------------------
2567 IMPLICIT NONE
2568 !----------------------------------------------------------------
2569 ! This program calculate the parameter for crystal growth rate
2570 ! in Bergeron process
2571 !----------------------------------------------------------------
2572
2573 REAL, INTENT (IN ) :: temp
2574 REAL, DIMENSION(32) :: a2
2575 INTEGER :: i1, i1p1
2576 REAL :: ratio
2577
2578 data a2/0.0100,0.4006,0.4831,0.5320,0.5307,0.5319,0.5249, &
2579 0.4888,0.3849,0.4047,0.4318,0.4771,0.5183,0.5463, &
2580 0.5651,0.5813,0.5655,0.5478,0.5203,0.4906,0.4447, &
2581 0.4126,0.3960,0.4149,0.4320,0.4506,0.4483,0.4460, &
2582 0.4433,0.4413,0.4382,0.4361/
2583 i1=int(-temp)+1
2584 i1p1=i1+1
2585 ratio=-(temp)-float(i1-1)
2586 parama2=a2(i1)+ratio*( a2(i1p1)-a2(i1) )
2587
2588 END FUNCTION parama2
2589
2590 !----------------------------------------------------------------
2591 REAL FUNCTION ggamma(X)
2592 !----------------------------------------------------------------
2593 IMPLICIT NONE
2594 !----------------------------------------------------------------
2595 REAL, INTENT(IN ) :: x
2596 REAL, DIMENSION(8) :: B
2597 INTEGER ::j, K1
2598 REAL ::PF, G1TO2 ,TEMP
2599
2600 DATA B/-.577191652,.988205891,-.897056937,.918206857, &
2601 -.756704078,.482199394,-.193527818,.035868343/
2602
2603 PF=1.
2604 TEMP=X
2605 DO 10 J=1,200
2606 IF (TEMP .LE. 2) GO TO 20
2607 TEMP=TEMP-1.
2608 10 PF=PF*TEMP
2609 100 FORMAT(//,5X,'module_mp_lin: INPUT TO GAMMA FUNCTION TOO LARGE, X=',E12.5)
2610 WRITE(wrf_err_message,100)X
2611 CALL wrf_error_fatal(wrf_err_message)
2612 20 G1TO2=1.
2613 TEMP=TEMP - 1.
2614 DO 30 K1=1,8
2615 30 G1TO2=G1TO2 + B(K1)*TEMP**K1
2616 ggamma=PF*G1TO2
2617
2618 END FUNCTION ggamma
2619
2620 !----------------------------------------------------------------
2621
2622 END MODULE module_mp_lin
2623