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 tmp2=abg*xnog*constg2
1687 pgdep(k)=amax1(0.0,tmp2)
1688 pgsub(k)=amin1(0.0,tmp2)
1689 pgsub(k)=amax1( pgsub(k),-qgzodt(k) )
1690
1691 4000 continue
1692 else
1693 !
1694 !***********************************************************************
1695 !********* graupel production processes for T > 0 C **********
1696 !***********************************************************************
1697 !
1698 !c
1699 !c (1) ACCRETION OF CLOUD WATER BY GRAUPEL (Pgacw): Lin (40)
1700 !c egw=1.0
1701 !c Cdrag=0.6 drag coefficients for hairstone
1702 !c constg=sqrt(4.*grav*rhograul*0.33334*orho(k)/Cdrag)
1703
1704 egw=1.0
1705 constg=sqrt(4.*grav*rhograul*0.33334*orho(k)*oCdrag)
1706 tmp1=pio4*xnog*gam3pt5*constg*olambdag(k)**3.5
1707 tmp2=qlz(k)*egw*tmp1
1708 Pgacw(k)=amin1( tmp2,qlzodt(k) )
1709
1710 !c
1711 !c (2) ACCRETION OF RAIN BY GRAUPEL (Pgacr): Lin (42)
1712 !c Compute processes (5) only when qrz > 0.0 and qgz > 0.0
1713 !c egr=1.
1714 !c
1715 egr=1.
1716 tmpa=olambdar(k)*olambdar(k)
1717 tmpb=olambdag(k)*olambdag(k)
1718 tmpc=olambdar(k)*olambdag(k)
1719 tmp1=pi*pi*xnor*xnog*abs( vtr(k)-vtg(k) )*orho(k)
1720 tmp2=tmpa*tmpa*olambdag(k)*(5.0*tmpa+2.0*tmpc+0.5*tmpb)
1721 tmp3=tmp1*egr*rhowater*tmp2
1722 pgacr(k)=amin1( tmp3,qrzodt(k) )
1723
1724
1725 !c
1726 !c (3) GRAUPEL MELTING TO FORM RAIN (Pgmlt): Lin (47)
1727 !c Pgmlt is negative value
1728 !c constg=sqrt(4.*grav*rhograul*0.33334*orho(k)/Cdrag)
1729 !c constg2=vf1s*olambdag(k)*olambdag(k)+
1730 !c vf2s*schmidt(k)**0.33334*gam2pt75*constg
1731 !c Cdrag=0.6 drag coefficients for hairstone
1732 !
1733 delrs=rs0(k)-qvz(k)
1734 term1=2.0*pi*orho(k)*( xlv*diffwv(k)*rho(k)*delrs- &
1735 xka(k)*temcc(k) )
1736 term0=sqrt(4.*grav*rhograul*0.33334*orho(k)*ocdrag) &
1737 *olambdag(k)**5.5/visc(k)
1738
1739 constg2=vf1s*olambdag(k)*olambdag(k)+ &
1740 vf2s*schmidt(k)**0.33334*gam2pt75*sqrt(term0)
1741 tmp2=xnog*constg2
1742 tmp3=term1*oxlf*tmp2-cwoxlf*temcc(k)*( pgacw(k)+pgacr(k) )
1743 tmp4=amin1(0.0,tmp3)
1744 pgmlt(k)=amax1( tmp4,-qgzodt(k) )
1745
1746
1747 !c
1748 !c (4) EVAPORATION OF MELTING GRAUPEL (Pgmltevp) : HR (A19)
1749 !c but use Lin et al. coefficience
1750 !c Pgmltevp is a negative value
1751 !c abg=2.0*pi*(qvoqsiz(k)-1.0)*tmpc/(tmpa+tmpb)
1752 !c
1753 tmpa=rvapor*xka(k)*tem(k)*tem(k)
1754 tmpb=xlv*xlv*rho(k)*qswz(k)*diffwv(k)
1755 tmpc=tmpa*qswz(k)*diffwv(k)
1756 tmpd=amin1( 0.0,(qvoqswz(k)-0.90)*qswz(k)*odtb )
1757
1758 !c
1759 !c abg=2*pi*(Si-1)/rho/(A"+B")
1760 !c
1761 abg=2.0*pi*(qvoqswz(k)-0.90)*tmpc/(tmpa+tmpb)
1762 !
1763 !**** allow evaporation to occur when RH less than 90%
1764 !**** here not using 100% because the evaporation cooling
1765 !**** of temperature is not taking into account yet; hence,
1766 !**** the qgw value is a little bit larger. This will avoid
1767 !**** evaporation can generate cloud.
1768 !
1769 !c vf1s,vf2s=ventilation factors for snow
1770 !c vf1s=0.78,vf2s=0.31 in LIN
1771 !c constg=sqrt(4.*grav*rhograul*0.33334*orho(k)/Cdrag)
1772 !c constg2=vf1s*olambdag(k)*olambdag(k)+
1773 !c vf2s*schmidt(k)**0.33334*gam2pt75*constg
1774 !
1775 tmp2=abg*xnog*constg2
1776 tmp3=amin1(0.0,tmp2)
1777 tmp3=amax1( tmp3,tmpd )
1778 pgmltevp(k)=amax1( tmp3,-qgzodt(k) )
1779
1780 !c
1781 !c (5) ACCRETION OF SNOW BY GRAUPEL (Pgacs) : Lin (29)
1782 !c Compute processes (3) only when qsz > 0.0 and qgz > 0.0
1783 !c egs=1.0
1784 !c
1785 egs=1.
1786 tmpa=olambdas(k)*olambdas(k)
1787 tmpb=olambdag(k)*olambdag(k)
1788 tmpc=olambdas(k)*olambdag(k)
1789 tmp1=pi*pi*xnos*xnog*abs( vts(k)-vtg(k) )*orho(k)
1790 tmp2=tmpa*tmpa*olambdag(k)*(5.0*tmpa+2.0*tmpc+0.5*tmpb)
1791 tmp3=tmp1*egs*rhosnow*tmp2
1792 Pgacs(k)=amin1( tmp3,qszodt(k) )
1793
1794 endif
1795
1796
1797 !
1798 900 continue
1799
1800 !cc
1801 !c
1802 !c**********************************************************************
1803 !c***** combine all processes together and avoid negative *****
1804 !c***** water substances
1805 !***********************************************************************
1806 !c
1807 if ( temcc(k) .lt. 0.0) then
1808 !,delta4,1.-delta4
1809 !c
1810 !c gdelta4=gindex*delta4
1811 !c g1sdelt4=gindex*(1.-delta4)
1812 !c
1813 gdelta4=gindex*delta4
1814 g1sdelt4=gindex*(1.-delta4)
1815 !c
1816 !c combined water vapor depletions
1817 !c
1818 !cc graupel
1819 tmp=psdep(k)+pgdep(k)*gindex
1820 if ( tmp .gt. qvzodt(k) ) then
1821 factor=qvzodt(k)/tmp
1822 psdep(k)=psdep(k)*factor
1823 pgdep(k)=pgdep(k)*factor*gindex
1824 end if
1825 !c
1826 !c combined cloud water depletions
1827 !c
1828 tmp=praut(k)+psacw(k)+psfw(k)+pracw(k)+gindex*Pgacw(k)
1829 if ( tmp .gt. qlzodt(k) ) then
1830 factor=qlzodt(k)/tmp
1831 praut(k)=praut(k)*factor
1832 psacw(k)=psacw(k)*factor
1833 psfw(k)=psfw(k)*factor
1834 pracw(k)=pracw(k)*factor
1835 !cc graupel
1836 Pgacw(k)=Pgacw(k)*factor*gindex
1837 end if
1838 !c
1839 !c combined cloud ice depletions
1840 !c
1841 tmp=psaut(k)+psaci(k)+praci(k)+psfi(k)+Pgaci(k)*gdelta4 &
1842 +Pgacip(k)*g1sdelt4
1843 if (tmp .gt. qizodt(k) ) then
1844 factor=qizodt(k)/tmp
1845 psaut(k)=psaut(k)*factor
1846 psaci(k)=psaci(k)*factor
1847 praci(k)=praci(k)*factor
1848 psfi(k)=psfi(k)*factor
1849 !cc graupel
1850 Pgaci(k)=Pgaci(k)*factor*gdelta4
1851 Pgacip(k)=Pgacip(k)*factor*g1sdelt4
1852 endif
1853 !c
1854 !c combined all rain processes
1855 !c
1856 tmp_r=piacr(k)+psacr(k)-prevp(k)-praut(k)-pracw(k) &
1857 +Pgfr(k)*gindex+Pgacr(k)*gdelta4 &
1858 +Pgacrp(k)*g1sdelt4
1859 if (tmp_r .gt. qrzodt(k) ) then
1860 factor=qrzodt(k)/tmp_r
1861 piacr(k)=piacr(k)*factor
1862 psacr(k)=psacr(k)*factor
1863 prevp(k)=prevp(k)*factor
1864 !cc graupel
1865 Pgfr(k)=Pgfr(k)*factor*gindex
1866 Pgacr(k)=Pgacr(k)*factor*gdelta4
1867 Pgacrp(k)=Pgacrp(k)*factor*g1sdelt4
1868 endif
1869
1870 !c
1871 !c if qrz < 1.0E-4 and qsz < 1.0E-4 then delta2=1.
1872 !c (all Pracs and Psacr become to snow)
1873 !c if qrz >= 1.0E-4 or qsz >= 1.0E-4 then delta2=0.
1874 !c (all Pracs and Psacr become to graupel)
1875 !c
1876 if (qrz(k) .lt. 1.0E-4 .and. qsz(k) .lt. 1.0E-4) then
1877 delta2=1.0
1878 else
1879 delta2=0.0
1880 endif
1881 !
1882 !cc graupel
1883
1884 !c
1885 !c if qrz(k) < 1.0e-4 then delta3=1. means praci(k) --> qs
1886 !c piacr(k) --> qs
1887 !c if qrz(k) > 1.0e-4 then delta3=0. means praci(k) --> qg
1888 !c piacr(k) --> qg : Lin (20)
1889
1890 if (qrz(k) .lt. 1.0e-4) then
1891 delta3=1.0
1892 else
1893 delta3=0.0
1894 endif
1895 !
1896 !c
1897 !c if gindex = 0.(no graupel) then delta2=1.0
1898 !c delta3=1.0
1899 !c
1900 if (gindex .eq. 0.) then
1901 delta2=1.0
1902 delta3=1.0
1903 endif
1904 !
1905 !c
1906 !c combined all snow processes
1907 !c
1908 tmp_s=-pssub(k)-(psaut(k)+psaci(k)+psacw(k)+psfw(k)+ &
1909 psfi(k)+praci(k)*delta3+piacr(k)*delta3+ &
1910 psdep(k))+Pgaut(k)*gindex+Pgacs(k)*gdelta4+ &
1911 Pgacsp(k)*g1sdelt4+Pracs(k)*(1.-delta2)- &
1912 Psacr(k)*delta2
1913 if ( tmp_s .gt. qszodt(k) ) then
1914 factor=qszodt(k)/tmp_s
1915 pssub(k)=pssub(k)*factor
1916 Pracs(k)=Pracs(k)*factor
1917 !cc graupel
1918 Pgaut(k)=Pgaut(k)*factor*gindex
1919 Pgacs(k)=Pgacs(k)*factor*gdelta4
1920 Pgacsp(k)=Pgacsp(k)*factor*g1sdelt4
1921 endif
1922
1923 !cc graupel
1924 !
1925
1926 ! if (gindex .eq. 0.) goto 998
1927 !c
1928 !c combined all graupel processes
1929 !c
1930 tmp_g=-pgaut(k)-pgfr(k)-Pgacw(k)*delta4-Pgaci(k)*delta4 &
1931 -Pgacr(k)*delta4-Pgacs(k)*delta4 &
1932 -pgwet(k)*(1.-delta4)-pgsub(k)-pgdep(k) &
1933 -psacr(k)*(1-delta2)-Pracs(k)*(1-delta2) &
1934 -praci(k)*(1-delta3)-piacr(k)*(1-delta3)
1935 if (tmp_g .gt. qgzodt(k)) then
1936 factor=qgzodt(k)/tmp_g
1937 pgsub(k)=pgsub(k)*factor
1938 endif
1939
1940 998 continue
1941 !c
1942 !c calculate new water substances, thetae, tem, and qvsbar
1943 !c
1944
1945 !cc graupel
1946 pvapor(k)=-pssub(k)-psdep(k)-prevp(k)-pgsub(k)*gindex &
1947 -pgdep(k)*gindex
1948 qvz(k)=amax1( qvmin,qvz(k)+dtb*pvapor(k) )
1949 pclw(k)=-praut(k)-pracw(k)-psacw(k)-psfw(k)-pgacw(k)*gindex
1950 if(flag_qndrop)then
1951 if( qlz(k) > 1e-20 ) &
1952 qndropz(k)=amax1( 0.0,qndropz(k)+dtb*pclw(k)*qndropz(k)/qlz(k) ) !sg
1953 endif
1954 qlz(k)=amax1( 0.0,qlz(k)+dtb*pclw(k) )
1955 pcli(k)=-psaut(k)-psfi(k)-psaci(k)-praci(k)-pgaci(k)*gdelta4 &
1956 -Pgacip(k)*g1sdelt4
1957 qiz(k)=amax1( 0.0,qiz(k)+dtb*pcli(k) )
1958 tmp_r=piacr(k)+psacr(k)-prevp(k)-praut(k)-pracw(k) &
1959 +Pgfr(k)*gindex+Pgacr(k)*gdelta4 &
1960 +Pgacrp(k)*g1sdelt4
1961 232 format(i2,1x,6(f9.3,1x))
1962 prain(k)=-tmp_r
1963 qrz(k)=amax1( 0.0,qrz(k)+dtb*prain(k) )
1964 tmp_s=-pssub(k)-(psaut(k)+psaci(k)+psacw(k)+psfw(k)+ &
1965 psfi(k)+praci(k)*delta3+piacr(k)*delta3+ &
1966 psdep(k))+Pgaut(k)*gindex+Pgacs(k)*gdelta4+ &
1967 Pgacsp(k)*g1sdelt4+Pracs(k)*(1.-delta2)- &
1968 Psacr(k)*delta2
1969 psnow(k)=-tmp_s
1970 qsz(k)=amax1( 0.0,qsz(k)+dtb*psnow(k) )
1971 qschg(k)=qschg(k)+psnow(k)
1972 qschg(k)=psnow(k)
1973 !cc graupel
1974 tmp_g=-pgaut(k)-pgfr(k)-Pgacw(k)*delta4-Pgaci(k)*delta4 &
1975 -Pgacr(k)*delta4-Pgacs(k)*delta4 &
1976 -pgwet(k)*(1.-delta4)-pgsub(k)-pgdep(k) &
1977 -psacr(k)*(1-delta2)-Pracs(k)*(1-delta2) &
1978 -praci(k)*(1-delta3)-piacr(k)*(1-delta3)
1979 252 format(i2,1x,6(f12.9,1x))
1980 262 format(i2,1x,7(f12.9,1x))
1981 pgraupel(k)=-tmp_g
1982 pgraupel(k)=pgraupel(k)*gindex
1983 qgz(k)=amax1( 0.0,qgz(k)+dtb*pgraupel(k))
1984 ! qgchg(k)=qgchg(k)+pgraupel(k)
1985 qgchg(k)=pgraupel(k)
1986 qgz(k)=qgz(k)*gindex
1987
1988 tmp=ocp/tothz(k)*xLf*(qschg(k)+qgchg(k))
1989 theiz(k)=theiz(k)+dtb*tmp
1990 thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
1991 tem(k)=thz(k)*tothz(k)
1992
1993 temcc(k)=tem(k)-273.15
1994
1995 if( temcc(k) .lt. -40.0 ) qswz(k)=qsiz(k)
1996 qlpqi=qlz(k)+qiz(k)
1997 if ( qlpqi .eq. 0.0 ) then
1998 qvsbar(k)=qsiz(k)
1999 else
2000 qvsbar(k)=( qiz(k)*qsiz(k)+qlz(k)*qswz(k) )/qlpqi
2001 endif
2002
2003 !
2004 else
2005 !c
2006 !c combined cloud water depletions
2007 !c
2008 tmp=praut(k)+psacw(k)+pracw(k)+pgacw(k)*gindex
2009 if ( tmp .gt. qlzodt(k) ) then
2010 factor=qlzodt(k)/tmp
2011 praut(k)=praut(k)*factor
2012 psacw(k)=psacw(k)*factor
2013 pracw(k)=pracw(k)*factor
2014 !cc graupel
2015 pgacw(k)=pgacw(k)*factor*gindex
2016 end if
2017 !c
2018 !c combined all snow processes
2019 !c
2020 tmp_s=-(psmlt(k)+psmltevp(k))+Pgacs(k)*gindex
2021 if (tmp_s .gt. qszodt(k) ) then
2022 factor=qszodt(k)/tmp_s
2023 psmlt(k)=psmlt(k)*factor
2024 psmltevp(k)=psmltevp(k)*factor
2025 !cc graupel
2026 Pgacs(k)=Pgacs(k)*factor*gindex
2027 endif
2028
2029 !c
2030 !c
2031 !cc graupel
2032 !c
2033 ! if (gindex .eq. 0.) goto 997
2034
2035 !c
2036 !c combined all graupel processes
2037 !c
2038 tmp_g=-pgmlt(k)-pgacs(k)-pgmltevp(k)
2039 if (tmp_g .gt. qgzodt(k)) then
2040 factor=qgzodt(k)/tmp_g
2041 pgmltevp(k)=pgmltevp(k)*factor
2042 pgmlt(k)=pgmlt(k)*factor
2043 endif
2044 !c
2045 997 continue
2046
2047 !c
2048 !c combined all rain processes
2049 !c
2050 tmp_r=-prevp(k)-(praut(k)+pracw(k)+psacw(k)-psmlt(k)) &
2051 +pgmlt(k)*gindex-pgacw(k)*gindex
2052 if (tmp_r .gt. qrzodt(k) ) then
2053 factor=qrzodt(k)/tmp_r
2054 prevp(k)=prevp(k)*factor
2055 endif
2056 !c
2057 !c
2058 !c calculate new water substances and thetae
2059 !c
2060
2061
2062 pvapor(k)=-psmltevp(k)-prevp(k)-pgmltevp(k)*gindex
2063 qvz(k)=amax1( qvmin,qvz(k)+dtb*pvapor(k))
2064 pclw(k)=-praut(k)-pracw(k)-psacw(k)-pgacw(k)*gindex
2065 if(flag_qndrop)then
2066 if( qlz(k) > 1e-20 ) &
2067 qndropz(k)=amax1( 0.0,qndropz(k)+dtb*pclw(k)*qndropz(k)/qlz(k) ) !sg
2068 endif
2069 qlz(k)=amax1( 0.0,qlz(k)+dtb*pclw(k) )
2070 pcli(k)=0.0
2071 qiz(k)=amax1( 0.0,qiz(k)+dtb*pcli(k) )
2072 tmp_r=-prevp(k)-(praut(k)+pracw(k)+psacw(k)-psmlt(k)) &
2073 +pgmlt(k)*gindex-pgacw(k)*gindex
2074 242 format(i2,1x,7(f9.6,1x))
2075 prain(k)=-tmp_r
2076 tmpqrz=qrz(k)
2077 qrz(k)=amax1( 0.0,qrz(k)+dtb*prain(k) )
2078 tmp_s=-(psmlt(k)+psmltevp(k))+Pgacs(k)*gindex
2079 psnow(k)=-tmp_s
2080 qsz(k)=amax1( 0.0,qsz(k)+dtb*psnow(k) )
2081 ! qschg(k)=qschg(k)+psnow(k)
2082 qschg(k)=psnow(k)
2083 !cc graupel
2084
2085 tmp_g=-pgmlt(k)-pgacs(k)-pgmltevp(k)
2086 ! write(*,272)k,pgmlt(k),pgacs(k),pgmltevp(k),
2087 272 format(i2,1x,3(f12.9,1x))
2088 pgraupel(k)=-tmp_g*gindex
2089 qgz(k)=amax1( 0.0,qgz(k)+dtb*pgraupel(k))
2090 ! qgchg(k)=qgchg(k)+pgraupel(k)
2091 qgchg(k)=pgraupel(k)
2092 qgz(k)=qgz(k)*gindex
2093 !
2094 tmp=ocp/tothz(k)*xLf*(qschg(k)+qgchg(k))
2095 theiz(k)=theiz(k)+dtb*tmp
2096 thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
2097
2098 tem(k)=thz(k)*tothz(k)
2099 temcc(k)=tem(k)-273.15
2100 ! qswz(k)=episp0k*oprez(k)* &
2101 ! exp( svp2*temcc(k)/(tem(k)-svp3) )
2102 es=1000.*svp1*exp( svp2*temcc(k)/(tem(k)-svp3) )
2103 qswz(k)=ep2*es/(prez(k)-es)
2104 qsiz(k)=qswz(k)
2105 qvsbar(k)=qswz(k)
2106 !
2107 end if
2108 preclw(k)=pclw(k) !sg
2109
2110 !
2111 !***********************************************************************
2112 !********** saturation adjustment **********
2113 !***********************************************************************
2114 !
2115 ! allow supersaturation exits linearly from 0% at 500 mb to 50%
2116 ! above 300 mb
2117 ! 5.0e-5=1.0/(500mb-300mb)
2118 !
2119 rsat=1.0+0.5*(50000.0-prez(k))*5.0e-5
2120 rsat=amax1(1.0,rsat)
2121 rsat=amin1(1.5,rsat)
2122 rsat=1.0
2123 if( qvz(k)+qlz(k)+qiz(k) .lt. rsat*qvsbar(k) ) then
2124
2125 !c
2126 !c unsaturated
2127 !c
2128 qvz(k)=qvz(k)+qlz(k)+qiz(k)
2129 qlz(k)=0.0
2130 qiz(k)=0.0
2131
2132 thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
2133 tem(k)=thz(k)*tothz(k)
2134 temcc(k)=tem(k)-273.15
2135
2136 go to 1800
2137 !
2138 else
2139 !c
2140 !c saturated
2141 !c
2142 !
2143 pladj(k)=qlz(k)
2144 piadj(k)=qiz(k)
2145 !
2146
2147 CALL satadj(qvz, qlz, qiz, prez, theiz, thz, tothz, kts, kte, &
2148 k, xLvocp, xLfocp, episp0k, EP2,SVP1,SVP2,SVP3,SVPT0 )
2149
2150 !
2151 pladj(k)=odtb*(qlz(k)-pladj(k))
2152 piadj(k)=odtb*(qiz(k)-piadj(k))
2153 !
2154 pclw(k)=pclw(k)+pladj(k)
2155 pcli(k)=pcli(k)+piadj(k)
2156 pvapor(k)=pvapor(k)-( pladj(k)+piadj(k) )
2157 !
2158 thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
2159 tem(k)=thz(k)*tothz(k)
2160
2161 temcc(k)=tem(k)-273.15
2162
2163 ! qswz(k)=episp0k*oprez(k)* &
2164 ! exp( svp2*temcc(k)/(tem(k)-svp3) )
2165 es=1000.*svp1*exp( svp2*temcc(k)/(tem(k)-svp3) )
2166 qswz(k)=ep2*es/(prez(k)-es)
2167 if (tem(k) .lt. 273.15 ) then
2168 ! qsiz(k)=episp0k*oprez(k)* &
2169 ! exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
2170 es=1000.*svp1*exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
2171 qsiz(k)=ep2*es/(prez(k)-es)
2172 if (temcc(k) .lt. -40.0) qswz(k)=qsiz(k)
2173 else
2174 qsiz(k)=qswz(k)
2175 endif
2176 qlpqi=qlz(k)+qiz(k)
2177 if ( qlpqi .eq. 0.0 ) then
2178 qvsbar(k)=qsiz(k)
2179 else
2180 qvsbar(k)=( qiz(k)*qsiz(k)+qlz(k)*qswz(k) )/qlpqi
2181 endif
2182
2183 end if
2184
2185 !
2186 !***********************************************************************
2187 !***** melting and freezing of cloud ice and cloud water *****
2188 !***********************************************************************
2189 qlpqi=qlz(k)+qiz(k)
2190 if(qlpqi .le. 0.0) go to 1800
2191 !
2192 !c
2193 !c (1) HOMOGENEOUS NUCLEATION WHEN T< -40 C (Pihom)
2194 !c
2195 if(temcc(k) .lt. -40.0) pihom(k)=qlz(k)*odtb
2196 !c
2197 !c (2) MELTING OF ICE CRYSTAL WHEN T> 0 C (Pimlt)
2198 !c
2199 if(temcc(k) .gt. 0.0) pimlt(k)=qiz(k)*odtb
2200 !c
2201 !c (3) PRODUCTION OF CLOUD ICE BY BERGERON PROCESS (Pidw): Hsie (p957)
2202 !c this process only considered when -31 C < T < 0 C
2203 !c
2204 if(temcc(k) .lt. 0.0 .and. temcc(k) .gt. -31.0) then
2205 !c!
2206 !c! parama1 and parama2 functions must be user supplied
2207 !c!
2208 a1=parama1( temcc(k) )
2209 a2=parama2( temcc(k) )
2210 !! change unit from cgs to mks
2211 a1=a1*0.001**(1.0-a2)
2212 xnin=xni0*exp(-bni*temcc(k))
2213 pidw(k)=xnin*orho(k)*(a1*xmnin**a2)
2214 end if
2215 !
2216 pcli(k)=pcli(k)+pihom(k)-pimlt(k)+pidw(k)
2217 pclw(k)=pclw(k)-pihom(k)+pimlt(k)-pidw(k)
2218 qlz(k)=amax1( 0.0,qlz(k)+dtb*(-pihom(k)+pimlt(k)-pidw(k)) )
2219 qiz(k)=amax1( 0.0,qiz(k)+dtb*(pihom(k)-pimlt(k)+pidw(k)) )
2220
2221 !
2222 CALL satadj(qvz, qlz, qiz, prez, theiz, thz, tothz, kts, kte, &
2223 k, xLvocp, xLfocp, episp0k ,EP2,SVP1,SVP2,SVP3,SVPT0)
2224
2225 thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
2226 tem(k)=thz(k)*tothz(k)
2227
2228 temcc(k)=tem(k)-273.15
2229
2230 ! qswz(k)=episp0k*oprez(k)* &
2231 ! exp( svp2*temcc(k)/(tem(k)-svp3) )
2232 es=1000.*svp1*exp( svp2*temcc(k)/(tem(k)-svp3) )
2233 qswz(k)=ep2*es/(prez(k)-es)
2234
2235 if (tem(k) .lt. 273.15 ) then
2236 ! qsiz(k)=episp0k*oprez(k)* &
2237 ! exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
2238 es=1000.*svp1*exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
2239 qsiz(k)=ep2*es/(prez(k)-es)
2240 if (temcc(k) .lt. -40.0) qswz(k)=qsiz(k)
2241 else
2242 qsiz(k)=qswz(k)
2243 endif
2244 qlpqi=qlz(k)+qiz(k)
2245 if ( qlpqi .eq. 0.0 ) then
2246 qvsbar(k)=qsiz(k)
2247 else
2248 qvsbar(k)=( qiz(k)*qsiz(k)+qlz(k)*qswz(k) )/qlpqi
2249 endif
2250
2251 1800 continue
2252 !
2253 !***********************************************************************
2254 !********** integrate the productions of rain and snow **********
2255 !***********************************************************************
2256 !c
2257
2258 2000 continue
2259
2260
2261 !---------------------------------------------------------------------
2262
2263 !
2264 !***********************************************************************
2265 !****** Write terms in cloud physics to time series dataset *****
2266 !***********************************************************************
2267 !
2268 ! open(unit=24,form='formatted',status='new',
2269 ! & file='cloud.dat')
2270
2271 !9030 format(10e12.6)
2272
2273 ! write(24,*)'tmp'
2274 ! write(24,9030) (tem(k),k=kts+1,kte)
2275 ! write(24,*)'qiz'
2276 ! write(24,9030) (qiz(k),k=kts+1,kte)
2277 ! write(24,*)'qsz'
2278 ! write(24,9030) (qsz(k),k=kts+1,kte)
2279 ! write(24,*)'qrz'
2280 ! write(24,9030) (qrz(k),k=kts+1,kte)
2281 ! write(24,*)'qgz'
2282 ! write(24,9030) (qgz(k),k=kts+1,kte)
2283 ! write(24,*)'qvoqsw'
2284 ! write(24,9030) (qvoqswz(k),k=kts+1,kte)
2285 ! write(24,*)'qvoqsi'
2286 ! write(24,9030) (qvoqsiz(k),k=kts+1,kte)
2287 ! write(24,*)'vtr'
2288 ! write(24,9030) (vtr(k),k=kts+1,kte)
2289 ! write(24,*)'vts'
2290 ! write(24,9030) (vts(k),k=kts+1,kte)
2291 ! write(24,*)'vtg'
2292 ! write(24,9030) (vtg(k),k=kts+1,kte)
2293 ! write(24,*)'pclw'
2294 ! write(24,9030) (pclw(k),k=kts+1,kte)
2295 ! write(24,*)'pvapor'
2296 ! write(24,9030) (pvapor(k),k=kts+1,kte)
2297 ! write(24,*)'pcli'
2298 ! write(24,9030) (pcli(k),k=kts+1,kte)
2299 ! write(24,*)'pimlt'
2300 ! write(24,9030) (pimlt(k),k=kts+1,kte)
2301 ! write(24,*)'pihom'
2302 ! write(24,9030) (pihom(k),k=kts+1,kte)
2303 ! write(24,*)'pidw'
2304 ! write(24,9030) (pidw(k),k=kts+1,kte)
2305 ! write(24,*)'prain'
2306 ! write(24,9030) (prain(k),k=kts+1,kte)
2307 ! write(24,*)'praut'
2308 ! write(24,9030) (praut(k),k=kts+1,kte)
2309 ! write(24,*)'pracw'
2310 ! write(24,9030) (pracw(k),k=kts+1,kte)
2311 ! write(24,*)'prevp'
2312 ! write(24,9030) (prevp(k),k=kts+1,kte)
2313 ! write(24,*)'psnow'
2314 ! write(24,9030) (psnow(k),k=kts+1,kte)
2315 ! write(24,*)'psaut'
2316 ! write(24,9030) (psaut(k),k=kts+1,kte)
2317 ! write(24,*)'psfw'
2318 ! write(24,9030) (psfw(k),k=kts+1,kte)
2319 ! write(24,*)'psfi'
2320 ! write(24,9030) (psfi(k),k=kts+1,kte)
2321 ! write(24,*)'praci'
2322 ! write(24,9030) (praci(k),k=kts+1,kte)
2323 ! write(24,*)'piacr'
2324 ! write(24,9030) (piacr(k),k=kts+1,kte)
2325 ! write(24,*)'psaci'
2326 ! write(24,9030) (psaci(k),k=kts+1,kte)
2327 ! write(24,*)'psacw'
2328 ! write(24,9030) (psacw(k),k=kts+1,kte)
2329 ! write(24,*)'psdep'
2330 ! write(24,9030) (psdep(k),k=kts+1,kte)
2331 ! write(24,*)'pssub'
2332 ! write(24,9030) (pssub(k),k=kts+1,kte)
2333 ! write(24,*)'pracs'
2334 ! write(24,9030) (pracs(k),k=kts+1,kte)
2335 ! write(24,*)'psacr'
2336 ! write(24,9030) (psacr(k),k=kts+1,kte)
2337 ! write(24,*)'psmlt'
2338 ! write(24,9030) (psmlt(k),k=kts+1,kte)
2339 ! write(24,*)'psmltevp'
2340 ! write(24,9030) (psmltevp(k),k=kts+1,kte)
2341 ! write(24,*)'pladj'
2342 ! write(24,9030) (pladj(k),k=kts+1,kte)
2343 ! write(24,*)'piadj'
2344 ! write(24,9030) (piadj(k),k=kts+1,kte)
2345 ! write(24,*)'pgraupel'
2346 ! write(24,9030) (pgraupel(k),k=kts+1,kte)
2347 ! write(24,*)'pgaut'
2348 ! write(24,9030) (pgaut(k),k=kts+1,kte)
2349 ! write(24,*)'pgfr'
2350 ! write(24,9030) (pgfr(k),k=kts+1,kte)
2351 ! write(24,*)'pgacw'
2352 ! write(24,9030) (pgacw(k),k=kts+1,kte)
2353 ! write(24,*)'pgaci'
2354 ! write(24,9030) (pgaci(k),k=kts+1,kte)
2355 ! write(24,*)'pgacr'
2356 ! write(24,9030) (pgacr(k),k=kts+1,kte)
2357 ! write(24,*)'pgacs'
2358 ! write(24,9030) (pgacs(k),k=kts+1,kte)
2359 ! write(24,*)'pgacip'
2360 ! write(24,9030) (pgacip(k),k=kts+1,kte)
2361 ! write(24,*)'pgacrP'
2362 ! write(24,9030) (pgacrP(k),k=kts+1,kte)
2363 ! write(24,*)'pgacsp'
2364 ! write(24,9030) (pgacsp(k),k=kts+1,kte)
2365 ! write(24,*)'pgwet'
2366 ! write(24,9030) (pgwet(k),k=kts+1,kte)
2367 ! write(24,*)'pdry'
2368 ! write(24,9030) (pdry(k),k=kts+1,kte)
2369 ! write(24,*)'pgsub'
2370 ! write(24,9030) (pgsub(k),k=kts+1,kte)
2371 ! write(24,*)'pgdep'
2372 ! write(24,9030) (pgdep(k),k=kts+1,kte)
2373 ! write(24,*)'pgmlt'
2374 ! write(24,9030) (pgmlt(k),k=kts+1,kte)
2375 ! write(24,*)'pgmltevp'
2376 ! write(24,9030) (pgmltevp(k),k=kts+1,kte)
2377
2378
2379
2380 !**** below if qv < qvmin then qv=qvmin, ql=0.0, and qi=0.0
2381 !
2382 do k=kts+1,kte
2383 if ( qvz(k) .lt. qvmin ) then
2384 qlz(k)=0.0
2385 qiz(k)=0.0
2386 qvz(k)=amax1( qvmin,qvz(k)+qlz(k)+qiz(k) )
2387 end if
2388 enddo
2389 !
2390 END SUBROUTINE clphy1d
2391
2392
2393 !---------------------------------------------------------------------
2394 ! SATURATED ADJUSTMENT
2395 !---------------------------------------------------------------------
2396 SUBROUTINE satadj(qvz, qlz, qiz, prez, theiz, thz, tothz, &
2397 kts, kte, k, xLvocp, xLfocp, episp0k, EP2,SVP1,SVP2,SVP3,SVPT0)
2398 !---------------------------------------------------------------------
2399 IMPLICIT NONE
2400 !---------------------------------------------------------------------
2401 ! This program use Newton's method for finding saturated temperature
2402 ! and saturation mixing ratio.
2403 !
2404 ! In this saturation adjustment scheme we assume
2405 ! (1) the saturation mixing ratio is the mass weighted average of
2406 ! saturation values over liquid water (qsw), and ice (qsi)
2407 ! following Lord et al., 1984 and Tao, 1989
2408 !
2409 ! (2) the percentage of cloud liquid and cloud ice will
2410 ! be fixed during the saturation calculation
2411 !---------------------------------------------------------------------
2412 !
2413
2414 INTEGER, INTENT(IN ) :: kts, kte, k
2415
2416 REAL, DIMENSION( kts:kte ), &
2417 INTENT(INOUT) :: qvz, qlz, qiz
2418 !
2419 REAL, DIMENSION( kts:kte ), &
2420 INTENT(IN ) :: prez, theiz, tothz
2421
2422 REAL, INTENT(IN ) :: xLvocp, xLfocp, episp0k
2423 REAL, INTENT(IN ) :: EP2,SVP1,SVP2,SVP3,SVPT0
2424
2425 ! LOCAL VARS
2426
2427 INTEGER :: n
2428
2429 REAL, DIMENSION( kts:kte ) :: thz, tem, temcc, qsiz, &
2430 qswz, qvsbar
2431
2432 REAL :: qsat, qlpqi, ratql, t0, t1, tmp1, ratqi, tsat, absft, &
2433 denom1, denom2, dqvsbar, ftsat, dftsat, qpz, &
2434 gindex, es
2435 !
2436 !---------------------------------------------------------------------
2437
2438 thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
2439
2440 tem(k)=tothz(k)*thz(k)
2441 if (tem(k) .gt. 273.15) then
2442 ! qsat=episp0k/prez(k)* &
2443 ! exp( svp2*(tem(k)-273.15)/(tem(k)-svp3) )
2444 es=1000.*svp1*exp( svp2*(tem(k)-svpt0)/(tem(k)-svp3) )
2445 qsat=ep2*es/(prez(k)-es)
2446 else
2447 qsat=episp0k/prez(k)* &
2448 exp( 21.8745584*(tem(k)-273.15)/(tem(k)-7.66) )
2449 end if
2450 qpz=qvz(k)+qlz(k)+qiz(k)
2451 if (qpz .lt. qsat) then
2452 qvz(k)=qpz
2453 qiz(k)=0.0
2454 qlz(k)=0.0
2455 go to 400
2456 end if
2457 qlpqi=qlz(k)+qiz(k)
2458 if( qlpqi .ge. 1.0e-5) then
2459 ratql=qlz(k)/qlpqi
2460 ratqi=qiz(k)/qlpqi
2461 else
2462 t0=273.15
2463 ! t1=233.15
2464 t1=248.15
2465 tmp1=( t0-tem(k) )/(t0-t1)
2466 tmp1=amin1(1.0,tmp1)
2467 tmp1=amax1(0.0,tmp1)
2468 ratqi=tmp1
2469 ratql=1.0-tmp1
2470 end if
2471 !
2472 !
2473 !-- saturation mixing ratios over water and ice
2474 !-- at the outset we will follow Bolton 1980 MWR for
2475 !-- the water and Murray JAS 1967 for the ice
2476 !
2477 !-- dqvsbar=d(qvsbar)/dT
2478 !-- ftsat=F(Tsat)
2479 !-- dftsat=d(F(T))/dT
2480 !
2481 ! First guess of tsat
2482
2483 tsat=tem(k)
2484 absft=1.0
2485 !
2486 do 200 n=1,20
2487 denom1=1.0/(tsat-svp3)
2488 denom2=1.0/(tsat-7.66)
2489 ! qswz(k)=episp0k/prez(k)* &
2490 ! exp( svp2*denom1*(tsat-273.15) )
2491 es=1000.*svp1*exp( svp2*denom1*(tsat-svpt0) )
2492 qswz(k)=ep2*es/(prez(k)-es)
2493 if (tem(k) .lt. 273.15) then
2494 ! qsiz(k)=episp0k/prez(k)* &
2495 ! exp( 21.8745584*denom2*(tsat-273.15) )
2496 es=1000.*svp1*exp( 21.8745584*denom2*(tsat-273.15) )
2497 qsiz(k)=ep2*es/(prez(k)-es)
2498 if (tem(k) .lt. 233.15) qswz(k)=qsiz(k)
2499 else
2500 qsiz(k)=qswz(k)
2501 endif
2502 qvsbar(k)=ratql*qswz(k)+ratqi*qsiz(k)
2503 !
2504 ! if( absft .lt. 0.01 .and. n .gt. 3 ) go to 300
2505 if( absft .lt. 0.01 ) go to 300
2506 !
2507 dqvsbar=ratql*qswz(k)*svp2*243.5*denom1*denom1+ &
2508 ratqi*qsiz(k)*21.8745584*265.5*denom2*denom2
2509 ftsat=tsat+(xlvocp+ratqi*xlfocp)*qvsbar(k)- &
2510 tothz(k)*theiz(k)-xlfocp*ratqi*(qvz(k)+qlz(k)+qiz(k))
2511 dftsat=1.0+(xlvocp+ratqi*xlfocp)*dqvsbar
2512 tsat=tsat-ftsat/dftsat
2513 absft=abs(ftsat)
2514
2515 200 continue
2516 9020 format(1x,'point can not converge, absft,n=',e12.5,i5)
2517 !
2518 300 continue
2519 if( qpz .gt. qvsbar(k) ) then
2520 qvz(k)=qvsbar(k)
2521 qiz(k)=ratqi*( qpz-qvz(k) )
2522 qlz(k)=ratql*( qpz-qvz(k) )
2523 else
2524 qvz(k)=qpz
2525 qiz(k)=0.0
2526 qlz(k)=0.0
2527 end if
2528 400 continue
2529
2530 END SUBROUTINE satadj
2531
2532
2533 !----------------------------------------------------------------
2534 REAL FUNCTION parama1(temp)
2535 !----------------------------------------------------------------
2536 IMPLICIT NONE
2537 !----------------------------------------------------------------
2538 ! This program calculate the parameter for crystal growth rate
2539 ! in Bergeron process
2540 !----------------------------------------------------------------
2541
2542 REAL, INTENT (IN ) :: temp
2543 REAL, DIMENSION(32) :: a1
2544 INTEGER :: i1, i1p1
2545 REAL :: ratio
2546
2547 data a1/0.100e-10,0.7939e-7,0.7841e-6,0.3369e-5,0.4336e-5, &
2548 0.5285e-5,0.3728e-5,0.1852e-5,0.2991e-6,0.4248e-6, &
2549 0.7434e-6,0.1812e-5,0.4394e-5,0.9145e-5,0.1725e-4, &
2550 0.3348e-4,0.1725e-4,0.9175e-5,0.4412e-5,0.2252e-5, &
2551 0.9115e-6,0.4876e-6,0.3473e-6,0.4758e-6,0.6306e-6, &
2552 0.8573e-6,0.7868e-6,0.7192e-6,0.6513e-6,0.5956e-6, &
2553 0.5333e-6,0.4834e-6/
2554
2555 i1=int(-temp)+1
2556 i1p1=i1+1
2557 ratio=-(temp)-float(i1-1)
2558 parama1=a1(i1)+ratio*( a1(i1p1)-a1(i1) )
2559
2560 END FUNCTION parama1
2561
2562 !----------------------------------------------------------------
2563 REAL FUNCTION parama2(temp)
2564 !----------------------------------------------------------------
2565 IMPLICIT NONE
2566 !----------------------------------------------------------------
2567 ! This program calculate the parameter for crystal growth rate
2568 ! in Bergeron process
2569 !----------------------------------------------------------------
2570
2571 REAL, INTENT (IN ) :: temp
2572 REAL, DIMENSION(32) :: a2
2573 INTEGER :: i1, i1p1
2574 REAL :: ratio
2575
2576 data a2/0.0100,0.4006,0.4831,0.5320,0.5307,0.5319,0.5249, &
2577 0.4888,0.3849,0.4047,0.4318,0.4771,0.5183,0.5463, &
2578 0.5651,0.5813,0.5655,0.5478,0.5203,0.4906,0.4447, &
2579 0.4126,0.3960,0.4149,0.4320,0.4506,0.4483,0.4460, &
2580 0.4433,0.4413,0.4382,0.4361/
2581 i1=int(-temp)+1
2582 i1p1=i1+1
2583 ratio=-(temp)-float(i1-1)
2584 parama2=a2(i1)+ratio*( a2(i1p1)-a2(i1) )
2585
2586 END FUNCTION parama2
2587
2588 !----------------------------------------------------------------
2589 REAL FUNCTION ggamma(X)
2590 !----------------------------------------------------------------
2591 IMPLICIT NONE
2592 !----------------------------------------------------------------
2593 REAL, INTENT(IN ) :: x
2594 REAL, DIMENSION(8) :: B
2595 INTEGER ::j, K1
2596 REAL ::PF, G1TO2 ,TEMP
2597
2598 DATA B/-.577191652,.988205891,-.897056937,.918206857, &
2599 -.756704078,.482199394,-.193527818,.035868343/
2600
2601 PF=1.
2602 TEMP=X
2603 DO 10 J=1,200
2604 IF (TEMP .LE. 2) GO TO 20
2605 TEMP=TEMP-1.
2606 10 PF=PF*TEMP
2607 100 FORMAT(//,5X,'module_mp_lin: INPUT TO GAMMA FUNCTION TOO LARGE, X=',E12.5)
2608 WRITE(wrf_err_message,100)X
2609 CALL wrf_error_fatal(wrf_err_message)
2610 20 G1TO2=1.
2611 TEMP=TEMP - 1.
2612 DO 30 K1=1,8
2613 30 G1TO2=G1TO2 + B(K1)*TEMP**K1
2614 ggamma=PF*G1TO2
2615
2616 END FUNCTION ggamma
2617
2618 !----------------------------------------------------------------
2619
2620 END MODULE module_mp_lin
2621