| blob | 3c0d716f740c428c375ffe1857ff0e1c38ee3d69 |
1 !WRF:MODEL_LAYER:PHYSICS
3 !--- The code is based on Lin and Colle (A New Bulk Microphysical Scheme
4 ! that Includes Riming Intensity and Temperature Dependent Ice Characteristics, 2011, MWR)
5 ! and Lin et al. (Parameterization of riming intensity and its impact on ice fall speed using ARM data, 2011, MWR)
6 !--- NOTE: 1) Prognose variables are: qi,PI(precipitating ice, qs, which includes snow, partially rimed snow and graupel),qw,qr
7 !--- 2) Sedimentation flux is based on Prudue Lin scheme
8 !--- 2) PI has varying properties depending on riming intensity (Ri, diagnosed currently following Lin et al. (2011, MWR) and T
9 !--- 3) Autoconverion is based on Liu and Daum (2004)
10 !--- 4) PI size distribution assuming Gamma distribution, but mu_s=0 (Exponential) currently
11 !--- 5) No density dependent fall speed since the V-D is derived using Best number approach, which already includes density effect
12 !--- 6) Future work will include radar equivalent reflectivity using the new PI property (A-D, M-D, N(D)). If you use RIP for reflectivity
13 !--- computation, please note that snow is (1-Ri)*qs and graupel is Ri*qs. Otherwise, reflectivity will be underestimated.
14 !--- 7) The Liu and Daum autoconverion is quite sensitive on Nt_c. For mixed-phase cloud and marine environment, Nt_c of 10 or 20 is suggested.
15 !--- default value is 10E.6. Change accordingly for your use.
18 MODULE module_mp_sbu_ylin
19 USE module_wrf_error
20 !
21 !..Parameters user might change based on their need
22 REAL, PARAMETER, PRIVATE :: RH = 1.0
23 REAL, PARAMETER, PRIVATE :: xnor = 8.0e6
24 REAL, PARAMETER, PRIVATE :: Nt_c = 10.E6
25 !..Water vapor and air gas constants at constant pressure
26 REAL, PARAMETER, PRIVATE :: Rvapor = 461.5
27 REAL, PARAMETER, PRIVATE :: oRv = 1./Rvapor
28 REAL, PARAMETER, PRIVATE :: Rair = 287.04
29 REAL, PARAMETER, PRIVATE :: Cp = 1004.0
30 REAL, PARAMETER, PRIVATE :: grav = 9.81
31 REAL, PARAMETER, PRIVATE :: rhowater = 1000.0
32 REAL, PARAMETER, PRIVATE :: rhosnow = 100.0
34 REAL, PARAMETER, PRIVATE :: SVP1=0.6112
35 REAL, PARAMETER, PRIVATE :: SVP2=17.67
36 REAL, PARAMETER, PRIVATE :: SVP3=29.65
37 REAL, PARAMETER, PRIVATE :: SVPT0=273.15
38 REAL, PARAMETER, PRIVATE :: EP1=Rvapor/Rair-1.
39 REAL, PARAMETER, PRIVATE :: EP2=Rair/Rvapor
40 !..Enthalpy of sublimation, vaporization, and fusion at 0C.
41 REAL, PARAMETER, PRIVATE :: XLS = 2.834E6
42 REAL, PARAMETER, PRIVATE :: XLV = 2.5E6
43 REAL, PARAMETER, PRIVATE :: XLF = XLS - XLV
45 !
46 REAL, PARAMETER, PRIVATE :: &
47 qi0 = 1.0e-3, & !--- ice aggregation to snow threshold
48 xmi50 = 4.8e-10, xmi40 = 2.46e-10, &
49 xni0 = 1.0e-2, xmnin = 1.05e-18, bni = 0.5, &
50 di50 = 1.0e-4, xmi = 4.19e-13, & !--- parameters used in BF process
51 bv_r = 0.8, bv_i = 0.25, &
52 o6 = 1./6., cdrag = 0.6, &
53 avisc = 1.49628e-6, adiffwv = 8.7602e-5, &
54 axka = 1.4132e3, cw = 4.187e3, ci = 2.093e3
55 CONTAINS
57 !-------------------------------------------------------------------
58 ! Lin et al., 1983, JAM, 1065-1092, and
59 ! Rutledge and Hobbs, 1984, JAS, 2949-2972
60 !-------------------------------------------------------------------
61 SUBROUTINE sbu_ylin(th &
62 ,qv, ql, qr &
63 ,qi, qs, Ri3D &
64 ,rho, pii, p &
65 ,dt_in &
66 ,z,ht, dz8w &
67 ,RAINNC, RAINNCV &
68 ,ids,ide, jds,jde, kds,kde &
69 ,ims,ime, jms,jme, kms,kme &
70 ,its,ite, jts,jte, kts,kte &
71 )
72 !-------------------------------------------------------------------
73 IMPLICIT NONE
74 !-------------------------------------------------------------------
75 !
76 !
77 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
78 ims,ime, jms,jme, kms,kme , &
79 its,ite, jts,jte, kts,kte
81 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
82 INTENT(INOUT) :: &
83 th, &
84 qv, &
85 qi,ql, &
86 qs,qr
88 ! YLIN
89 ! Adding RI3D as a variable to the interface
90 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
91 INTENT(INOUT) :: Ri3D
94 !
95 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
96 INTENT(IN ) :: &
97 rho, &
98 pii, &
99 z,p, &
100 dz8w
103 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: ht
104 REAL, INTENT(IN ) :: dt_in
106 REAL, DIMENSION( ims:ime , jms:jme ), &
107 INTENT(INOUT) :: RAINNC, &
108 RAINNCV
110 ! LOCAL VAR
112 INTEGER :: min_q, max_q
114 REAL, DIMENSION( its:ite , jts:jte ) &
115 :: rain, snow,ice
117 REAL, DIMENSION( kts:kte ) :: qvz, qlz, qrz, &
118 qiz, qsz, qgz, &
119 thz, &
120 tothz, rhoz, &
121 orhoz, sqrhoz, &
122 prez, zz, &
123 dzw
125 ! Added vertical profile of Ri (riz) as a variable
126 REAL, DIMENSION( kts:kte ) :: riz
128 !
129 REAL :: dt, pptice, pptrain, pptsnow, pptgraul, rhoe_s
131 INTEGER :: i,j,k
133 dt=dt_in
134 rhoe_s=1.29
136 j_loop: DO j = jts, jte
137 i_loop: DO i = its, ite
138 !
139 !- write data from 3-D to 1-D
140 !
141 DO k = kts, kte
142 qvz(k)=qv(i,k,j)
143 qlz(k)=ql(i,k,j)
144 qrz(k)=qr(i,k,j)
145 qiz(k)=qi(i,k,j)
146 qsz(k)=qs(i,k,j)
147 thz(k)=th(i,k,j)
148 rhoz(k)=rho(i,k,j)
149 orhoz(k)=1./rhoz(k)
150 prez(k)=p(i,k,j)
151 ! sqrhoz(k)=sqrt(rhoe_s*orhoz(k))
152 ! no density dependence of fall speed as Note #5, you can turn it on to increase fall speed at low pressure.
153 sqrhoz(k)=1.0
154 tothz(k)=pii(i,k,j)
155 zz(k)=z(i,k,j)
156 dzw(k)=dz8w(i,k,j)
157 END DO
158 !
159 pptrain=0.
160 pptsnow=0.
161 pptice =0.
163 ! CALL wrf_debug ( 100 , 'microphysics_driver: calling clphy1d_ylin' )
165 CALL clphy1d_ylin( dt, qvz, qlz, qrz, qiz, qsz, &
166 thz, tothz, rhoz, orhoz, sqrhoz, &
167 prez, zz, dzw, ht(I,J), &
168 pptrain, pptsnow, pptice, &
169 kts, kte, i, j, riz )
171 !
172 ! Precipitation from cloud microphysics -- only for one time step
173 !
174 ! unit is transferred from m to mm
175 !
176 rain(i,j)= pptrain
177 snow(i,j)= pptsnow
178 ice(i,j) = pptice
179 !
180 RAINNCV(i,j)= pptrain + pptsnow + pptice
181 RAINNC(i,j) = RAINNC(i,j) + pptrain + pptsnow + pptice
183 !
184 !- update data from 1-D back to 3-D
185 !
186 DO k = kts, kte
187 qv(i,k,j)=qvz(k)
188 ql(i,k,j)=qlz(k)
189 qr(i,k,j)=qrz(k)
190 th(i,k,j)=thz(k)
191 qi(i,k,j)=qiz(k)
192 qs(i,k,j)=qsz(k)
193 ri3d(i,k,j)=riz(k)
194 END DO
195 !
196 ENDDO i_loop
197 ENDDO j_loop
199 END SUBROUTINE sbu_ylin
202 !-----------------------------------------------------------------------
203 SUBROUTINE clphy1d_ylin(dt, qvz, qlz, qrz, qiz, qsz, &
204 thz, tothz, rho, orho, sqrho, &
205 prez, zz, dzw, zsfc, pptrain, pptsnow,pptice, &
206 kts, kte, i, j,riz )
207 !-----------------------------------------------------------------------
208 IMPLICIT NONE
209 !-----------------------------------------------------------------------
210 ! This program handles the vertical 1-D cloud micphysics
211 !-----------------------------------------------------------------------
212 ! avisc: constant in empirical formular for dynamic viscosity of air
213 ! =1.49628e-6 [kg/m/s] = 1.49628e-5 [g/cm/s]
214 ! adiffwv: constant in empirical formular for diffusivity of water
215 ! vapor in air
216 ! = 8.7602e-5 [kgm/s3] = 8.7602 [gcm/s3]
217 ! axka: constant in empirical formular for thermal conductivity of air
218 ! = 1.4132e3 [m2/s2/K] = 1.4132e7 [cm2/s2/K]
219 ! qi0: mixing ratio threshold for cloud ice aggregation [kg/kg]
220 ! xmi50: mass of a 50 micron ice crystal
221 ! = 4.8e-10 [kg] =4.8e-7 [g]
222 ! xmi40: mass of a 40 micron ice crystal
223 ! = 2.46e-10 [kg] = 2.46e-7 [g]
224 ! di50: diameter of a 50 micro (radius) ice crystal
225 ! =1.0e-4 [m]
226 ! xmi: mass of one cloud ice crystal
227 ! =4.19e-13 [kg] = 4.19e-10 [g]
228 ! oxmi=1.0/xmi
229 !
230 ! xni0=1.0e-2 [m-3] The value given in Lin et al. is wrong.(see
231 ! Hsie et al.(1980) and Rutledge and Hobbs(1983) )
232 ! bni=0.5 [K-1]
233 ! xmnin: mass of a natural ice nucleus
234 ! = 1.05e-18 [kg] = 1.05e-15 [g] This values is suggested by
235 ! Hsie et al. (1980)
236 ! = 1.0e-12 [kg] suggested by Rutlegde and Hobbs (1983)
238 ! av_r: av_r in empirical formular for terminal
239 ! velocity of raindrop
240 ! =2115.0 [cm**(1-b)/s] = 2115.0*0.01**(1-b) [m**(1-b)/s]
241 ! bv_r: bv_r in empirical formular for terminal
242 ! velocity of raindrop
243 ! =0.8
244 ! av_i: av_i in empirical formular for terminal
245 ! velocity of snow
246 ! =152.93 [cm**(1-d)/s] = 152.93*0.01**(1-d) [m**(1-d)/s]
247 ! bv_i: bv_i in empirical formular for terminal
248 ! velocity of snow
249 ! =0.25
250 ! vf1r: ventilation factors for rain =0.78
251 ! vf2r: ventilation factors for rain =0.31
252 ! vf1s: ventilation factors for snow =0.65
253 ! vf2s: ventilation factors for snow =0.44
254 !
255 !----------------------------------------------------------------------
257 INTEGER, INTENT(IN ) :: kts, kte, i, j
259 REAL, DIMENSION( kts:kte ), &
260 INTENT(INOUT) :: qvz, qlz, qrz, qiz, qsz, &
261 thz
263 REAL, DIMENSION( kts:kte ), &
264 INTENT(IN ) :: tothz, rho, orho, sqrho, &
265 prez, zz, dzw
268 REAL, INTENT(INOUT) :: pptrain, pptsnow, pptice
270 REAL, INTENT(IN ) :: dt, zsfc
272 ! local vars
274 REAL :: obp4, bp3, bp5, bp6, odp4, &
275 dp3, dp5, dp5o2
277 ! temperary vars
279 REAL :: tmp, tmp0, tmp1, tmp2,tmp3, &
280 tmp4, tmpa,tmpb,tmpc,tmpd,alpha1, &
281 qic, abi,abr, abg, odtberg, &
282 vti50,eiw,eri,esi,esr, esw, &
283 erw,delrs,term0,term1, &
284 Ap, Bp, &
285 factor, tmp_r, tmp_s,tmp_g, &
286 qlpqi, rsat, a1, a2, xnin
288 !
289 REAL, DIMENSION( kts:kte ) :: oprez, tem, temcc, theiz, qswz, &
290 qsiz, qvoqswz, qvoqsiz, qvzodt, &
291 qlzodt, qizodt, qszodt, qrzodt
293 !--- microphysical processes
295 REAL, DIMENSION( kts:kte ) :: psnow, psaut, psfw, psfi, praci, &
296 piacr, psaci, psacw, psdep, pssub, &
297 pracs, psacr, psmlt, psmltevp, &
298 prain, praut, pracw, prevp, pvapor, &
299 pclw, pladj, pcli, pimlt, pihom, &
300 pidw, piadj, pgfr, &
301 qschg
302 !
304 REAL, DIMENSION( kts:kte ) :: qvsbar, rs0, viscmu, visc, diffwv, &
305 schmidt, xka
306 !---- new snow parameters
308 REAL, DIMENSION( kts:kte ):: ab_s,ab_r,ab_riming,lamc
309 REAL, DIMENSION( kts:kte ):: cap_s !---- capacitance of snow
311 REAL, PARAMETER :: vf1s = 0.65, vf2s = 0.44, vf1r =0.78 , vf2r = 0.31
313 REAL, PARAMETER :: am_c1=0.004, am_c2= 6e-5, am_c3=0.15
314 REAL, PARAMETER :: bm_c1=1.85, bm_c2= 0.003, bm_c3=1.25
315 REAL, PARAMETER :: aa_c1=1.28, aa_c2= -0.012, aa_c3=-0.6
316 REAL, PARAMETER :: ba_c1=1.5, ba_c2= 0.0075, ba_c3=0.5
318 REAL, PARAMETER :: best_a=1.08 , best_b = 0.499
319 REAL, DIMENSION(kts:kte):: am_s,bm_s,av_s,bv_s,Ri,N0_s,tmp_ss,lams
320 REAL, DIMENSION(kts:kte):: aa_s,ba_s,tmp_sa
321 REAL, PARAMETER :: mu_s=0.,mu_i=0.,mu_r=0.
323 REAL :: tc0, disp, Dc_liu, eta, mu_c, R6c !--- for Liu's autoconversion
325 ! Adding variable Riz, which will duplicate Ri but be a copy passed upward
327 REAL, DIMENSION(kts:kte) :: Riz
329 REAL, DIMENSION( kts:kte ) :: vtr, vts, &
330 vtrold, vtsold, vtiold, &
331 xlambdar, xlambdas, &
332 olambdar, olambdas
334 REAL :: episp0k, dtb, odtb, pi, pio4, &
335 pio6, oxLf, xLvocp, xLfocp, av_r, &
336 av_i, ocdrag, gambp4, gamdp4, &
337 gam4pt5, Cpor, oxmi, gambp3, gamdp3,&
338 gambp6, gam3pt5, gam2pt75, gambp5o2,&
339 gamdp5o2, cwoxlf, ocp, xni50, es
340 !
341 REAL :: qvmin=1.e-20
342 REAL :: temc1,save1,save2,xni50mx
344 ! for terminal velocity flux
346 INTEGER :: min_q, max_q, max_ri_k, k
347 REAL :: max_ri
348 REAL :: t_del_tv, del_tv, flux, fluxin, fluxout ,tmpqrz
349 LOGICAL :: notlast
350 !
351 mu_c = AMIN1(15., (1000.E6/Nt_c + 2.))
352 R6c = 10.0E-6 !---- 10 micron, threshold radius of cloud droplet
353 dtb=dt
354 odtb=1./dtb
355 pi =acos(-1.)
356 pio4=acos(-1.)/4.
357 pio6=acos(-1.)/6.
358 ocp=1./cp
359 oxLf=1./xLf
360 xLvocp=xLv/cp
361 xLfocp=xLf/cp
362 Cpor=cp/Rair
363 oxmi=1.0/xmi
364 cwoxlf=cw/xlf
365 av_r=2115.0*0.01**(1-bv_r)
366 av_i=152.93*0.01**(1-bv_i)
367 ocdrag=1./Cdrag
368 episp0k=RH*ep2*1000.*svp1
369 !
370 gambp4=ggamma(bv_r+4.)
371 gamdp4=ggamma(bv_i+4.)
372 gambp3=ggamma(bv_r+3.)
373 gambp6=ggamma(bv_r+6)
374 gambp5o2=ggamma((bv_r+5.)/2.)
375 gamdp5o2=ggamma((bv_i+5.)/2.)
376 !
377 ! oprez 1./prez ( prez : pressure)
378 ! qsw saturated mixing ratio on water surface
379 ! qsi saturated mixing ratio on ice surface
380 ! episp0k RH*e*saturated pressure at 273.15 K = 611.2 hPa (Rogers and Yau 1989)
381 ! qvoqsw qv/qsw
382 ! qvoqsi qv/qsi
383 ! qvzodt qv/dt
384 ! qlzodt ql/dt
385 ! qizodt qi/dt
386 ! qszodt qs/dt
387 ! qrzodt qr/dt
388 ! temcc temperature in dregee C
389 !
391 obp4=1.0/(bv_r+4.0)
392 bp3=bv_r+3.0
393 bp5=bv_r+5.0
394 bp6=bv_r+6.0
395 odp4=1.0/(bv_i+4.0)
396 dp3=bv_i+3.0
397 dp5=bv_i+5.0
398 dp5o2=0.5*(bv_i+5.0)
399 !
400 do k=kts,kte
401 oprez(k)=1./prez(k)
402 qlz(k)=amax1( 0.0,qlz(k) )
403 qiz(k)=amax1( 0.0,qiz(k) )
404 qvz(k)=amax1( qvmin,qvz(k) )
405 qsz(k)=amax1( 0.0,qsz(k) )
406 qrz(k)=amax1( 0.0,qrz(k) )
407 tem(k)=thz(k)*tothz(k)
408 temcc(k)=tem(k)-273.15
409 es=1000.*svp1*exp( svp2*temcc(k)/(tem(k)-svp3) ) !--- RY89 Eq(2.17)
410 qswz(k)=ep2*es/(prez(k)-es)
411 if (tem(k) .lt. 273.15 ) then
412 es=1000.*svp1*exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
413 qsiz(k)=ep2*es/(prez(k)-es)
414 if (temcc(k) .lt. -40.0) qswz(k)=qsiz(k)
415 else
416 qsiz(k)=qswz(k)
417 endif
418 !
419 qvoqswz(k)=qvz(k)/qswz(k)
420 qvoqsiz(k)=qvz(k)/qsiz(k)
421 qvzodt(k)=amax1( 0.0,odtb*qvz(k) )
422 qlzodt(k)=amax1( 0.0,odtb*qlz(k) )
423 qizodt(k)=amax1( 0.0,odtb*qiz(k) )
424 qszodt(k)=amax1( 0.0,odtb*qsz(k) )
425 qrzodt(k)=amax1( 0.0,odtb*qrz(k) )
426 theiz(k)=thz(k)+(xlvocp*qvz(k)-xlfocp*qiz(k))/tothz(k)
427 enddo
429 do k=kts,kte
431 psnow(k)=0.0
432 psaut(k)=0.0
433 psfw(k)=0.0
434 psfi(k)=0.0
435 praci(k)=0.0
436 piacr(k)=0.0
437 psaci(k)=0.0
438 psacw(k)=0.0
439 psdep(k)=0.0
440 pssub(k)=0.0
441 pracs(k)=0.0
442 psacr(k)=0.0
443 psmlt(k)=0.0
444 psmltevp(k)=0.0
446 prain(k)=0.0
447 praut(k)=0.0
448 pracw(k)=0.0
449 prevp(k)=0.0
450 pgfr(k)=0.0
452 pvapor(k)=0.0
454 pclw(k)=0.0
455 pladj(k)=0.0
457 pcli(k)=0.0
458 pimlt(k)=0.0
459 pihom(k)=0.0
460 pidw(k)=0.0
461 piadj(k)=0.0
463 qschg(k)=0.
465 enddo
467 !***********************************************************************
468 !***** compute viscosity,difusivity,thermal conductivity, and ******
469 !***** Schmidt number ******
470 !***********************************************************************
471 !c------------------------------------------------------------------
472 !c viscmu: dynamic viscosity of air kg/m/s
473 !c visc: kinematic viscosity of air = viscmu/rho (m2/s)
474 !c avisc=1.49628e-6 kg/m/s=1.49628e-5 g/cm/s
475 !c viscmu=1.718e-5 kg/m/s in RH
476 !c diffwv: Diffusivity of water vapor in air
477 !c adiffwv = 8.7602e-5 (8.794e-5 in MM5) kgm/s3
478 !c = 8.7602 (8.794 in MM5) gcm/s3
479 !c diffwv(k)=2.26e-5 m2/s
480 !c schmidt: Schmidt number=visc/diffwv
481 !c xka: thermal conductivity of air J/m/s/K (Kgm/s3/K)
482 !c xka(k)=2.43e-2 J/m/s/K in RH
483 !c axka=1.4132e3 (1.414e3 in MM5) m2/s2/k = 1.4132e7 cm2/s2/k
484 !c------------------------------------------------------------------
485 DO k=kts,kte
486 viscmu(k)=avisc*tem(k)**1.5/(tem(k)+120.0)
487 visc(k)=viscmu(k)*orho(k)
488 diffwv(k)=adiffwv*tem(k)**1.81*oprez(k)
489 schmidt(k)=visc(k)/diffwv(k)
490 xka(k)=axka*viscmu(k)
491 rs0(k)=ep2*1000.*svp1/(prez(k)-1000.*svp1)
492 END DO
493 !
494 ! ---- YLIN, set snow variables
495 !
496 !---- A+B in depositional growth, the first try just take from Rogers and Yau(1989)
497 ! ab_s(k) = lsub*lsub*orv/(tcond(k)*temp(k))+&
498 ! rv*temp(k)/(diffu(k)*qvsi(k))
500 do k = kts, kte
501 tc0 = tem(k)-273.15
502 if (rho(k)*qlz(k) .gt. 1e-5 .AND. rho(k)*qsz(k) .gt. 1e-5) then
503 Ri(k) = 1.0/(1.0+6e-5/(rho(k)**1.170*qlz(k)*qsz(k)**0.170))
504 else
505 Ri(k) = 0
506 endif
507 enddo
508 !
509 !--- make sure Ri does not decrease downward in a column
510 !
511 max_ri_k = MAXLOC(Ri,dim=1)
512 max_ri = MAXVAL(Ri)
514 do k = kts, max_ri_k
515 Ri(k) = max_ri
516 enddo
518 !--- YLIN, get PI properties
519 do k = kts, kte
520 Ri(k) = AMAX1(0.,AMIN1(Ri(k),1.0))
521 ! Store the value of Ri(k) as Riz(k)
522 Riz(k) = Ri(k)
524 cap_s(k)= 0.25*(1+Ri(k))
525 tc0 = AMIN1(-0.1, tem(k)-273.15)
526 N0_s(k) = amin1(2.0E8, 2.0E6*exp(-0.12*tc0))
527 am_s(k) = am_c1+am_c2*tc0+am_c3*Ri(k)*Ri(k) !--- Heymsfield 2007
528 am_s(k) = AMAX1(0.000023,am_s(k)) !--- use the a_min in table 1 of Heymsfield
529 bm_s(k) = bm_c1+bm_c2*tc0+bm_c3*Ri(k)
530 bm_s(k) = AMIN1(bm_s(k),3.0) !---- capped by 3
531 !--- converting from cgs to SI unit
532 am_s(k) = 10**(2*bm_s(k)-3.0)*am_s(k)
533 aa_s(k) = aa_c1 + aa_c2*tc0 + aa_c3*Ri(k)
534 ba_s(k) = ba_c1 + ba_c2*tc0 + ba_c3*Ri(k)
535 !--- convert to SI unit as in paper
536 aa_s(k) = (1e-2)**(2.0-ba_s(k))*aa_s(k)
537 !---- get v from Mitchell 1996
538 av_s(k) = best_a*viscmu(k)*(2*grav*am_s(k)/rho(k)/aa_s(k)/(viscmu(k)**2))**best_b
539 bv_s(k) = best_b*(bm_s(k)-ba_s(k)+2)-1
541 tmp_ss(k)= bm_s(k)+mu_s+1
542 tmp_sa(k)= ba_s(k)+mu_s+1
544 enddo
546 !
547 !***********************************************************************
548 ! Calculate precipitation fluxes due to terminal velocities.
549 !***********************************************************************
550 !
551 !- Calculate termianl velocity (vt?) of precipitation q?z
552 !- Find maximum vt? to determine the small delta t
553 !
554 !-- rain
555 !
556 ! CALL wrf_debug ( 100 , 'module_ylin, start precip fluxes' )
558 t_del_tv=0.
559 del_tv=dtb
560 notlast=.true.
562 DO while (notlast)
563 !
564 min_q=kte
565 max_q=kts-1
566 !
567 do k=kts,kte-1
568 if (qrz(k) .gt. 1.0e-8) then
569 min_q=min0(min_q,k)
570 max_q=max0(max_q,k)
571 tmp1=sqrt(pi*rhowater*xnor/rho(k)/qrz(k))
572 tmp1=sqrt(tmp1)
573 vtrold(k)=o6*av_r*gambp4*sqrho(k)/tmp1**bv_r
574 if (k .eq. 1) then
575 del_tv=amin1(del_tv,0.9*(zz(k)-zsfc)/vtrold(k))
576 else
577 del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtrold(k))
578 endif
579 else
580 vtrold(k)=0.
581 endif
582 enddo
584 if (max_q .ge. min_q) then
585 !
586 !- Check if the summation of the small delta t >= big delta t
587 ! (t_del_tv) (del_tv) (dtb)
589 t_del_tv=t_del_tv+del_tv
590 !
591 if ( t_del_tv .ge. dtb ) then
592 notlast=.false.
593 del_tv=dtb+del_tv-t_del_tv
594 endif
595 !
596 fluxin=0.
597 do k=max_q,min_q,-1
598 fluxout=rho(k)*vtrold(k)*qrz(k)
599 flux=(fluxin-fluxout)/rho(k)/dzw(k)
600 tmpqrz=qrz(k)
601 qrz(k)=qrz(k)+del_tv*flux
602 fluxin=fluxout
603 enddo
604 if (min_q .eq. 1) then
605 pptrain=pptrain+fluxin*del_tv
606 else
607 qrz(min_q-1)=qrz(min_q-1)+del_tv* &
608 fluxin/rho(min_q-1)/dzw(min_q-1)
609 endif
610 !
611 else
612 notlast=.false.
613 endif
614 ENDDO
616 !
617 !-- snow
618 !
619 t_del_tv=0.
620 del_tv=dtb
621 notlast=.true.
623 DO while (notlast)
624 !
625 min_q=kte
626 max_q=kts-1
627 !
628 do k=kts,kte-1
629 if (qsz(k) .gt. 1.0e-8) then
630 min_q=min0(min_q,k)
631 max_q=max0(max_q,k)
633 tmp1= (am_s(k)*N0_s(k)*ggamma(tmp_ss(k))*orho(k)/qsz(k))&
634 **(1./tmp_ss(k))
636 vtsold(k)= sqrho(k)*av_s(k)*ggamma(bv_s(k)+tmp_ss(k))/ &
637 ggamma(tmp_ss(k))/(tmp1**bv_s(k))
639 if (k .eq. 1) then
640 del_tv=amin1(del_tv,0.9*(zz(k)-zsfc)/vtsold(k))
641 else
642 del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtsold(k))
643 endif
644 else
645 vtsold(k)=0.
646 endif
647 enddo
649 if (max_q .ge. min_q) then
650 !
651 !
652 !- Check if the summation of the small delta t >= big delta t
653 ! (t_del_tv) (del_tv) (dtb)
655 t_del_tv=t_del_tv+del_tv
657 if ( t_del_tv .ge. dtb ) then
658 notlast=.false.
659 del_tv=dtb+del_tv-t_del_tv
660 endif
661 !
662 fluxin=0.
663 do k=max_q,min_q,-1
664 fluxout=rho(k)*vtsold(k)*qsz(k)
665 flux=(fluxin-fluxout)/rho(k)/dzw(k)
666 qsz(k)=qsz(k)+del_tv*flux
667 qsz(k)=amax1(0.,qsz(k))
668 fluxin=fluxout
669 enddo
670 if (min_q .eq. 1) then
671 pptsnow=pptsnow+fluxin*del_tv
672 else
673 qsz(min_q-1)=qsz(min_q-1)+del_tv* &
674 fluxin/rho(min_q-1)/dzw(min_q-1)
675 endif
676 !
677 else
678 notlast=.false.
679 endif
681 ENDDO
683 !
684 !-- cloud ice (03/21/02) using Heymsfield and Donner (1990) Vi=3.29*qi^0.16
685 !
686 t_del_tv=0.
687 del_tv=dtb
688 notlast=.true.
689 !
690 DO while (notlast)
691 !
692 min_q=kte
693 max_q=kts-1
694 !
695 do k=kts,kte-1
696 if (qiz(k) .gt. 1.0e-8) then
697 min_q=min0(min_q,k)
698 max_q=max0(max_q,k)
699 vtiold(k)= 3.29 * (rho(k)* qiz(k))** 0.16 ! Heymsfield and Donner
700 if (k .eq. 1) then
701 del_tv=amin1(del_tv,0.9*(zz(k)-zsfc)/vtiold(k))
702 else
703 del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtiold(k))
704 endif
705 else
706 vtiold(k)=0.
707 endif
708 enddo
710 if (max_q .ge. min_q) then
711 !
712 !- Check if the summation of the small delta t >= big delta t
713 ! (t_del_tv) (del_tv) (dtb)
715 t_del_tv=t_del_tv+del_tv
717 if ( t_del_tv .ge. dtb ) then
718 notlast=.false.
719 del_tv=dtb+del_tv-t_del_tv
720 endif
722 fluxin=0.
723 do k=max_q,min_q,-1
724 fluxout=rho(k)*vtiold(k)*qiz(k)
725 flux=(fluxin-fluxout)/rho(k)/dzw(k)
726 qiz(k)=qiz(k)+del_tv*flux
727 qiz(k)=amax1(0.,qiz(k))
728 fluxin=fluxout
729 enddo
730 if (min_q .eq. 1) then
731 pptice=pptice+fluxin*del_tv
732 else
733 qiz(min_q-1)=qiz(min_q-1)+del_tv* &
734 fluxin/rho(min_q-1)/dzw(min_q-1)
735 endif
736 !
737 else
738 notlast=.false.
739 endif
740 !
741 ENDDO
743 ! CALL wrf_debug ( 100 , 'module_ylin: end precip flux' )
745 ! Microphpysics processes
747 DO 2000 k=kts,kte
748 !
749 qvzodt(k)=amax1( 0.0,odtb*qvz(k) )
750 qlzodt(k)=amax1( 0.0,odtb*qlz(k) )
751 qizodt(k)=amax1( 0.0,odtb*qiz(k) )
752 qszodt(k)=amax1( 0.0,odtb*qsz(k) )
753 qrzodt(k)=amax1( 0.0,odtb*qrz(k) )
755 !***********************************************************************
756 !***** diagnose mixing ratios (qrz,qsz), terminal *****
757 !***** velocities (vtr,vts), and slope parameters in size *****
758 !***** distribution(xlambdar,xlambdas) of rain and snow *****
759 !***** follows Nagata and Ogura, 1991, MWR, 1309-1337. Eq (A7) *****
760 !***********************************************************************
761 !
762 !**** assuming no cloud water can exist in the top two levels due to
763 !**** radiation consideration
764 !
765 !! if
766 !! unsaturated,
767 !! no cloud water, rain, ice, snow
768 !! then
769 !! skip these processes and jump to line 2000
770 !
771 !
772 tmp=qiz(k)+qlz(k)+qsz(k)+qrz(k)
773 if( qvz(k)+qlz(k)+qiz(k) .lt. qsiz(k) &
774 .and. tmp .eq. 0.0 ) go to 2000
775 !
776 !! calculate terminal velocity of rain
777 !
778 if (qrz(k) .gt. 1.0e-8) then
779 tmp1=sqrt(pi*rhowater*xnor*orho(k)/qrz(k))
780 xlambdar(k)=sqrt(tmp1)
781 olambdar(k)=1.0/xlambdar(k)
782 vtrold(k)=o6*av_r*gambp4*sqrho(k)*olambdar(k)**bv_r
783 else
784 vtrold(k)=0.
785 olambdar(k)=0.
786 endif
787 !
788 if (qrz(k) .gt. 1.0e-8) then
789 tmp1=sqrt(pi*rhowater*xnor*orho(k)/qrz(k))
790 xlambdar(k)=sqrt(tmp1)
791 olambdar(k)=1.0/xlambdar(k)
792 vtr(k)=o6*av_r*gambp4*sqrho(k)*olambdar(k)**bv_r
793 else
794 vtr(k)=0.
795 olambdar(k)=0.
796 endif
797 !
798 !! calculate terminal velocity of snow
799 !
800 if (qsz(k) .gt. 1.0e-8) then
801 tmp1= (am_s(k)*N0_s(k)*ggamma(tmp_ss(k))*orho(k)/qsz(k))&
802 **(1./tmp_ss(k))
803 olambdas(k)=1.0/tmp1
804 vtsold(k)= sqrho(k)*av_s(k)*ggamma(bv_s(k)+tmp_ss(k))/ &
805 ggamma(tmp_ss(k))/(tmp1**bv_s(k))
807 else
808 vtsold(k)=0.
809 olambdas(k)=0.
810 endif
811 !
812 if (qsz(k) .gt. 1.0e-8) then
813 tmp1= (am_s(k)*N0_s(k)*ggamma(tmp_ss(k))*orho(k)/qsz(k))&
814 **(1./tmp_ss(k))
815 olambdas(k)=1.0/tmp1
816 vts(k)= sqrho(k)*av_s(k)*ggamma(bv_s(k)+tmp_ss(k))/ &
817 ggamma(tmp_ss(k))/(tmp1**bv_s(k))
819 else
820 vts(k)=0.
821 olambdas(k)=0.
822 endif
824 !---------- start of snow/ice processes below freezing
826 if (tem(k) .lt. 273.15) then
828 !
829 !***********************************************************************
830 !********* snow production processes for T < 0 C **********
831 !***********************************************************************
832 !c
833 !c (1) ICE CRYSTAL AGGREGATION TO SNOW (Psaut): Lin (21)
834 !c! psaut=alpha1*(qi-qi0)
835 !c! alpha1=1.0e-3*exp(0.025*(T-T0))
836 !c
837 alpha1=1.0e-3*exp( 0.025*temcc(k) )
838 !
839 if(temcc(k) .lt. -20.0) then
840 tmp1=-7.6+4.0*exp( -0.2443e-3*(abs(temcc(k))-20)**2.455 )
841 qic=1.0e-3*exp(tmp1)*orho(k)
842 else
843 qic=qi0
844 end if
846 tmp1=odtb*(qiz(k)-qic)*(1.0-exp(-alpha1*dtb))
847 psaut(k)=amax1( 0.0,tmp1 )
849 !c
850 !c (2) BERGERON PROCESS TRANSFER OF CLOUD WATER TO SNOW (Psfw)
851 !c this process only considered when -31 C < T < 0 C
852 !c Lin (33) and Hsie (17)
853 !c
854 !c!
855 !c! parama1 and parama2 functions must be user supplied
856 !c!
858 if( qlz(k) .gt. 1.0e-10 ) then
859 temc1=amax1(-30.99,temcc(k))
860 a1=parama1( temc1 )
861 a2=parama2( temc1 )
862 tmp1=1.0-a2
863 !! change unit from cgs to mks
864 a1=a1*0.001**tmp1
865 !! dtberg is the time needed for a crystal to grow from 40 to 50 um
866 !! odtberg=1.0/dtberg
867 odtberg=(a1*tmp1)/(xmi50**tmp1-xmi40**tmp1)
868 !
869 !! compute terminal velocity of a 50 micron ice cystal
870 !
871 vti50=av_i*di50**bv_i*sqrho(k)
872 !
873 eiw=1.0
874 save1=a1*xmi50**a2
875 save2=0.25*pi*eiw*rho(k)*di50*di50*vti50
876 !
877 tmp2=( save1 + save2*qlz(k) )
878 !
879 !! maximum number of 50 micron crystals limited by the amount
880 !! of supercool water
881 !
882 xni50mx=qlzodt(k)/tmp2
883 !
884 !! number of 50 micron crystals produced
885 !
886 xni50=qiz(k)*( 1.0-exp(-dtb*odtberg) )/xmi50
887 xni50=amin1(xni50,xni50mx)
888 !
889 tmp3=odtb*tmp2/save2*( 1.0-exp(-save2*xni50*dtb) )
890 psfw(k)=amin1( tmp3,qlzodt(k) )
891 !c
892 !c (3) REDUCTION OF CLOUD ICE BY BERGERON PROCESS (Psfi): Lin (34)
893 !c this process only considered when -31 C < T < 0 C
894 !c
895 tmp1=xni50*xmi50-psfw(k)
896 psfi(k)=amin1(tmp1,qizodt(k))
897 end if
898 !
899 !
900 if(qrz(k) .le. 0.0) go to 1000
901 !
902 ! Processes (4) and (5) only need when qrz > 0.0
903 !
904 !c
905 !c (4) CLOUD ICE ACCRETION BY RAIN (Praci): Lin (25)
906 !c produce PI
907 !c
908 eri=1.0
909 save1=pio4*eri*xnor*av_r*sqrho(k)
910 tmp1=save1*gambp3*olambdar(k)**bp3
911 praci(k)=qizodt(k)*( 1.0-exp(-tmp1*dtb) )
913 !c
914 !c (5) RAIN ACCRETION BY CLOUD ICE (Piacr): Lin (26)
915 !c
916 tmp2=qiz(k)*save1*rho(k)*pio6*rhowater*gambp6*oxmi* &
917 olambdar(k)**bp6
918 piacr(k)=amin1( tmp2,qrzodt(k) )
920 !
921 1000 continue
922 !
923 if(qsz(k) .le. 0.0) go to 1200
924 !
925 ! Compute the following processes only when qsz > 0.0
926 !
927 !c
928 !c (6) ICE CRYSTAL ACCRETION BY SNOW (Psaci): Lin (22)
929 !c
930 esi=exp( 0.025*temcc(k) )
931 save1 = aa_s(k)*sqrho(k)*N0_s(k)* &
932 ggamma(bv_s(k)+tmp_sa(k))*olambdas(k)**(bv_s(k)+tmp_sa(k))
934 tmp1=esi*save1
935 psaci(k)=qizodt(k)*( 1.0-exp(-tmp1*dtb) )
937 !c
938 !c (7) CLOUD WATER ACCRETION BY SNOW (Psacw): Lin (24)
939 !c
940 esw=1.0
941 tmp1=esw*save1
942 psacw(k)=qlzodt(K)*( 1.0-exp(-tmp1*dtb) )
944 !c
945 !c (8) DEPOSITION/SUBLIMATION OF SNOW (Psdep/Pssub): Lin (31)
946 !c includes consideration of ventilation effect
947 !c
948 tmpa=rvapor*xka(k)*tem(k)*tem(k)
949 tmpb=xls*xls*rho(k)*qsiz(k)*diffwv(k)
950 tmpc=tmpa*qsiz(k)*diffwv(k)
951 abi=4.0*pi*cap_s(k)*(qvoqsiz(k)-1.0)*tmpc/(tmpa+tmpb)
952 tmp1=av_s(k)*sqrho(k)*olambdas(k)**(5+bv_s(k)+2*mu_s)/visc(k)
954 !---- YLIN, here there is some approximation assuming mu_s =1, so gamma(2)=1, etc.
956 tmp2= abi*N0_s(k)*( vf1s*olambdas(k)*olambdas(k)+ &
957 vf2s*schmidt(k)**0.33334* &
958 ggamma(2.5+0.5*bv_s(k)+mu_s)*sqrt(tmp1) )
960 tmp3=odtb*( qvz(k)-qsiz(k) )
961 !
962 if( tmp2 .le. 0.0) then
963 tmp2=amax1( tmp2,tmp3)
964 pssub(k)=amax1( tmp2,-qszodt(k) )
965 psdep(k)=0.0
966 else
967 psdep(k)=amin1( tmp2,tmp3 )
968 pssub(k)=0.0
969 end if
971 !
972 if(qrz(k) .le. 0.0) go to 1200
973 !
974 ! Compute processes (9) and (10) only when qsz > 0.0 and qrz > 0.0
975 ! these two terms need to be refined in the future, they should be equal
976 !c
977 !c (9) ACCRETION OF SNOW BY RAIN (Pracs): Lin (27)
978 !c
979 esr=1.0
980 tmpa=olambdar(k)*olambdar(k)
981 tmpb=olambdas(k)*olambdas(k)
982 tmpc=olambdar(k)*olambdas(k)
983 tmp1=pi*pi*esr*xnor*N0_s(k)*abs( vtr(k)-vts(k) )*orho(k)
984 tmp2=tmpb*tmpb*olambdar(k)*(5.0*tmpb+2.0*tmpc+0.5*tmpa)
985 tmp3=tmp1*rhosnow*tmp2
986 pracs(k)=amin1( tmp3,qszodt(k) )
987 pracs(k)=0.0
988 !c
989 !c (10) ACCRETION OF RAIN BY SNOW (Psacr): Lin (28)
990 !c
991 tmp3=tmpa*tmpa*olambdas(k)*(5.0*tmpa+2.0*tmpc+0.5*tmpb)
992 tmp4=tmp1*rhowater*tmp3
993 psacr(k)=amin1( tmp4,qrzodt(k) )
994 !
995 !c
996 !c (2) FREEZING OF RAIN TO FORM GRAUPEL (pgfr): Lin (45), added to PI
997 !c positive value
998 !c Constant in Bigg freezing Aplume=Ap=0.66 /k
999 !c Constant in raindrop freezing equ. Bplume=Bp=100./m/m/m/s
1000 !
1002 if (qrz(k) .gt. 1.e-8 ) then
1003 Bp=100.
1004 Ap=0.66
1005 tmp1=olambdar(k)*olambdar(k)*olambdar(k)
1006 tmp2=20.*pi*pi*Bp*xnor*rhowater*orho(k)* &
1007 (exp(-Ap*temcc(k))-1.0)*tmp1*tmp1*olambdar(k)
1008 pgfr(k)=amin1( tmp2,qrzodt(k) )
1009 else
1010 pgfr(k)=0
1011 endif
1013 1200 continue
1014 !
1016 else
1018 !
1019 !***********************************************************************
1020 !********* snow production processes for T > 0 C **********
1021 !***********************************************************************
1022 !
1023 if (qsz(k) .le. 0.0) go to 1400
1024 !c
1025 !c (1) CLOUD WATER ACCRETION BY SNOW (Psacw): Lin (24)
1026 !c
1027 esw=1.0
1029 save1 =aa_s(k)*sqrho(k)*N0_s(k)* &
1030 ggamma(bv_s(k)+tmp_sa(k))*olambdas(k)**(bv_s(k)+tmp_sa(k))
1032 tmp1=esw*save1
1033 psacw(k)=qlzodt(k)*( 1.0-exp(-tmp1*dtb) )
1035 !c
1036 !c (2) ACCRETION OF RAIN BY SNOW (Psacr): Lin (28)
1037 !c
1038 esr=1.0
1039 tmpa=olambdar(k)*olambdar(k)
1040 tmpb=olambdas(k)*olambdas(k)
1041 tmpc=olambdar(k)*olambdas(k)
1042 tmp1=pi*pi*esr*xnor*N0_s(k)*abs( vtr(k)-vts(k) )*orho(k)
1043 tmp2=tmpa*tmpa*olambdas(k)*(5.0*tmpa+2.0*tmpc+0.5*tmpb)
1044 tmp3=tmp1*rhowater*tmp2
1045 psacr(k)=amin1( tmp3,qrzodt(k) )
1046 !c
1047 !c (3) MELTING OF SNOW (Psmlt): Lin (32)
1048 !c Psmlt is negative value
1049 !
1050 delrs=rs0(k)-qvz(k)
1051 term1=2.0*pi*orho(k)*( xlv*diffwv(k)*rho(k)*delrs- &
1052 xka(k)*temcc(k) )
1053 tmp1= av_s(k)*sqrho(k)*olambdas(k)**(5+bv_s(k)+2*mu_s)/visc(k)
1054 tmp2= N0_s(k)*( vf1s*olambdas(k)*olambdas(k)+ &
1055 vf2s*schmidt(k)**0.33334* &
1056 ggamma(2.5+0.5*bv_s(k)+mu_s)*sqrt(tmp1) )
1057 tmp3=term1*oxlf*tmp2-cwoxlf*temcc(k)*( psacw(k)+psacr(k) )
1058 tmp4=amin1(0.0,tmp3)
1059 psmlt(k)=amax1( tmp4,-qszodt(k) )
1060 !c
1061 !c (4) EVAPORATION OF MELTING SNOW (Psmltevp): HR (A27)
1062 !c but use Lin et al. coefficience
1063 !c Psmltevp is a negative value
1064 !c
1065 tmpa=rvapor*xka(k)*tem(k)*tem(k)
1066 tmpb=xlv*xlv*rho(k)*qswz(k)*diffwv(k)
1067 tmpc=tmpa*qswz(k)*diffwv(k)
1068 tmpd=amin1( 0.0,(qvoqswz(k)-0.90)*qswz(k)*odtb )
1070 abr=2.0*pi*(qvoqswz(k)-0.90)*tmpc/(tmpa+tmpb)
1071 !
1072 !**** allow evaporation to occur when RH less than 90%
1073 !**** here not using 100% because the evaporation cooling
1074 !**** of temperature is not taking into account yet; hence,
1075 !**** the qsw value is a little bit larger. This will avoid
1076 !**** evaporation can generate cloud.
1077 !
1078 tmp1=av_s(k)*sqrho(k)*olambdas(k)**(5+bv_s(k)+2*mu_s)/visc(k)
1079 tmp2= N0_s(k)*( vf1s*olambdas(k)*olambdas(k)+ &
1080 vf2s*schmidt(k)**0.33334* &
1081 ggamma(2.5+0.5*bv_s(k)+mu_s)*sqrt(tmp1) )
1082 tmp3=amin1(0.0,tmp2)
1083 tmp3=amax1( tmp3,tmpd )
1084 psmltevp(k)=amax1( tmp3,-qszodt(k) )
1085 1400 continue
1086 !
1087 end if !---- end of snow/ice processes
1089 ! CALL wrf_debug ( 100 , 'module_ylin: finish ice/snow processes' )
1092 !***********************************************************************
1093 !********* rain production processes **********
1094 !***********************************************************************
1096 !c
1097 !c (1) AUTOCONVERSION OF RAIN (Praut): using Liu and Daum (2004)
1098 !c
1100 !---- YLIN, autoconversion use Liu and Daum (2004), unit = g cm-3 s-1, in the scheme kg/kg s-1, so
1102 if (qlz(k) .gt. 1e-6) then
1103 lamc(k) = (Nt_c*rhowater*pi*ggamma(4.+mu_c)/(6.*rho(k)*qlz(k))/ & !--- N(D) = N0*D^mu*exp(-lamc*D)
1104 ggamma(1+mu_c))**0.3333
1105 Dc_liu = (ggamma(6+1+mu_c)/ggamma(1+mu_c))**(1./6.)/lamc(k) !----- R6 in m
1107 if (Dc_liu .gt. R6c) then
1108 disp = 1./(mu_c+1.) !--- square of relative dispersion
1109 eta = (0.75/pi/(1e-3*rhowater))**2*1.9e11*((1+3*disp)*(1+4*disp)*&
1110 (1+5*disp)/(1+disp)/(1+2*disp))
1111 praut(k) = eta*(1e-3*rho(k)*qlz(k))**3/(1e-6*Nt_c) !--- g cm-3 s-1
1112 praut(k) = praut(k)/(1e-3*rho(k)) !--- kg kg-1 s-1
1113 else
1114 praut(k) = 0.0
1115 endif
1116 else
1117 praut(k) = 0.0
1118 endif
1120 !c
1121 !c (2) ACCRETION OF CLOUD WATER BY RAIN (Pracw): Lin (51)
1122 !c
1123 erw=1.0
1125 tmp1=pio4*erw*xnor*av_r*sqrho(k)* &
1126 gambp3*olambdar(k)**bp3
1127 pracw(k)=qlzodt(k)*( 1.0-exp(-tmp1*dtb) )
1129 !c
1130 !c (3) EVAPORATION OF RAIN (Prevp): Lin (52)
1131 !c Prevp is negative value
1132 !c
1133 !c Sw=qvoqsw : saturation ratio
1134 !c
1135 tmpa=rvapor*xka(k)*tem(k)*tem(k)
1136 tmpb=xlv*xlv*rho(k)*qswz(k)*diffwv(k)
1137 tmpc=tmpa*qswz(k)*diffwv(k)
1138 tmpd=amin1(0.0,(qvoqswz(k)-0.90)*qswz(k)*odtb)
1140 abr=2.0*pi*(qvoqswz(k)-0.90)*tmpc/(tmpa+tmpb)
1141 tmp1=av_r*sqrho(k)*olambdar(k)**bp5/visc(k)
1142 tmp2=abr*xnor*( vf1r*olambdar(k)*olambdar(k)+ &
1143 vf2r*schmidt(k)**0.33334*gambp5o2*sqrt(tmp1) )
1144 tmp3=amin1( 0.0,tmp2 )
1145 tmp3=amax1( tmp3,tmpd )
1146 prevp(k)=amax1( tmp3,-qrzodt(k) )
1148 ! CALL wrf_debug ( 100 , 'module_ylin: finish rain processes' )
1150 !c
1151 !c**********************************************************************
1152 !c***** combine all processes together and avoid negative *****
1153 !c***** water substances
1154 !***********************************************************************
1155 !c
1156 if ( temcc(k) .lt. 0.0) then
1157 !c
1158 !c combined water vapor depletions
1159 !c
1160 tmp=psdep(k)
1161 if ( tmp .gt. qvzodt(k) ) then
1162 factor=qvzodt(k)/tmp
1163 psdep(k)=psdep(k)*factor
1164 end if
1165 !c
1166 !c combined cloud water depletions
1167 !c
1168 tmp=praut(k)+psacw(k)+psfw(k)+pracw(k)
1169 if ( tmp .gt. qlzodt(k) ) then
1170 factor=qlzodt(k)/tmp
1171 praut(k)=praut(k)*factor
1172 psacw(k)=psacw(k)*factor
1173 psfw(k)=psfw(k)*factor
1174 pracw(k)=pracw(k)*factor
1175 end if
1176 !c
1177 !c combined cloud ice depletions
1178 !c
1179 tmp=psaut(k)+psaci(k)+praci(k)+psfi(k)
1180 if (tmp .gt. qizodt(k) ) then
1181 factor=qizodt(k)/tmp
1182 psaut(k)=psaut(k)*factor
1183 psaci(k)=psaci(k)*factor
1184 praci(k)=praci(k)*factor
1185 psfi(k)=psfi(k)*factor
1186 endif
1187 !c
1188 !c combined all rain processes
1189 !c
1190 tmp_r=piacr(k)+psacr(k)-prevp(k)-praut(k)-pracw(k)+pgfr(k)
1191 if (tmp_r .gt. qrzodt(k) ) then
1192 factor=qrzodt(k)/tmp_r
1193 piacr(k)=piacr(k)*factor
1194 psacr(k)=psacr(k)*factor
1195 prevp(k)=prevp(k)*factor
1196 pgfr(k)=pgfr(k)*factor
1197 endif
1198 !c
1199 !c combined all snow processes
1200 !c
1201 tmp_s=-pssub(k)-(psaut(k)+psaci(k)+psacw(k)+psfw(k)+pgfr(k)+ &
1202 psfi(k)+praci(k)+piacr(k)+ &
1203 psdep(k)+psacr(k)-pracs(k))
1204 if ( tmp_s .gt. qszodt(k) ) then
1205 factor=qszodt(k)/tmp_s
1206 pssub(k)=pssub(k)*factor
1207 Pracs(k)=Pracs(k)*factor
1208 endif
1210 !c
1211 !c calculate new water substances, thetae, tem, and qvsbar
1212 !c
1214 pvapor(k)=-pssub(k)-psdep(k)-prevp(k)
1215 qvz(k)=amax1( qvmin,qvz(k)+dtb*pvapor(k) )
1216 pclw(k)=-praut(k)-pracw(k)-psacw(k)-psfw(k)
1217 qlz(k)=amax1( 0.0,qlz(k)+dtb*pclw(k) )
1218 pcli(k)=-psaut(k)-psfi(k)-psaci(k)-praci(k)
1219 qiz(k)=amax1( 0.0,qiz(k)+dtb*pcli(k) )
1220 tmp_r=piacr(k)+psacr(k)-prevp(k)-praut(k)-pracw(k)+pgfr(k)-pracs(k)
1221 prain(k)=-tmp_r
1222 qrz(k)=amax1( 0.0,qrz(k)+dtb*prain(k) )
1223 tmp_s=-pssub(k)-(psaut(k)+psaci(k)+psacw(k)+psfw(k)+pgfr(k)+ &
1224 psfi(k)+praci(k)+piacr(k)+ &
1225 psdep(k)+psacr(k)-pracs(k))
1226 psnow(k)=-tmp_s
1227 qsz(k)=amax1( 0.0,qsz(k)+dtb*psnow(k) )
1229 qschg(k)=qschg(k)+psnow(k)
1230 qschg(k)=psnow(k)
1232 tmp=ocp/tothz(k)*xLf*qschg(k)
1233 theiz(k)=theiz(k)+dtb*tmp
1234 thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
1235 tem(k)=thz(k)*tothz(k)
1237 temcc(k)=tem(k)-273.15
1239 if( temcc(k) .lt. -40.0 ) qswz(k)=qsiz(k)
1240 qlpqi=qlz(k)+qiz(k)
1241 if ( qlpqi .eq. 0.0 ) then
1242 qvsbar(k)=qsiz(k)
1243 else
1244 qvsbar(k)=( qiz(k)*qsiz(k)+qlz(k)*qswz(k) )/qlpqi
1245 endif
1247 !
1248 else !>0 C
1249 !c
1250 !c combined cloud water depletions
1251 !c
1252 tmp=praut(k)+psacw(k)+pracw(k)
1253 if ( tmp .gt. qlzodt(k) ) then
1254 factor=qlzodt(k)/tmp
1255 praut(k)=praut(k)*factor
1256 psacw(k)=psacw(k)*factor
1257 pracw(k)=pracw(k)*factor
1258 end if
1259 !c
1260 !c combined all snow processes
1261 !c
1262 tmp_s=-(psmlt(k)+psmltevp(k))
1263 if (tmp_s .gt. qszodt(k) ) then
1264 factor=qszodt(k)/tmp_s
1265 psmlt(k)=psmlt(k)*factor
1266 psmltevp(k)=psmltevp(k)*factor
1267 endif
1268 !c
1269 !c combined all rain processes
1270 !c
1271 tmp_r=-prevp(k)-(praut(k)+pracw(k)+psacw(k)-psmlt(k))
1272 if (tmp_r .gt. qrzodt(k) ) then
1273 factor=qrzodt(k)/tmp_r
1274 prevp(k)=prevp(k)*factor
1275 endif
1276 !c
1277 !c calculate new water substances and thetae
1278 !c
1279 pvapor(k)=-psmltevp(k)-prevp(k)
1280 qvz(k)=amax1( qvmin,qvz(k)+dtb*pvapor(k))
1281 pclw(k)=-praut(k)-pracw(k)-psacw(k)
1282 qlz(k)=amax1( 0.0,qlz(k)+dtb*pclw(k) )
1283 pcli(k)=0.0
1284 qiz(k)=amax1( 0.0,qiz(k)+dtb*pcli(k) )
1285 tmp_r=-prevp(k)-(praut(k)+pracw(k)+psacw(k)-psmlt(k))
1286 prain(k)=-tmp_r
1287 tmpqrz=qrz(k)
1288 qrz(k)=amax1( 0.0,qrz(k)+dtb*prain(k) )
1289 tmp_s=-(psmlt(k)+psmltevp(k))
1290 psnow(k)=-tmp_s
1291 qsz(k)=amax1( 0.0,qsz(k)+dtb*psnow(k) )
1292 qschg(k)=psnow(k)
1293 !
1294 tmp=ocp/tothz(k)*xLf*qschg(k)
1295 theiz(k)=theiz(k)+dtb*tmp
1296 thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
1298 tem(k)=thz(k)*tothz(k)
1299 temcc(k)=tem(k)-273.15
1300 es=1000.*svp1*exp( svp2*temcc(k)/(tem(k)-svp3) )
1301 qswz(k)=ep2*es/(prez(k)-es)
1302 qsiz(k)=qswz(k)
1303 qvsbar(k)=qswz(k)
1304 !
1305 end if
1306 ! CALL wrf_debug ( 100 , 'module_ylin: finish sum of all processes' )
1308 !
1309 !***********************************************************************
1310 !********** saturation adjustment **********
1311 !***********************************************************************
1312 !
1313 ! allow supersaturation exits linearly from 0% at 500 mb to 50%
1314 ! above 300 mb
1315 ! 5.0e-5=1.0/(500mb-300mb)
1316 !
1317 rsat=1.0
1318 if( qvz(k)+qlz(k)+qiz(k) .lt. rsat*qvsbar(k) ) then
1320 !c
1321 !c unsaturated
1322 !c
1323 qvz(k)=qvz(k)+qlz(k)+qiz(k)
1324 qlz(k)=0.0
1325 qiz(k)=0.0
1327 thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
1328 tem(k)=thz(k)*tothz(k)
1329 temcc(k)=tem(k)-273.15
1331 go to 1800
1332 !
1333 else
1334 !c
1335 !c saturated
1336 !c
1337 pladj(k)=qlz(k)
1338 piadj(k)=qiz(k)
1339 !
1341 CALL satadj(qvz, qlz, qiz, prez, theiz, thz, tothz, kts, kte, &
1342 k, xLvocp, xLfocp, episp0k, EP2,SVP1,SVP2,SVP3,SVPT0)
1344 !
1345 pladj(k)=odtb*(qlz(k)-pladj(k))
1346 piadj(k)=odtb*(qiz(k)-piadj(k))
1347 !
1348 pclw(k)=pclw(k)+pladj(k)
1349 pcli(k)=pcli(k)+piadj(k)
1350 pvapor(k)=pvapor(k)-( pladj(k)+piadj(k) )
1351 !
1352 thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
1353 tem(k)=thz(k)*tothz(k)
1355 temcc(k)=tem(k)-273.15
1357 es=1000.*svp1*exp( svp2*temcc(k)/(tem(k)-svp3) )
1358 qswz(k)=ep2*es/(prez(k)-es)
1359 if (tem(k) .lt. 273.15 ) then
1360 es=1000.*svp1*exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
1361 qsiz(k)=ep2*es/(prez(k)-es)
1362 if (temcc(k) .lt. -40.0) qswz(k)=qsiz(k)
1363 else
1364 qsiz(k)=qswz(k)
1365 endif
1366 qlpqi=qlz(k)+qiz(k)
1367 if ( qlpqi .eq. 0.0 ) then
1368 qvsbar(k)=qsiz(k)
1369 else
1370 qvsbar(k)=( qiz(k)*qsiz(k)+qlz(k)*qswz(k) )/qlpqi
1371 endif
1373 end if
1375 !
1376 !***********************************************************************
1377 !***** melting and freezing of cloud ice and cloud water *****
1378 !***********************************************************************
1379 qlpqi=qlz(k)+qiz(k)
1380 if(qlpqi .le. 0.0) go to 1800
1381 !
1382 !c
1383 !c (1) HOMOGENEOUS NUCLEATION WHEN T< -40 C (Pihom)
1384 !c
1385 if(temcc(k) .lt. -40.0) pihom(k)=qlz(k)*odtb
1386 !c
1387 !c (2) MELTING OF ICE CRYSTAL WHEN T> 0 C (Pimlt)
1388 !c
1389 if(temcc(k) .gt. 0.0) pimlt(k)=qiz(k)*odtb
1390 !c
1391 !c (3) PRODUCTION OF CLOUD ICE BY BERGERON PROCESS (Pidw): Hsie (p957)
1392 !c this process only considered when -31 C < T < 0 C
1393 !c
1394 if(temcc(k) .lt. 0.0 .and. temcc(k) .gt. -31.0) then
1395 !c!
1396 !c! parama1 and parama2 functions must be user supplied
1397 !c!
1398 a1=parama1( temcc(k) )
1399 a2=parama2( temcc(k) )
1400 !! change unit from cgs to mks
1401 a1=a1*0.001**(1.0-a2)
1402 xnin=xni0*exp(-bni*temcc(k))
1403 pidw(k)=xnin*orho(k)*(a1*xmnin**a2)
1404 end if
1405 !
1406 pcli(k)=pcli(k)+pihom(k)-pimlt(k)+pidw(k)
1407 pclw(k)=pclw(k)-pihom(k)+pimlt(k)-pidw(k)
1408 qlz(k)=amax1( 0.0,qlz(k)+dtb*(-pihom(k)+pimlt(k)-pidw(k)) )
1409 qiz(k)=amax1( 0.0,qiz(k)+dtb*(pihom(k)-pimlt(k)+pidw(k)) )
1411 !
1412 CALL satadj(qvz, qlz, qiz, prez, theiz, thz, tothz, kts, kte, &
1413 k, xLvocp, xLfocp, episp0k ,EP2,SVP1,SVP2,SVP3,SVPT0)
1415 thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
1416 tem(k)=thz(k)*tothz(k)
1418 temcc(k)=tem(k)-273.15
1420 es=1000.*svp1*exp( svp2*temcc(k)/(tem(k)-svp3) )
1421 qswz(k)=ep2*es/(prez(k)-es)
1423 if (tem(k) .lt. 273.15 ) then
1424 es=1000.*svp1*exp( 21.8745584*(tem(k)-273.16)/(tem(k)-7.66) )
1425 qsiz(k)=ep2*es/(prez(k)-es)
1426 if (temcc(k) .lt. -40.0) qswz(k)=qsiz(k)
1427 else
1428 qsiz(k)=qswz(k)
1429 endif
1430 qlpqi=qlz(k)+qiz(k)
1432 if ( qlpqi .eq. 0.0 ) then
1433 qvsbar(k)=qsiz(k)
1434 else
1435 qvsbar(k)=( qiz(k)*qsiz(k)+qlz(k)*qswz(k) )/qlpqi
1436 endif
1438 1800 continue
1439 !
1440 !***********************************************************************
1441 !********** integrate the productions of rain and snow **********
1442 !***********************************************************************
1443 !
1444 2000 continue
1446 !
1447 !**** below if qv < qvmin then qv=qvmin, ql=0.0, and qi=0.0
1448 !
1449 do k=kts+1,kte
1450 if ( qvz(k) .lt. qvmin ) then
1451 qlz(k)=0.0
1452 qiz(k)=0.0
1453 qvz(k)=amax1( qvmin,qvz(k)+qlz(k)+qiz(k) )
1454 end if
1455 enddo
1456 !
1458 ! CALL wrf_debug ( 100 , 'module_ylin: finish saturation adjustment' )
1460 END SUBROUTINE clphy1d_ylin
1465 !---------------------------------------------------------------------
1466 ! SATURATED ADJUSTMENT
1467 !---------------------------------------------------------------------
1468 SUBROUTINE satadj(qvz, qlz, qiz, prez, theiz, thz, tothz, &
1469 kts, kte, k, xLvocp, xLfocp, episp0k, EP2,SVP1,SVP2,SVP3,SVPT0)
1470 !---------------------------------------------------------------------
1471 IMPLICIT NONE
1472 !---------------------------------------------------------------------
1473 ! This program use Newton's method for finding saturated temperature
1474 ! and saturation mixing ratio.
1475 !
1476 ! In this saturation adjustment scheme we assume
1477 ! (1) the saturation mixing ratio is the mass weighted average of
1478 ! saturation values over liquid water (qsw), and ice (qsi)
1479 ! following Lord et al., 1984 and Tao, 1989
1480 !
1481 ! (2) the percentage of cloud liquid and cloud ice will
1482 ! be fixed during the saturation calculation
1483 !---------------------------------------------------------------------
1484 !
1486 INTEGER, INTENT(IN ) :: kts, kte, k
1488 REAL, DIMENSION( kts:kte ), &
1489 INTENT(INOUT) :: qvz, qlz, qiz
1490 !
1491 REAL, DIMENSION( kts:kte ), &
1492 INTENT(IN ) :: prez, theiz, tothz
1494 REAL, INTENT(IN ) :: xLvocp, xLfocp, episp0k
1495 REAL, INTENT(IN ) :: EP2,SVP1,SVP2,SVP3,SVPT0
1497 ! LOCAL VARS
1499 INTEGER :: n
1501 REAL, DIMENSION( kts:kte ) :: thz, tem, temcc, qsiz, &
1502 qswz, qvsbar
1504 REAL :: qsat, qlpqi, ratql, t0, t1, tmp1, ratqi, tsat, absft, &
1505 denom1, denom2, dqvsbar, ftsat, dftsat, qpz,es
1506 !
1507 !---------------------------------------------------------------------
1509 thz(k)=theiz(k)-(xLvocp*qvz(k)-xLfocp*qiz(k))/tothz(k)
1511 tem(k)=tothz(k)*thz(k)
1512 if (tem(k) .gt. 273.15) then
1513 ! qsat=episp0k/prez(k)* &
1514 ! exp( svp2*(tem(k)-273.15)/(tem(k)-svp3) )
1515 es=1000.*svp1*exp( svp2*(tem(k)-svpt0)/(tem(k)-svp3) )
1516 qsat=ep2*es/(prez(k)-es)
1517 else
1518 qsat=episp0k/prez(k)* &
1519 exp( 21.8745584*(tem(k)-273.15)/(tem(k)-7.66) )
1520 end if
1521 qpz=qvz(k)+qlz(k)+qiz(k)
1522 if (qpz .lt. qsat) then
1523 qvz(k)=qpz
1524 qiz(k)=0.0
1525 qlz(k)=0.0
1526 go to 400
1527 end if
1528 qlpqi=qlz(k)+qiz(k)
1529 if( qlpqi .ge. 1.0e-5) then
1530 ratql=qlz(k)/qlpqi
1531 ratqi=qiz(k)/qlpqi
1532 else
1533 t0=273.15
1534 ! t1=233.15
1535 t1=248.15
1536 tmp1=( t0-tem(k) )/(t0-t1)
1537 tmp1=amin1(1.0,tmp1)
1538 tmp1=amax1(0.0,tmp1)
1539 ratqi=tmp1
1540 ratql=1.0-tmp1
1541 end if
1542 !
1543 !
1544 !-- saturation mixing ratios over water and ice
1545 !-- at the outset we will follow Bolton 1980 MWR for
1546 !-- the water and Murray JAS 1967 for the ice
1547 !
1548 !-- dqvsbar=d(qvsbar)/dT
1549 !-- ftsat=F(Tsat)
1550 !-- dftsat=d(F(T))/dT
1551 !
1552 ! First guess of tsat
1554 tsat=tem(k)
1555 absft=1.0
1556 !
1557 do 200 n=1,20
1558 denom1=1.0/(tsat-svp3)
1559 denom2=1.0/(tsat-7.66)
1560 ! qswz(k)=episp0k/prez(k)* &
1561 ! exp( svp2*denom1*(tsat-273.15) )
1562 es=1000.*svp1*exp( svp2*denom1*(tsat-svpt0) )
1563 qswz(k)=ep2*es/(prez(k)-es)
1564 if (tem(k) .lt. 273.15) then
1565 ! qsiz(k)=episp0k/prez(k)* &
1566 ! exp( 21.8745584*denom2*(tsat-273.15) )
1567 es=1000.*svp1*exp( 21.8745584*denom2*(tsat-273.15) )
1568 qsiz(k)=ep2*es/(prez(k)-es)
1569 if (tem(k) .lt. 233.15) qswz(k)=qsiz(k)
1570 else
1571 qsiz(k)=qswz(k)
1572 endif
1573 qvsbar(k)=ratql*qswz(k)+ratqi*qsiz(k)
1574 !
1575 ! if( absft .lt. 0.01 .and. n .gt. 3 ) go to 300
1576 if( absft .lt. 0.01 ) go to 300
1577 !
1578 dqvsbar=ratql*qswz(k)*svp2*243.5*denom1*denom1+ &
1579 ratqi*qsiz(k)*21.8745584*265.5*denom2*denom2
1580 ftsat=tsat+(xlvocp+ratqi*xlfocp)*qvsbar(k)- &
1581 tothz(k)*theiz(k)-xlfocp*ratqi*(qvz(k)+qlz(k)+qiz(k))
1582 dftsat=1.0+(xlvocp+ratqi*xlfocp)*dqvsbar
1583 tsat=tsat-ftsat/dftsat
1584 absft=abs(ftsat)
1586 200 continue
1587 9020 format(1x,'point can not converge, absft,n=',e12.5,i5)
1588 300 continue
1590 if( qpz .gt. qvsbar(k) ) then
1591 qvz(k)=qvsbar(k)
1592 qiz(k)=ratqi*( qpz-qvz(k) )
1593 qlz(k)=ratql*( qpz-qvz(k) )
1594 else
1595 qvz(k)=qpz
1596 qiz(k)=0.0
1597 qlz(k)=0.0
1598 end if
1599 400 continue
1601 END SUBROUTINE satadj
1604 !----------------------------------------------------------------
1605 REAL FUNCTION parama1(temp)
1606 !----------------------------------------------------------------
1607 IMPLICIT NONE
1608 !----------------------------------------------------------------
1609 ! This program calculate the parameter for crystal growth rate
1610 ! in Bergeron process
1611 !----------------------------------------------------------------
1613 REAL, INTENT (IN ) :: temp
1614 REAL, DIMENSION(32) :: a1
1615 INTEGER :: i1, i1p1
1616 REAL :: ratio
1618 data a1/0.100e-10,0.7939e-7,0.7841e-6,0.3369e-5,0.4336e-5, &
1619 0.5285e-5,0.3728e-5,0.1852e-5,0.2991e-6,0.4248e-6, &
1620 0.7434e-6,0.1812e-5,0.4394e-5,0.9145e-5,0.1725e-4, &
1621 0.3348e-4,0.1725e-4,0.9175e-5,0.4412e-5,0.2252e-5, &
1622 0.9115e-6,0.4876e-6,0.3473e-6,0.4758e-6,0.6306e-6, &
1623 0.8573e-6,0.7868e-6,0.7192e-6,0.6513e-6,0.5956e-6, &
1624 0.5333e-6,0.4834e-6/
1626 i1=int(-temp)+1
1627 i1p1=i1+1
1628 ratio=-(temp)-float(i1-1)
1629 parama1=a1(i1)+ratio*( a1(i1p1)-a1(i1) )
1631 END FUNCTION parama1
1633 !----------------------------------------------------------------
1634 REAL FUNCTION parama2(temp)
1635 !----------------------------------------------------------------
1636 IMPLICIT NONE
1637 !----------------------------------------------------------------
1638 ! This program calculate the parameter for crystal growth rate
1639 ! in Bergeron process
1640 !----------------------------------------------------------------
1642 REAL, INTENT (IN ) :: temp
1643 REAL, DIMENSION(32) :: a2
1644 INTEGER :: i1, i1p1
1645 REAL :: ratio
1647 data a2/0.0100,0.4006,0.4831,0.5320,0.5307,0.5319,0.5249, &
1648 0.4888,0.3849,0.4047,0.4318,0.4771,0.5183,0.5463, &
1649 0.5651,0.5813,0.5655,0.5478,0.5203,0.4906,0.4447, &
1650 0.4126,0.3960,0.4149,0.4320,0.4506,0.4483,0.4460, &
1651 0.4433,0.4413,0.4382,0.4361/
1652 i1=int(-temp)+1
1653 i1p1=i1+1
1654 ratio=-(temp)-float(i1-1)
1655 parama2=a2(i1)+ratio*( a2(i1p1)-a2(i1) )
1657 END FUNCTION parama2
1659 !+---+-----------------------------------------------------------------+
1660 ! THIS FUNCTION CALCULATES THE LIQUID SATURATION VAPOR MIXING RATIO AS
1661 ! A FUNCTION OF TEMPERATURE AND PRESSURE
1662 !
1663 REAL FUNCTION RSLF(P,T)
1665 IMPLICIT NONE
1666 REAL, INTENT(IN):: P, T
1667 REAL:: ESL,X
1668 REAL, PARAMETER:: C0= .611583699E03
1669 REAL, PARAMETER:: C1= .444606896E02
1670 REAL, PARAMETER:: C2= .143177157E01
1671 REAL, PARAMETER:: C3= .264224321E-1
1672 REAL, PARAMETER:: C4= .299291081E-3
1673 REAL, PARAMETER:: C5= .203154182E-5
1674 REAL, PARAMETER:: C6= .702620698E-8
1675 REAL, PARAMETER:: C7= .379534310E-11
1676 REAL, PARAMETER:: C8=-.321582393E-13
1678 X=MAX(-80.,T-273.16)
1680 ! ESL=612.2*EXP(17.67*X/(T-29.65))
1681 ESL=C0+X*(C1+X*(C2+X*(C3+X*(C4+X*(C5+X*(C6+X*(C7+X*C8)))))))
1682 RSLF=.622*ESL/(P-ESL)
1684 END FUNCTION RSLF
1685 !
1686 !+---+-----------------------------------------------------------------+
1687 ! THIS FUNCTION CALCULATES THE ICE SATURATION VAPOR MIXING RATIO AS A
1688 ! FUNCTION OF TEMPERATURE AND PRESSURE
1689 !
1690 REAL FUNCTION RSIF(P,T)
1692 IMPLICIT NONE
1693 REAL, INTENT(IN):: P, T
1694 REAL:: ESI,X
1695 REAL, PARAMETER:: C0= .609868993E03
1696 REAL, PARAMETER:: C1= .499320233E02
1697 REAL, PARAMETER:: C2= .184672631E01
1698 REAL, PARAMETER:: C3= .402737184E-1
1699 REAL, PARAMETER:: C4= .565392987E-3
1700 REAL, PARAMETER:: C5= .521693933E-5
1701 REAL, PARAMETER:: C6= .307839583E-7
1702 REAL, PARAMETER:: C7= .105785160E-9
1703 REAL, PARAMETER:: C8= .161444444E-12
1705 X=MAX(-80.,T-273.16)
1706 ESI=C0+X*(C1+X*(C2+X*(C3+X*(C4+X*(C5+X*(C6+X*(C7+X*C8)))))))
1707 RSIF=.622*ESI/(P-ESI)
1709 END FUNCTION RSIF
1710 !+---+-----------------------------------------------------------------+
1712 !----------------------------------------------------------------
1713 REAL FUNCTION ggamma(X)
1714 !----------------------------------------------------------------
1715 IMPLICIT NONE
1716 !----------------------------------------------------------------
1717 REAL, INTENT(IN ) :: x
1718 REAL, DIMENSION(8) :: B
1719 INTEGER ::j, K1
1720 REAL ::PF, G1TO2 ,TEMP
1722 DATA B/-.577191652,.988205891,-.897056937,.918206857, &
1723 -.756704078,.482199394,-.193527818,.035868343/
1725 PF=1.
1726 TEMP=X
1727 DO 10 J=1,200
1728 IF (TEMP .LE. 2) GO TO 20
1729 TEMP=TEMP-1.
1730 10 PF=PF*TEMP
1731 ! 100 FORMAT(//,5X,'module_mp_lin: INPUT TO GAMMA FUNCTION TOO LARGE, X=',E12.5)
1732 ! WRITE(wrf_err_message,100)X
1733 ! CALL wrf_error_fatal(wrf_err_message)
1734 20 G1TO2=1.
1735 TEMP=TEMP - 1.
1736 DO 30 K1=1,8
1737 30 G1TO2=G1TO2 + B(K1)*TEMP**K1
1738 ggamma=PF*G1TO2
1740 END FUNCTION ggamma
1742 !----------------------------------------------------------------
1744 END MODULE module_mp_sbu_ylin