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