| blob | 8816a75fb539ccb595324da573bd35116f806d1f |
1 !added
2 !WRF:MODEL_LAYER:PHYSICS
3 !based on module_mp_gsfcgce_311_new_20101115_newdbz.F
4 !implemented conv/stra separation scheme
5 !--->cal_sepa
6 !on 06/11/2011
7 !
8 !Implemented conv/stra diagnostic variable rainncv_sepa and rainnc_sepa
9 !--->var_sepa
10 !on 11/18/2011
11 !Di
12 !
13 !added small fix to prevent underflow and overflow issues
14 !on 12/02/2011
15 !EMK
16 !
17 !implemented 4ice writen by Steve Lang
18 !added tervrh.F_v3.0
19 !on 08/06/2012
20 !
21 !adding saticerh.F_v3.0
22 !on 08/07/2012
23 !
24 !from GCE saticerh.F_v3.08j, sgmap.F_v3.08t.clean
25 !consatrh.F_v3.08t, tervrh.F_v3.08t, and vqrqi.F_v3.05
26 !on 12/09/2014
27 !Di
29 MODULE module_mp_gsfcgce_4ice_nuwrf
31 #if (WRF_CHEM == 1)
32 use module_gocart_coupling
33 #endif
34 USE module_mp_radar
36 INTEGER, PARAMETER, PRIVATE:: chunk = 16
38 LOGICAL, EXTERNAL :: wrf_dm_on_monitor
40 !JJS 20140117 vvvvv
41 PRIVATE ! privatize all variables/subroutines in this module excepting public parameter below
42 PUBLIC :: gsfcgce_4ice_nuwrf
43 !JJS 20140117 ^^^^^
45 !JJS 1/3/2008 vvvvv
47 ! common /bt/
48 REAL, PRIVATE :: rd1, rd2, al, cp
50 ! common /cont/
51 REAL, PRIVATE :: c38, c358, c610, c149, &
52 c879, c172, c409, c76, &
53 c218, c580, c141
54 ! common /b3cs/
55 REAL, PRIVATE :: ag, bg, as, bs, &
56 aw, bw, bgh, bgq, &
57 bsh, bsq, bwh, bwq, &
58 ah, bh, bh3, bhh, &
59 bhh5, bhq, bh3_2
61 ! common /size/
62 REAL, PRIVATE :: tnw, tns, tng, tnh, &
63 roqs, roqg, roqr, roqh
65 ! common /rterv/
66 REAL, PRIVATE :: zrc, zgc, zsc, zhc, vrc, &
67 vrc0, vrc1, vrc2, vrc3, &
68 vgc, vsc, vhc
69 ! common /rterv_2/
70 REAL, PRIVATE :: zgc2, vgc2
72 ! common /bsnw/
73 REAL, PRIVATE :: rn11a
75 REAL, PRIVATE :: rn17, rn19b, &
76 bnd3, rn23a, &
77 rn23b, rn30b, &
78 rn30c
80 ! common /rsnw/
81 REAL, PRIVATE :: alv, alf, als, t0, t00, &
82 avc, afc, asc, esi, rn1, rn2, &
83 bnd2, rn3, rn4, rn5, rn50, &
84 rn51, rn52, rn53, rn6, rn60, &
85 rn61, rn62, rn63, rn7, rn8, &
86 rn9, rn10, rn101, rn102, rn10a, &
87 rn10b, rn10c, rn11, rn12, rn14, &
88 rn15, rn15a, rn16, rn171, rn172, &
89 rn17a, rn17b, rn17c, rn18, rn18a, &
90 rn19, rn191, rn192, rn19a, rn20, &
91 rn20a, rn20b, rn30, rn30a, rn21, &
92 bnd21, rn22, rn23, rn231, rn232, &
93 rn25, rn31, beta, rn32, rn33, &
94 rn331, rn332, rn34, rn35,rnn30a, &
95 rnn191,rnn192
96 REAL, PRIVATE, DIMENSION( 31 ) :: rn12a, rn12b, rn13, rn25a
98 ! common /rsnw2h/
99 REAL, PRIVATE :: hn9, hn10, hn10a, hn14, hn15a, hn16, &
100 hn17, hn17a, hn19, hn19a, hn20, hn20b
101 REAL, PRIVATE :: gn17,gn17a,gn17a2 !4ice revised
103 ! common /rainmap/
104 REAL, PRIVATE :: draimax
106 ! common /snomap/
107 real, PRIVATE :: xs,sno11,sno00,dsno11,dsno00,sexp11,sexp00,stt, &
108 stexp,sbase,tslopes,dsnomin,dsnomin4,slim
109 ! common /grpmap/
110 real, PRIVATE :: xg,grp11,grp00,dgrp11,dgrp00,gexp11,gexp00,gtt, &
111 gtexp,gbase,tslopeg,dgrpmin,dgrpmin4,glim
112 ! common /haimap/
113 real, PRIVATE :: hai00,hai11,htt0,htt1,haixp
115 ! common /b3cg_2/
116 REAL, PRIVATE :: ag2, bg2, bgh2, bgq2, roqg2, qrog2
118 ! common /rsnw2/
119 REAL, PRIVATE :: rn142, rn152, rn15a2, rn17a2, rn192_2
121 ! common /icemass/
122 REAL, PRIVATE :: ami50, ami40, ami100
124 ! common /BergCon/
125 REAL, PRIVATE, DIMENSION( 31 ) :: BergCon1, BergCon2, &
126 BergCon3, BergCon4
127 REAL, PRIVATE :: cmin
128 REAL, PRIVATE :: cpi
129 !
130 REAL, PRIVATE, DIMENSION( 31 ) :: aa1, aa2
131 DATA aa1/.7939e-7, .7841e-6, .3369e-5, .4336e-5, .5285e-5, &
132 .3728e-5, .1852e-5, .2991e-6, .4248e-6, .7434e-6, &
133 .1812e-5, .4394e-5, .9145e-5, .1725e-4, .3348e-4, &
134 .1725e-4, .9175e-5, .4412e-5, .2252e-5, .9115e-6, &
135 .4876e-6, .3473e-6, .4758e-6, .6306e-6, .8573e-6, &
136 .7868e-6, .7192e-6, .6513e-6, .5956e-6, .5333e-6, &
137 .4834e-6/
138 DATA aa2/.4006, .4831, .5320, .5307, .5319, &
139 .5249, .4888, .3894, .4047, .4318, &
140 .4771, .5183, .5463, .5651, .5813, &
141 .5655, .5478, .5203, .4906, .4447, &
142 .4126, .3960, .4149, .4320, .4506, &
143 .4483, .4460, .4433, .4413, .4382, &
144 .4361/
146 !+---+-----------------------------------------------------------------+
147 !..The following 6 variables moved here to facilitate reflectivity
148 !.. calculation similar to other MP schemes, because when they get
149 !.. declared later in the code (now commented out), it makes things
150 !.. more difficult to integreate with the radar code.
151 REAL :: xnor, xnos, xnoh, xnog
152 REAL :: rhohail, rhograul
153 !.. Values will be defined in subroutine fall_flux --- JJS 20150731
154 ! REAL , PARAMETER :: xnor = 8.0e6
155 ! REAL , PARAMETER :: xnos = 1.6e7
156 ! REAL , PARAMETER :: xnoh = 2.0e5
157 ! REAL , PARAMETER :: xnog = 4.0e6
158 ! REAL , PARAMETER :: rhohail = 917.
159 ! REAL , PARAMETER :: rhograul = 400.
160 !+---+-----------------------------------------------------------------+
162 !JJS 1/3/2008 ^^^^^
164 CONTAINS
166 !--------------------------------------------------------------------
167 ! NASA/GSFC GCE
168 ! Tao et al, 2001, Meteo. & Atmos. Phy., 97-137
169 !--------------------------------------------------------------------
170 SUBROUTINE gsfcgce_4ice_nuwrf( th &
171 ,qv, ql &
172 ,qr, qi &
173 ,qs, qh & ! 4ice
174 ,rho, pii, p, dt_in, z &
175 ,ht, dz8w, grav, w &
176 ,rhowater, rhosnow &
177 ,itimestep, xland, dx &
178 ,ids,ide, jds,jde, kds,kde & ! domain dims
179 ,ims,ime, jms,jme, kms,kme & ! memory dims
180 ,its,ite, jts,jte, kts,kte & ! tile dims
181 ,rainnc, rainncv &
182 ,snownc, snowncv, sr &
183 ,graupelnc, graupelncv &
184 ,refl_10cm, diagflag, do_radar_ref &
185 ,hailnc, hailncv & !Hail
186 ,f_qg, qg &
187 ,physc, physe, physd, physs, physm, physf &
188 ,acphysc, acphyse, acphysd, acphyss, acphysm, acphysf &
189 ,re_cloud_gsfc, re_rain_gsfc, re_ice_gsfc &
190 ,re_snow_gsfc, re_graupel_gsfc, re_hail_gsfc & ! cloud effective radius
191 ,preci3d, precs3d, precg3d, prech3d, precr3d &
192 #if ( WRF_CHEM == 1)
193 !JJS 20110525 vvvvv
194 ,aero, icn_diag, nc_diag, gid &
195 !JJS 20110525 ^^^^^
196 ! EMK
197 ,chem_opt &
198 ,gsfcgce_gocart_coupling &
199 #endif
200 !NUWRF END
201 )
204 !-------------------------------------------------------------------
205 IMPLICIT NONE
206 !-------------------------------------------------------------------
207 !
208 ! JJS 2/15/2005
209 !
210 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde , &
211 ims,ime, jms,jme, kms,kme , &
212 its,ite, jts,jte, kts,kte
213 INTEGER, INTENT(IN ) :: itimestep
215 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
216 INTENT(INOUT) :: &
217 th, &
218 qv, &
219 ql, &
220 qr, &
221 qi, &
222 qs, &
223 qg, &
224 qh
226 !NUWRF BEGIN
227 #if ( WRF_CHEM == 1)
228 ! JJS 20110525 vvvvv
229 ! for inline Gocart coupling
230 INTEGER, PARAMETER :: num_go = 14 ! number of the gocart aerosol species
231 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_go), intent(in) :: aero
232 REAL, DIMENSION( ims:ime, kms:kme, jms:jme), intent(out) :: icn_diag, nc_diag
233 INTEGER, INTENT(IN ) :: gid
234 ! JJS 20110525 ^^^^^
235 integer,intent(in) :: chem_opt ! EMK
236 integer,intent(in) :: gsfcgce_gocart_coupling ! EMK
237 #endif
238 !NUWRF END
240 !
241 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
242 INTENT(IN ) :: &
243 rho, &
244 pii, &
245 p, &
246 dz8w, &
247 z, &
248 w
250 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
251 INTENT(INOUT) :: &
252 physc, &
253 physe, &
254 physd, &
255 physs, &
256 physm, &
257 physf, &
258 acphysc, &
259 acphyse, &
260 acphysd, &
261 acphyss, &
262 acphysm, &
263 acphysf, &
264 preci3d, &
265 precs3d, &
266 precg3d, &
267 prech3d, &
268 precr3d
270 REAL, DIMENSION( ims:ime , jms:jme ), &
271 INTENT(INOUT) :: rainnc, &
272 rainncv, &
273 snownc, &
274 snowncv, &
275 sr, &
276 graupelnc, &
277 graupelncv,&
278 hailnc, &
279 hailncv
281 !JJS 20140225 for calculation of effective radius of cloud species
282 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: XLAND
283 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
284 INTENT(INOUT) :: re_cloud_gsfc, &
285 re_rain_gsfc, &
286 re_ice_gsfc, &
287 re_snow_gsfc, &
288 re_graupel_gsfc, &
289 re_hail_gsfc
290 !JJS 20140225 ^^^^^
292 !+---+-----------------------------------------------------------------+
293 REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: & ! GT
294 refl_10cm
295 LOGICAL, OPTIONAL, INTENT(IN) :: diagflag
296 INTEGER, OPTIONAL, INTENT(IN) :: do_radar_ref
297 !+---+-----------------------------------------------------------------+
299 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: ht
301 REAL, INTENT(IN ) :: dt_in, &
302 grav, &
303 rhowater, &
304 rhosnow, &
305 dx
307 LOGICAL, INTENT(IN), OPTIONAL :: F_QG
309 ! LOCAL VAR
311 !
312 INTEGER :: itaobraun, istatmin, new_ice_sat, id
313 INTEGER :: improve
315 INTEGER :: i, j, k, ip, ii, ic
316 INTEGER :: iskip, ih, icount, ibud, i24h
317 REAL :: hour
318 REAL :: dth, dqv, dqrest, dqall, dqall1, rhotot, a1, a2
320 REAL, DIMENSION(CHUNK, kms:kme):: th2d, qv2d, ql2d, qr2d
321 REAL, DIMENSION(CHUNK, kms:kme):: qi2d, qs2d, qg2d, qh2d
322 REAL, DIMENSION(CHUNK, kms:kme):: rho2d, pii2d, p2d, w2d
323 REAL, DIMENSION(CHUNK, kms:kme):: refc2d, refr2d, refi2d
324 REAL, DIMENSION(CHUNK, kms:kme):: refs2d, refg2d, refh2d
325 REAL, DIMENSION(CHUNK, kms:kme):: physc2d, physe2d, physd2d
326 REAL, DIMENSION(CHUNK, kms:kme):: physs2d, physm2d, physf2d
327 REAL, DIMENSION(CHUNK, kms:kme):: acphysc2d, acphyse2d, acphysd2d
328 REAL, DIMENSION(CHUNK, kms:kme):: acphyss2d, acphysm2d, acphysf2d
329 REAL, DIMENSION(CHUNK) :: xland1d
330 REAL, DIMENSION(CHUNK, kms:kme):: refl_10cm2d
331 #if ( WRF_CHEM == 1)
332 REAL, DIMENSION(CHUNK, kms:kme, num_go):: aero3d
333 REAL, DIMENSION(CHUNK, kms:kme):: icn_diag2d, nc_diag2d
334 #endif
336 !-----------------------------------------------------------------------
337 !WRF radar reflectivity initialization only need to run once
339 INTEGER:: NCALL=0
341 !-----------------------------------------------------------------------
344 ! itaobraun: 0 for Tao's constantis, 1 for Braun's constants
345 !c if ( itaobraun.eq.1 ) --> betah=0.5*beta=-.46*0.5=-0.23; cn0=1.e-6
346 !c if ( itaobraun.eq.0 ) --> betah=0.5*beta=-.6*0.5=-0.30; cn0=1.e-8
347 itaobraun = 0
349 ! Use Steve's new improvement 9/18/2009
351 improve = 8
353 !c new_ice_sat = 0, 1, 2, or 3
354 new_ice_sat = 9
356 !c istatmin
357 istatmin = 180
359 !c id = 0 without in-line staticstics
360 !c id = 1 with in-line staticstics
361 id = 0
363 !c ibud = 0 no calculation of dth, dqv, dqrest and dqall
364 !c ibud = 1 yes
365 ibud = 0
367 !c set up constants used internally in GCE
369 call consat_s ( itaobraun)
371 ! calculte fallflux and precipiation in MKS system
373 call fall_flux( dt_in,ql, qr, qi, qs, qg, qh, p, &
374 rho, th, pii, z, dz8w, ht, rainnc, &
375 rainncv, grav,itimestep, &
376 preci3d, precs3d, precg3d, prech3d, precr3d, &
377 snownc, snowncv, sr, &
378 graupelnc, graupelncv, &
379 hailnc, hailncv, &
380 vgc, vgc2,vhc,bhq, &
381 improve, &
382 ims,ime, jms,jme, kms,kme, & ! memory dims
383 its,ite, jts,jte, kts,kte ) ! tile dims
384 !-----------------------------------------------------------------------
385 ! EMK NUWRF...Moved this WRF radar reflectivity initialization to after
386 ! fall_flux, as the rhohail and rhograul variables are set in that
387 ! subroutine.
389 IF (NCALL .EQ. 0) THEN
390 !..Set these variables needed for computing radar reflectivity. These
391 !.. get used within radar_init to create other variables used in the
392 !.. radar module.
393 xam_r = 3.14159*rhowater/6.
394 xbm_r = 3.
395 xmu_r = 0.
396 xam_s = 3.14159*rhosnow/6.
397 xbm_s = 3.
398 xmu_s = 0.
399 xam_g = 3.14159*rhograul/6.
400 xbm_g = 3.
401 xmu_g = 0.
403 call radar_init
404 NCALL = 1
405 ENDIF
406 !+---+-----------------------------------------------------------------+
408 !c Negative values correction
410 ! iskip = 1
411 !
412 ! if (iskip.eq.0) then
413 ! call negcor(qv,rho,dz8w,ims,ime,jms,jme,kms,kme, &
414 ! itimestep,1, &
415 ! its,ite,jts,jte,kts,kte)
416 ! call negcor(ql,rho,dz8w,ims,ime,jms,jme,kms,kme, &
417 ! itimestep,2, &
418 ! its,ite,jts,jte,kts,kte)
419 ! call negcor(qr,rho,dz8w,ims,ime,jms,jme,kms,kme, &
420 ! itimestep,3, &
421 ! its,ite,jts,jte,kts,kte)
422 ! call negcor(qi,rho,dz8w,ims,ime,jms,jme,kms,kme, &
423 ! itimestep,4, &
424 ! its,ite,jts,jte,kts,kte)
425 ! call negcor(qs,rho,dz8w,ims,ime,jms,jme,kms,kme, &
426 ! itimestep,5, &
427 ! its,ite,jts,jte,kts,kte)
428 ! call negcor(qg,rho,dz8w,ims,ime,jms,jme,kms,kme, &
429 ! itimestep,6, &
430 ! its,ite,jts,jte,kts,kte)
431 !! else if (mod(itimestep,i24h).eq.1) then
432 !! print *,'no neg correction in mp at timestep=',itimestep
433 ! endif ! iskip
435 !c microphysics in GCE
437 !$OMP PARALLEL DO &
438 !$OMP PRIVATE ( ic, j, ii, i, k ) &
439 !$OMP PRIVATE ( th2d, qv2d, ql2d, qr2d ) &
440 !$OMP PRIVATE ( qi2d, qs2d, qg2d, qh2d ) &
441 !$OMP PRIVATE ( rho2d, pii2d, p2d, w2d ) &
442 !$OMP PRIVATE ( refc2d, refr2d, refi2d, refs2d, refg2d, refh2d ) &
443 !$OMP PRIVATE ( physc2d, physe2d, physd2d, physs2d, physm2d, physf2d ) &
444 !$OMP PRIVATE ( acphysc2d, acphyse2d, acphysd2d, acphyss2d, acphysm2d, acphysf2d ) &
445 !$OMP PRIVATE ( xland1d, refl_10cm2d ) &
446 #if ( WRF_CHEM == 1)
447 !$OMP PRIVATE ( aero3d, icn_diag2d, nc_diag2d) &
448 #endif
449 !$OMP SCHEDULE(dynamic,1)
450 DO ip = 1,((1+(ite-its+1)/CHUNK)*CHUNK)*(jte-jts+1),CHUNK ! i-dim contains '(1+(ite-its+1)/CHUNK)' blocks of size 'CHUNK'
451 j = jts+(ip-1)/((1+(ite-its+1)/CHUNK)*CHUNK)
452 IF ( j .ge. jts .and. j .le. jte ) THEN ! j: [jts, jte]
453 ii = its+mod((ip-1),((1+(ite-its+1)/CHUNK)*CHUNK)) ! ii: [its, ((1+(ite-its+1)/CHUNK)*CHUNK)]
455 DO ic=1,min(CHUNK,ite-ii+1)
456 i = ii+ic -1
457 xland1d(ic) = xland(i,j)
458 ENDDO
460 do k = kts, kte
461 DO ic=1,min(CHUNK,ite-ii+1)
462 i = ii+ic -1
463 th2d(ic,k) = th(i,k,j)
464 qv2d(ic,k) = qv(i,k,j)
465 ql2d(ic,k) = ql(i,k,j)
466 qr2d(ic,k) = qr(i,k,j)
467 qi2d(ic,k) = qi(i,k,j)
468 qs2d(ic,k) = qs(i,k,j)
469 qg2d(ic,k) = qg(i,k,j)
470 qh2d(ic,k) = qh(i,k,j)
471 rho2d(ic,k) = rho(i,k,j)
472 pii2d(ic,k) = pii(i,k,j)
473 p2d(ic,k) = p(i,k,j)
474 w2d(ic,k) = w(i,k,j)
475 refl_10cm2d(ic,k)=refl_10cm(i,k,j)
476 refc2d(ic,k)=re_cloud_gsfc(i,k,j)
477 refr2d(ic,k)=re_rain_gsfc(i,k,j)
478 refi2d(ic,k)=re_ice_gsfc(i,k,j)
479 refs2d(ic,k)=re_snow_gsfc(i,k,j)
480 refg2d(ic,k)=re_graupel_gsfc(i,k,j)
481 refh2d(ic,k)=re_hail_gsfc(i,k,j)
482 physc2d(ic,k)=physc(i,k,j)
483 physe2d(ic,k)=physe(i,k,j)
484 physd2d(ic,k)=physd(i,k,j)
485 physs2d(ic,k)=physs(i,k,j)
486 physm2d(ic,k)=physm(i,k,j)
487 physf2d(ic,k)=physf(i,k,j)
488 acphysc2d(ic,k)=acphysc(i,k,j)
489 acphyse2d(ic,k)=acphyse(i,k,j)
490 acphysd2d(ic,k)=acphysd(i,k,j)
491 acphyss2d(ic,k)=acphyss(i,k,j)
492 acphysm2d(ic,k)=acphysm(i,k,j)
493 acphysf2d(ic,k)=acphysf(i,k,j)
494 #if ( WRF_CHEM == 1)
495 aero3d(ic,k,:)=aero(i,k,j,:)
496 #endif
497 ENDDO
498 enddo
500 IF ( min(CHUNK,ite-ii+1) .gt. 0 ) THEN
501 call saticel_s( dt_in, dx, itaobraun, istatmin, &
502 new_ice_sat, id, improve, &
503 th2d, qv2d, ql2d, qr2d, &
504 qi2d, qs2d, qg2d, qh2d, &
505 rho2d, pii2d, p2d, w2d, &
506 itimestep, xland1d, &
507 refl_10cm2d, diagflag, do_radar_ref, & ! GT added for reflectivity calcs
508 ids,ide, jds,jde, kds,kde, & ! domain dims
509 ims,ime, jms,jme, kms,kme, & ! memory dims
510 its,ite, jts,jte, kts,kte, & ! tile dims
511 !NUWRF BEGIN
512 refc2d, refr2d, refi2d, refs2d, refg2d, refh2d, & ! cloud effective radius
513 physc2d, physe2d, physd2d, physs2d, physm2d, physf2d, &
514 acphysc2d, acphyse2d, acphysd2d, acphyss2d, acphysm2d, acphysf2d &
516 #if ( WRF_CHEM == 1)
517 !JJS 20110525 vvvvv
518 ,aero3d, icn_diag2d, nc_diag2d, gid, &
519 !JJS 20110525 ^^^^^
520 !EMK
521 chem_opt, &
522 gsfcgce_gocart_coupling &
523 #endif
524 ,ii, j, min(CHUNK,ite-ii+1))
525 !NUWRF END
526 ENDIF
528 do k = kts, kte
529 DO ic=1,min(CHUNK,ite-ii+1)
530 i = ii+ic -1
531 th(i,k,j) = th2d(ic,k)
532 qv(i,k,j) = qv2d(ic,k)
533 ql(i,k,j) = ql2d(ic,k)
534 qr(i,k,j) = qr2d(ic,k)
535 qi(i,k,j) = qi2d(ic,k)
536 qs(i,k,j) = qs2d(ic,k)
537 qg(i,k,j) = qg2d(ic,k)
538 qh(i,k,j) = qh2d(ic,k)
539 re_cloud_gsfc(i,k,j)=refc2d(ic,k)
540 re_rain_gsfc(i,k,j)=refr2d(ic,k)
541 re_ice_gsfc(i,k,j)=refi2d(ic,k)
542 re_snow_gsfc(i,k,j)=refs2d(ic,k)
543 re_graupel_gsfc(i,k,j)=refg2d(ic,k)
544 re_hail_gsfc(i,k,j)=refh2d(ic,k)
545 refl_10cm(i,k,j)=refl_10cm2d(ic,k)
546 physc(i,k,j)=physc2d(ic,k)
547 physe(i,k,j)=physe2d(ic,k)
548 physd(i,k,j)=physd2d(ic,k)
549 physs(i,k,j)=physs2d(ic,k)
550 physm(i,k,j)=physm2d(ic,k)
551 physf(i,k,j)=physf2d(ic,k)
552 acphysc(i,k,j)=acphysc2d(ic,k)
553 acphyse(i,k,j)=acphyse2d(ic,k)
554 acphysd(i,k,j)=acphysd2d(ic,k)
555 acphyss(i,k,j)=acphyss2d(ic,k)
556 acphysm(i,k,j)=acphysm2d(ic,k)
557 acphysf(i,k,j)=acphysf2d(ic,k)
558 #if ( WRF_CHEM == 1)
559 icn_diag(i,k,j)=icn_diag2d(ic,k)
560 nc_diag(i,k,j)=nc_diag2d(ic,k)
561 #endif
562 ENDDO
563 enddo
565 ENDIF
566 ENDDO ! ip_loop
568 END SUBROUTINE gsfcgce_4ice_nuwrf
570 SUBROUTINE fall_flux ( dt, ql, qr, qi, qs, qg, qh, p, &
571 rho, th, pi_mks, z, dz8w, topo, rainnc, &
572 rainncv, grav, itimestep, &
573 preci3d, precs3d, precg3d, prech3d, precr3d, &
574 snownc, snowncv, sr, &
575 graupelnc, graupelncv, &
576 hailnc, hailncv, &
577 vgc, vgc2, vhc, bhq, &
578 improve, &
579 ims,ime, jms,jme, kms,kme, & ! memory dims
580 its,ite, jts,jte, kts,kte ) ! tile dims
581 !-----------------------------------------------------------------------
582 ! adopted from Jiun-Dar Chern's codes for Purdue Regional Model
583 ! adopted by Jainn J. Shi, 6/10/2005
584 ! modified by Goddard 7/24/2010
585 ! modified by Tao 11/12/2010
586 !-----------------------------------------------------------------------
588 IMPLICIT NONE
589 INTEGER, INTENT(IN ) :: improve, &
590 ims,ime, jms,jme, kms,kme, &
591 its,ite, jts,jte, kts,kte
592 INTEGER, INTENT(IN ) :: itimestep
593 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
594 INTENT(INOUT) :: ql, qr, qi, qs, qg, qh
595 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
596 INTENT(IN) :: th, pi_mks
598 REAL, DIMENSION( ims:ime , jms:jme ), &
599 INTENT(INOUT) :: rainnc, rainncv, &
600 snownc, snowncv, sr, &
601 graupelnc, graupelncv, &
602 hailnc, hailncv
603 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
604 INTENT(IN ) :: rho, z, dz8w, p
606 REAL, INTENT(IN ) :: dt, grav, vgc, vgc2, vhc, bhq
609 REAL, DIMENSION( ims:ime , jms:jme ), &
610 INTENT(IN ) :: topo
611 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
612 INTENT(OUT) :: preci3d, precs3d, precg3d, prech3d, precr3d
614 ! temperary vars
616 REAL, DIMENSION( kts:kte ) :: fv
617 REAL :: tmp1, term0
618 REAL :: pptrain, pptsnow, &
619 pptgraul, pptice, ppthail
620 REAL :: qrz, qiz
621 REAL, DIMENSION( kts:kte ) :: qcz, qsz, qgz, qhz, &
622 zz, dzw, prez, rhoz, &
623 orhoz,r00
624 REAL, DIMENSION( kts:kte ) :: csed, rsed, ised, ssed, gsed, hsed
626 REAL, DIMENSION( kts:kte ) :: thz, piz
628 INTEGER :: k, i, j
629 !
631 REAL, DIMENSION( kts:kte ) :: vtr, vts, vtg, vth, vti
633 REAL :: dtb, pi, consta, constc, gambp4, &
634 gamdp4, gam4pt5, gam4bbar
636 ! New local variable
637 REAL :: y1, y2, vr, vs, vg
638 REAL :: vgcr, vgcr2, vscf, &
639 vhcr, vhcr2
640 REAL :: tair, tairc, fexp
641 REAL :: ftns, ftnsQ, ftng, ftngQ
642 REAL :: const_vt, const_d, const_m, bb1, bb2
643 REAL, DIMENSION(7) :: aice, vice
644 REAL :: ftns0, ftng0
645 REAL :: fros, fros0
646 REAL :: qhz2,qgz2
648 ! DATA tslopes/0./, tslopeg/0./
649 ! DATA const_vt//, const_d//, const_m//
650 DATA aice/1.e-6, 1.e-5, 1.e-4, 1.e-3, 0.01, 0.1, 1./ ! g/m**3
651 DATA vice/5.,15.,30.,35.,40.,45.,50./ ! cm/s
652 DATA ftns/1./, ftng/1./
655 ! will be defined later using consat values
656 REAL :: rhowater
657 REAL :: rhosnow
658 !JJS 20140116 These variables are declared in the beginning of the module.
659 !JJS 20140116 No need to declared again here.
660 ! REAL :: xnor
661 ! REAL :: xnos
662 ! REAL :: xnoh, rhohail
663 ! REAL :: xnog, rhograul
664 !
666 REAL , PARAMETER :: &
667 ! constb = 0.8, constd = 0.25, o6 = 1./6., &
668 constb = 0.8, constd = 0.11, o6 = 1./6., &
669 cdrag = 0.6
670 REAL , PARAMETER :: abar = 19.3, bbar = 0.37, &
671 p0 = 1.0e5
672 REAL , PARAMETER :: rhoe_s = 1.29
674 ! for terminal velocity flux
675 INTEGER :: min_q, max_q
676 REAL :: t_del_tv, del_tv, flux, fluxin, fluxout
677 LOGICAL :: notlast
679 !-----------------------------------------------------------------------
680 ! This program calculates precipitation fluxes due to terminal velocities.
681 !-----------------------------------------------------------------------
683 dtb=dt
684 pi=acos(-1.)
686 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
687 xnor = tnw*1.0e8
688 rhowater = roqr*1000.
689 xnos = tns*1.0e8 ! Consistent with ConSat
690 rhosnow = roqs*1000.
692 ! tng and roqg are assigned with hail numbers in consat if (ihail.eq.1)
693 xnog = tng*1.0e8
694 rhograul = roqg*1000.
695 xnoh = tnh*1.0e8
696 rhohail = roqh*1000.
698 consta=2115.0*0.01**(1-constb)
699 ! constc=152.93*0.01**(1-constd)
700 constc=78.63*0.01**(1-constd)
702 ! Gamma function
703 gambp4=gammagce(constb+4.)
704 gamdp4=gammagce(constd+4.)
705 gam4pt5=gammagce(4.5)
706 gam4bbar=gammagce(4.+bbar)
707 !
708 ! cmin=1.e-10
709 !
710 !***********************************************************************
711 ! Calculate precipitation fluxes due to terminal velocities.
712 !***********************************************************************
713 !
714 !- Calculate termianl velocity (vt?) of precipitation q?z
715 !- Find maximum vt? to determine the small delta t
717 !$OMP PARALLEL DO &
718 !$OMP FIRSTPRIVATE(dtb, pi, rhowater, rhosnow, gambp4) &
719 !$OMP FIRSTPRIVATE(consta, constc, gamdp4, gam4pt5, gam4bbar) &
720 !$OMP PRIVATE(i,k,fv,tmp1,term0,pptrain, pptsnow,pptgraul, pptice, ppthail) &
721 !$OMP PRIVATE(qcz, qrz, qiz, qsz, qgz, qhz, zz, dzw, prez, rhoz, orhoz,r00) &
722 !$OMP PRIVATE(csed, rsed, ised, ssed, gsed, hsed) &
723 !$OMP PRIVATE(thz, piz) &
724 !$OMP PRIVATE(vtr, vts, vtg, vth, vti) &
725 !$OMP PRIVATE(y1, y2, vr, vs, vg) &
726 !$OMP PRIVATE(vgcr, vgcr2, vscf, vhcr, vhcr2) &
727 !$OMP PRIVATE(tair, tairc, fexp) &
728 !$OMP PRIVATE(ftns, ftnsQ, ftng, ftngQ) &
729 !$OMP PRIVATE(const_vt, const_d, const_m, bb1, bb2) &
730 !$OMP PRIVATE(aice, vice, ftns0, ftng0, fros, fros0, qhz2,qgz2) &
731 !$OMP PRIVATE(min_q, max_q) &
732 !$OMP PRIVATE(t_del_tv, del_tv, flux, fluxin, fluxout) &
733 !$OMP PRIVATE(notlast) &
734 !$OMP SCHEDULE(dynamic)
736 j_loop: do j = jts, jte
737 i_loop: do i = its, ite
739 do k = kts, kte
740 preci3d(i,k,j)=0.
741 precs3d(i,k,j)=0.
742 precg3d(i,k,j)=0.
743 prech3d(i,k,j)=0.
744 precr3d(i,k,j)=0.
745 ised(k)=0.
746 ssed(k)=0.
747 gsed(k)=0.
748 hsed(k)=0.
749 rsed(k)=0.
750 end do
752 pptrain = 0.
753 pptsnow = 0.
754 pptgraul = 0.
755 ppthail = 0.
756 pptice = 0.
758 ! in MKS system
759 do k = kts, kte
760 qcz(k)=ql(i,k,j) !Di
761 qsz(k)=qs(i,k,j)
762 qhz(k)=qh(i,k,j)
763 rhoz(k)=rho(i,k,j)
764 r00(k)=rhoz(k)*0.001 !rho in cgs
765 thz(k)=th(i,k,j)
766 piz(k)=pi_mks(i,k,j)
767 orhoz(k)=1./rhoz(k)
768 prez(k)=p(i,k,j)
769 fv(k)=sqrt(rhoe_s/rhoz(k))
770 ! fv(k)=sqrt(rho(i,1,j)/rhoz(k))
771 zz(k)=z(i,k,j)
772 dzw(k)=dz8w(i,k,j)
773 enddo !k
775 DO k = kts, kte
776 qgz(k)=qg(i,k,j)
777 ENDDO
779 !
780 !-- rain
781 !
782 t_del_tv=0.
783 del_tv=dtb
784 notlast=.true.
785 DO while (notlast)
786 !
787 min_q=kte
788 max_q=kts-1
789 !
792 do k=kts,kte-1
794 vtr(k)=0.
795 qrz=qr(i,k,j)
796 if (qrz .gt. cmin) then
797 min_q=min0(min_q,k)
798 max_q=max0(max_q,k)
800 ! old codes from Chern's in MKS
801 ! min_q=min0(min_q,k)
802 ! max_q=max0(max_q,k)
803 tmp1=sqrt(pi*rhowater*xnor/rhoz(k)/qrz)
804 tmp1=sqrt(tmp1)
805 vtr(k)=consta*gambp4*fv(k)/tmp1**constb
806 vtr(k)=vtr(k)/6.
807 ! new codes from Steve's in cgs
808 tair=thz(k)*piz(k)
809 tairc=tair-t0
810 y1=qrz ! rhoz(k) need to be in CGS
811 y2=qcz(k) ! rhoz(k) need to be in CGS
812 call vqrqi(1,r00(k),fv(k),y1,y2,tair,vtr(k)) !Di
813 vtr(k)=vtr(k) * 0.01 ! convert back to MKS !Di
815 if (.not. vtr(k) .gt. 0.0) cycle ! EMK NUWRF Bug fix
817 if (k .eq. 1) then
818 del_tv=amin1(del_tv,0.9*(zz(k)-topo(i,j))/vtr(k))
819 else
820 del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtr(k))
821 endif
822 endif
823 enddo !do K
825 if (max_q .ge. min_q) then
826 !
827 !- Check if the summation of the small delta t >= big delta t
828 ! (t_del_tv) (del_tv) (dtb)
830 t_del_tv=t_del_tv+del_tv
831 !
832 if ( t_del_tv .ge. dtb ) then
833 notlast=.false.
834 del_tv=dtb+del_tv-t_del_tv
835 endif
837 ! use small delta t to calculate the qrz flux
838 ! termi is the qrz flux pass in the grid box through the upper boundary
839 ! termo is the qrz flux pass out the grid box through the lower boundary
840 !
841 fluxin=0.
842 do k=max_q,min_q,-1
843 qrz=qr(i,k,j)
844 fluxout=rhoz(k)*vtr(k)*qrz
845 flux=(fluxin-fluxout)/rhoz(k)/dzw(k)
846 qrz=qrz+del_tv*flux
847 qrz=amax1(0.,qrz)
848 qr(i,k,j)=qrz
849 fluxin=fluxout
850 rsed(k)=rsed(k)+fluxin
851 enddo
852 if (min_q .eq. 1) then
853 pptrain=pptrain+fluxin*del_tv
854 else
855 qrz=qr(i,min_q-1,j)
856 qrz=qrz+del_tv* &
857 fluxin/rhoz(min_q-1)/dzw(min_q-1)
858 qrz=amax1(0.,qrz) !Di 10/23/2012
859 qr(i,min_q-1,j)=qrz
860 endif
861 !
862 else
863 notlast=.false.
864 endif
865 ENDDO ! DO WHILE
867 !
868 !-- snow
869 !
870 t_del_tv=0.
871 del_tv=dtb
872 notlast=.true.
874 DO while (notlast)
875 !
876 min_q=kte
877 max_q=kts-1
879 !
880 do k=kts,kte-1
881 vts(k)=0.
883 if (qsz(k) .gt. cmin) then
884 min_q=min0(min_q,k)
885 max_q=max0(max_q,k)
887 ! old codes from Chern's in MKS
888 tmp1=sqrt(pi*rhosnow*xnos/rhoz(k)/qsz(k))
889 tmp1=sqrt(tmp1)
890 vts(k)=constc*gamdp4*fv(k)/tmp1**constd
891 vts(k)=vts(k)/6.
893 ! new codes from Steve's in cgs
894 y1 = qsz(k)
895 vscf=vsc*fv(k)
897 ftns=1.
898 ftns0=1.
899 tair=thz(k)*piz(k)
900 tairc=tair-t0
902 qhz2=qhz(k)
903 qgz2=qgz(k)
904 if (k .lt. kte-2 .and. tairc .ge. -5) then
905 qhz2=qhz(k+1)
906 qgz2=qgz(k+1)
907 endif
908 call sgmap(1,qsz(k),qgz(k),qgz2,qhz(k),qhz2,r00(k),tairc,ftns0) !snow intercept
909 ftns=ftns0**bsq
910 fros=1.
911 call sgmap(3,qsz(k),qgz(k),qgz2,qhz(k),qhz2,r00(k),tairc,fros0) !snow density
912 fros=fros0**bsq
913 vts(k)=max(vscf*(r00(k)*y1)**bsq/ftns/fros,0.e0)
914 ! vts(k)=max(vscf*(r00*y1)**bsq/ftns, 0.0)
915 ! ! bs, bsq, vscf are defined in new consat_s
916 vts(k)=vts(k) * 0.01 ! convert back to MKS
918 if (k .eq. 1) then
919 del_tv=amin1(del_tv,0.9*(zz(k)-topo(i,j))/vts(k))
920 else
921 del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vts(k))
922 endif !k
923 endif !cmin
924 enddo ! do k
926 if (max_q .ge. min_q) then
927 !
928 !
929 !- Check if the summation of the small delta t >= big delta t
930 ! (t_del_tv) (del_tv) (dtb)
932 t_del_tv=t_del_tv+del_tv
934 if ( t_del_tv .ge. dtb ) then
935 notlast=.false.
936 del_tv=dtb+del_tv-t_del_tv
937 endif
939 ! use small delta t to calculate the qsz flux
940 ! termi is the qsz flux pass in the grid box through the upper boundary
941 ! termo is the qsz flux pass out the grid box through the lower boundary
942 !
943 fluxin=0.
944 do k=max_q,min_q,-1
945 fluxout=rhoz(k)*vts(k)*qsz(k)
946 flux=(fluxin-fluxout)/rhoz(k)/dzw(k)
947 qsz(k)=qsz(k)+del_tv*flux
948 qsz(k)=amax1(0.,qsz(k))
949 qs(i,k,j)=qsz(k)
950 fluxin=fluxout
951 ssed(k)=ssed(k)+fluxin
952 enddo
953 if (min_q .eq. 1) then
954 pptsnow=pptsnow+fluxin*del_tv
955 else
956 qsz(min_q-1)=qsz(min_q-1)+del_tv* &
957 fluxin/rhoz(min_q-1)/dzw(min_q-1)
958 qsz(min_q-1)=amax1(0.,qsz(min_q-1)) !Di 10/23/2012
959 qs(i,min_q-1,j)=qsz(min_q-1)
960 endif
961 !
962 else
963 notlast=.false.
964 endif
966 ENDDO
968 !
969 !--- graupel
970 !
972 t_del_tv=0.
973 del_tv=dtb
974 notlast=.true.
975 !
976 DO while (notlast)
977 !
978 min_q=kte
979 max_q=kts-1
980 !
981 do k=kts,kte-1
982 vtg(k)=0.
984 if (qgz(k) .gt. cmin) then
985 min_q=min0(min_q,k)
986 max_q=max0(max_q,k)
988 ! for graupel, based on RH (1984)
989 ! new codes from Steve's in cgs
990 ! y1=rhoz(k) * 0.001 * qgz(k) ! rhoz(k) need to be in CGS
991 y1 = qgz(k)
992 vgcr=vgc*fv(k)
993 vgcr2=vgc2*fv(k)
994 ftng=1.
995 ftng0=1.
996 tair=thz(k)*piz(k)
997 tairc=tair-t0
998 qhz2=qhz(k)
999 qgz2=qgz(k)
1000 if (k .lt. kte-2 .and. tairc .ge. -5) then
1001 qhz2=qhz(k+1)
1002 qgz2=qgz(k+1)
1003 endif
1004 call sgmap(2,qsz(k),qgz(k),qgz2,qhz(k),qhz2,r00(k),tairc,ftng0) !ftng0 is graupel intercept
1005 ftng=ftng0**bgq
1006 vtg(k)=amax1(vgcr*(r00(k)*y1)**bgq/ftng, 0.0)
1007 ! bg, vgcr, bgq are defined in new consat_s
1008 vtg(k)=vtg(k) * 0.01 ! convert back to MKS
1009 if (y1.gt.qrog2)then !Di
1010 ftng=ftng0**bgq2 !Di
1011 vtg(k)=amax1(vgcr2*(r00(k)*y1)**bgq2/ftng, 0.0) !Di
1012 ! bg, vgcr, bgq are defined in new consat_s !Di
1013 vtg(k)=vtg(k) * 0.01 ! convert back to MKS
1014 endif
1016 if (k .eq. 1) then
1017 del_tv=amin1(del_tv,0.9*(zz(k)-topo(i,j))/vtg(k))
1018 else
1019 del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vtg(k))
1020 endif
1021 !
1022 endif !qgz
1023 enddo !k
1025 if (max_q .ge. min_q) then
1026 !
1027 !
1028 !- Check if the summation of the small delta t >= big delta t
1029 ! (t_del_tv) (del_tv) (dtb)
1031 t_del_tv=t_del_tv+del_tv
1033 if ( t_del_tv .ge. dtb ) then
1034 notlast=.false.
1035 del_tv=dtb+del_tv-t_del_tv
1036 endif
1038 ! use small delta t to calculate the qgz flux
1039 ! termi is the qgz flux pass in the grid box through the upper boundary
1040 ! termo is the qgz flux pass out the grid box through the lower boundary
1041 !
1042 fluxin=0.
1043 do k=max_q,min_q,-1
1044 fluxout=rhoz(k)*vtg(k)*qgz(k)
1045 flux=(fluxin-fluxout)/rhoz(k)/dzw(k)
1046 qgz(k)=qgz(k)+del_tv*flux
1047 qgz(k)=amax1(0.,qgz(k))
1048 qg(i,k,j)=qgz(k)
1049 fluxin=fluxout
1050 gsed(k)=gsed(k)+fluxin
1051 enddo
1052 if (min_q .eq. 1) then
1053 pptgraul=pptgraul+fluxin*del_tv
1054 else
1055 qgz(min_q-1)=qgz(min_q-1)+del_tv* &
1056 fluxin/rhoz(min_q-1)/dzw(min_q-1)
1057 qgz(min_q-1)=amax1(0.,qgz(min_q-1)) !Di 10/23/2012
1058 qg(i,min_q-1,j)=qgz(min_q-1)
1059 endif
1060 !
1061 else
1062 notlast=.false.
1063 endif
1064 !
1065 ENDDO
1066 !
1067 !--- hail
1068 !
1070 t_del_tv=0.
1071 del_tv=dtb
1072 notlast=.true.
1073 !
1074 DO while (notlast)
1075 !
1076 min_q=kte
1077 max_q=kts-1
1079 do k=kts,kte-1
1080 vth(k)=0.
1082 if (qhz(k) .gt. cmin) then
1083 min_q=min0(min_q,k)
1084 max_q=max0(max_q,k)
1086 ! new codes from Steve's in cgs
1087 y1 = qhz(k)
1088 vhcr=vhc/sqrt(r00(k)) !Di
1089 vth(k)=amax1(vhcr*(y1*r00(k))**bhq, 0.e0) !Di
1090 vth(k)=vth(k) * 0.01 ! convert back to MKS
1092 if (k .eq. 1) then
1093 del_tv=amin1(del_tv,0.9*(zz(k)-topo(i,j))/vth(k))
1094 else
1095 del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vth(k))
1096 endif
1097 !
1098 endif !qhz
1099 enddo !k
1101 if (max_q .ge. min_q) then
1102 !
1103 !
1104 !- Check if the summation of the small delta t >= big delta t
1105 ! (t_del_tv) (del_tv) (dtb)
1107 t_del_tv=t_del_tv+del_tv
1109 if ( t_del_tv .ge. dtb ) then
1110 notlast=.false.
1111 del_tv=dtb+del_tv-t_del_tv
1112 endif
1114 ! use small delta t to calculate the qhz flux
1115 ! termi is the qhz flux pass in the grid box through the upper boundary
1116 ! termo is the qhz flux pass out the grid box through the lower boundary
1117 !
1118 fluxin=0.
1119 do k=max_q,min_q,-1
1120 fluxout=rhoz(k)*vth(k)*qhz(k)
1121 flux=(fluxin-fluxout)/rhoz(k)/dzw(k)
1122 qhz(k)=qhz(k)+del_tv*flux
1123 qhz(k)=amax1(0.,qhz(k))
1124 qh(i,k,j)=qhz(k)
1125 fluxin=fluxout
1126 hsed(k)=hsed(k)+fluxin
1127 enddo
1128 if (min_q .eq. 1) then
1129 ppthail=ppthail+fluxin*del_tv
1130 else
1131 qhz(min_q-1)=qhz(min_q-1)+del_tv* &
1132 fluxin/rhoz(min_q-1)/dzw(min_q-1)
1133 qhz(min_q-1)=amax1(0.,qhz(min_q-1)) !Di 10/23/2012
1134 qh(i,min_q-1,j)=qhz(min_q-1)
1135 endif
1136 !
1137 else
1138 notlast=.false.
1139 endif
1140 !
1141 ENDDO
1143 !
1144 !-- cloud ice (03/21/02) follow Vaughan T.J. Phillips at GFDL
1145 !
1147 t_del_tv=0.
1148 del_tv=dtb
1149 notlast=.true.
1150 !
1151 DO while (notlast)
1152 !
1153 min_q=kte
1154 max_q=kts-1
1155 !
1156 do k=kts,kte-1
1157 qiz=qi(i,k,j)
1158 vti(k)=0.
1159 if (qiz .gt. cmin) then
1160 min_q=min0(min_q,k)
1161 max_q=max0(max_q,k)
1163 ! new codes from Steve's in cgs.
1164 vti(k)=0.
1165 y1=rhoz(k) * 1000. * qiz ! y1 in g/m**3
1166 if (y1 .ge. 1.e-6) then
1167 y1=qiz
1168 y2=qcz(k)
1169 tair=thz(k)*piz(k)
1170 tairc=tair-t0
1171 call vqrqi(2,r00(k),fv(k),y1,y2,tair,vti(k)) ! Di
1172 endif !y1
1173 vti(k)=vti(k) * 0.01 ! convert back to MKS
1175 ! EMK: prevent divsion by zero and underflow value
1176 if (vti(k) .gt. 1.e-20) then
1177 if (k .eq. 1) then
1178 del_tv=amin1(del_tv,0.9*(zz(k)-topo(i,j))/vti(k))
1179 else
1180 del_tv=amin1(del_tv,0.9*(zz(k)-zz(k-1))/vti(k))
1181 endif
1182 endif
1183 ! else
1184 ! vti(k)=0.
1185 endif
1186 enddo
1188 if (max_q .ge. min_q) then
1189 !
1190 !
1191 !- Check if the summation of the small delta t >= big delta t
1192 ! (t_del_tv) (del_tv) (dtb)
1194 t_del_tv=t_del_tv+del_tv
1196 if ( t_del_tv .ge. dtb ) then
1197 notlast=.false.
1198 del_tv=dtb+del_tv-t_del_tv
1199 endif
1201 ! use small delta t to calculate the qiz flux
1202 ! termi is the qiz flux pass in the grid box through the upper boundary
1203 ! termo is the qiz flux pass out the grid box through the lower boundary
1204 !
1206 fluxin=0.
1207 do k=max_q,min_q,-1
1208 qiz=qi(i,k,j)
1209 fluxout=rhoz(k)*vti(k)*qiz
1210 flux=(fluxin-fluxout)/rhoz(k)/dzw(k)
1211 qiz=qiz+del_tv*flux
1212 qiz=amax1(0.,qiz)
1213 qi(i,k,j)=qiz
1214 fluxin=fluxout
1215 ised(k)=ised(k)+fluxin
1216 enddo
1217 if (min_q .eq. 1) then
1218 pptice=pptice+fluxin*del_tv
1219 else
1220 qiz=qi(i,min_q-1,j)
1221 qiz=qiz+del_tv* &
1222 fluxin/rhoz(min_q-1)/dzw(min_q-1)
1223 qiz=amax1(0.,qiz) !Di 10/23/2012
1224 qi(i,min_q-1,j)=qiz
1225 endif
1226 !
1227 else
1228 notlast=.false.
1229 endif
1230 !
1231 ENDDO !notlast
1233 do k = kts, kte
1234 preci3d(i,k,j)=ised(k)
1235 precs3d(i,k,j)=ssed(k)
1236 precg3d(i,k,j)=gsed(k)
1237 prech3d(i,k,j)=hsed(k)
1238 precr3d(i,k,j)=rsed(k)
1239 end do
1241 ! prnc(i,j)=prnc(i,j)+pptrain
1242 ! psnowc(i,j)=psnowc(i,j)+pptsnow
1243 ! pgrauc(i,j)=pgrauc(i,j)+pptgraul
1244 ! picec(i,j)=picec(i,j)+pptice
1245 !
1247 ! write(6,*) 'i=',i,' j=',j,' ', pptrain, pptsnow, pptgraul, pptice
1248 ! call flush(6)
1250 snowncv(i,j) = pptsnow
1251 snownc(i,j) = snownc(i,j) + pptsnow
1252 graupelncv(i,j) = pptgraul
1253 graupelnc(i,j) = graupelnc(i,j) + pptgraul
1254 hailncv(i,j) = ppthail
1255 hailnc(i,j) = hailnc(i,j) + ppthail
1256 RAINNCV(i,j) = pptrain + pptsnow + pptgraul + pptice + ppthail
1257 RAINNC(i,j) = RAINNC(i,j) + pptrain + pptsnow + pptgraul + pptice + ppthail
1258 sr(i,j) = 0.
1259 if (RAINNCV(i,j) .gt. 0.) sr(i,j) = (pptsnow + pptgraul + pptice + ppthail) / RAINNCV(i,j)
1261 ENDDO i_loop
1262 ENDDO j_loop
1265 END SUBROUTINE fall_flux
1267 !-----------------------------------------------------------------------
1268 !c Correction of negative values
1269 SUBROUTINE negcor ( X, rho, dz8w, &
1270 ims,ime, jms,jme, kms,kme, & ! memory dims
1271 itimestep, ics, &
1272 its,ite, jts,jte, kts,kte ) ! tile dims
1273 !-----------------------------------------------------------------------
1274 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
1275 INTENT(INOUT) :: X
1276 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
1277 INTENT(IN ) :: rho, dz8w
1278 integer, INTENT(IN ) :: itimestep, ics
1280 !c Local variables
1281 ! REAL, DIMENSION( kts:kte ) :: Y1, Y2
1282 REAL :: A0, A1, A2
1284 A1=0.
1285 A2=0.
1286 do k=kts,kte
1287 do j=jts,jte
1288 do i=its,ite
1289 A1=A1+max(X(i,k,j), 0.)*rho(i,k,j)*dz8w(i,k,j)
1290 A2=A2+max(-X(i,k,j), 0.)*rho(i,k,j)*dz8w(i,k,j)
1291 enddo
1292 enddo
1293 enddo
1295 ! A1=0.0
1296 ! A2=0.0
1297 ! do k=kts,kte
1298 ! A1=A1+Y1(k)
1299 ! A2=A2+Y2(k)
1300 ! enddo
1302 A0=0.0
1304 if (A1.NE.0.0.and.A1.GT.A2) then
1305 A0=(A1-A2)/A1
1307 if (mod(itimestep,540).eq.0) then
1308 if (ics.eq.1) then
1309 write(61,*) 'kms=',kms,' kme=',kme,' kts=',kts,' kte=',kte
1310 write(61,*) 'jms=',jms,' jme=',jme,' jts=',jts,' jte=',jte
1311 write(61,*) 'ims=',ims,' ime=',ime,' its=',its,' ite=',ite
1312 endif
1313 if (ics.eq.1) then
1314 write(61,*) 'qv timestep=',itimestep
1315 write(61,*) ' A1=',A1,' A2=',A2,' A0=',A0
1316 else if (ics.eq.2) then
1317 write(61,*) 'ql timestep=',itimestep
1318 write(61,*) ' A1=',A1,' A2=',A2,' A0=',A0
1319 else if (ics.eq.3) then
1320 write(61,*) 'qr timestep=',itimestep
1321 write(61,*) ' A1=',A1,' A2=',A2,' A0=',A0
1322 else if (ics.eq.4) then
1323 write(61,*) 'qi timestep=',itimestep
1324 write(61,*) ' A1=',A1,' A2=',A2,' A0=',A0
1325 else if (ics.eq.5) then
1326 write(61,*) 'qs timestep=',itimestep
1327 write(61,*) ' A1=',A1,' A2=',A2,' A0=',A0
1328 else if (ics.eq.6) then
1329 write(61,*) 'qg timestep=',itimestep
1330 write(61,*) ' A1=',A1,' A2=',A2,' A0=',A0
1331 else
1332 write(61,*) 'wrong cloud specieis number'
1333 endif
1334 endif
1336 do k=kts,kte
1337 do j=jts,jte
1338 do i=its,ite
1339 X(i,k,j)=A0*AMAX1(X(i,k,j), 0.0)
1340 enddo
1341 enddo
1342 enddo
1343 endif
1345 END SUBROUTINE negcor
1347 SUBROUTINE consat_s ( itaobraun)
1349 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1350 ! c
1351 ! Tao, W.-K., and J. Simpson, 1989: Modeling study of a tropical c
1352 ! squall-type convective line. J. Atmos. Sci., 46, 177-202. c
1353 ! c
1354 ! Tao, W.-K., J. Simpson and M. McCumber, 1989: An ice-water c
1355 ! saturation adjustment. Mon. Wea. Rev., 117, 231-235. c
1357 ! c
1358 ! Tao, W.-K., and J. Simpson, 1993: The Goddard Cumulus Ensemble c
1359 ! Model. Part I: Model description. Terrestrial, Atmospheric and c
1360 ! Oceanic Sciences, 4, 35-72. c
1361 ! c
1362 ! Tao, W.-K., J. Simpson, D. Baker, S. Braun, M.-D. Chou, B. c
1363 ! Ferrier,D. Johnson, A. Khain, S. Lang, B. Lynn, C.-L. Shie, c
1364 ! D. Starr, C.-H. Sui, Y. Wang and P. Wetzel, 2003: Microphysics, c
1365 ! radiation and surface processes in the Goddard Cumulus Ensemble c
1366 ! (GCE) model, A Special Issue on Non-hydrostatic Mesoscale c
1367 ! Modeling, Meteorology and Atmospheric Physics, 82, 97-137. c
1368 ! c
1369 ! Lang, S., W.-K. Tao, R. Cifelli, W. Olson, J. Halverson, S. c
1370 ! Rutledge, and J. Simpson, 2007: Improving simulations of c
1371 ! convective system from TRMM LBA: Easterly and Westerly regimes. c
1372 ! J. Atmos. Sci., 64, 1141-1164. c
1373 ! c
1374 ! Coded by Tao (1989-2003), modified by S. Lang (2006/07) c
1375 ! c
1376 ! Implemented into WRF by Roger Shi 2006/2007 c
1377 ! Additional modifications by Tao, Roger and Steve 2009 c
1378 ! July 25 2010 c
1379 ! Tao November 12 2010 c
1380 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1381 IMPLICIT NONE
1383 ! itaobraun=0 ! see Tao and Simpson (1993)
1384 ! itaobraun=1 ! see Tao et al. (2003)
1386 integer :: itaobraun
1387 real :: cn0
1389 integer k
1390 real :: ga3, ga4, ga5, ga7, ga8, ga9, ga3g2, ga4g2, ga5g2, ga6d
1391 real :: ga3h, ga4h, ga5hh, bc1, dc1, esc, egs, erc, amc, ehs, ehg
1392 real :: ehw, ehi, ehr, sc13, ga6, ga5gh2, egc, cpi2, grvt, tca, dwv
1393 real :: dva, amw, ars, rw, cw, ci, cd1, cd2, ga3b, ga4b, ga6b, ga5bh
1394 real :: ga3g, ga4g, ga5gh, ga3d, ga4d, ga5dh, ac1, ac2, ac3, cc1, eri
1395 real :: ami, ESR, eiw, ui50, ri50, cmn, y1, egi, egr, apri, bpri
1397 !JJS 1/3/2008 vvvvv
1398 !JJS the following common blocks have been moved to the top of
1399 !JJS module_mp_gsfcgce.F
1400 !
1401 ! real, dimension (1:31) :: a1, a2
1402 ! data a1/.7939e-7,.7841e-6,.3369e-5,.4336e-5,.5285e-5,.3728e-5, &
1403 ! .1852e-5,.2991e-6,.4248e-6,.7434e-6,.1812e-5,.4394e-5,.9145e-5, &
1404 ! .1725e-4,.3348e-4,.1725e-4,.9175e-5,.4412e-5,.2252e-5,.9115e-6, &
1405 ! .4876e-6,.3473e-6,.4758e-6,.6306e-6,.8573e-6,.7868e-6,.7192e-6, &
1406 ! .6513e-6,.5956e-6,.5333e-6,.4834e-6/
1407 ! data a2/.4006,.4831,.5320,.5307,.5319,.5249,.4888,.3894,.4047, &
1408 ! .4318,.4771,.5183,.5463,.5651,.5813,.5655,.5478,.5203,.4906, &
1409 ! .4447,.4126,.3960,.4149,.4320,.4506,.4483,.4460,.4433,.4413, &
1410 ! .4382,.4361/
1411 !JJS 1/3/2008 ^^^^^
1413 !23456789012345678901234567890123456789012345678901234567890123456789012
1414 ! ******************************************************************
1415 !JJS
1416 cmin=1.e-20
1417 al = 2.5e10
1418 cp = 1.004e7
1419 rd1 = 1.e-3
1420 rd2 = 2.2
1421 !JJS
1422 cpi=4.*atan(1.)
1423 cpi2=cpi*cpi
1424 grvt=980. !hail
1425 !
1426 c38=3.799052e3
1427 c358=35.86
1428 c610=6.1078e3
1429 c149=1.496286e-5
1430 c879=8.794142
1431 c172=17.26939
1432 c409=4098.026
1433 c76=7.66
1434 c218=21.87456
1435 c580=5807.695
1436 c141=1.414435e7
1438 !
1439 tca=2.43e3
1440 dwv=.226
1441 dva=1.718e-4
1442 amw=18.016
1443 ars=8.314e7
1444 ! scv=2.2904487
1445 !
1446 t0=273.16
1447 t00=238.16
1448 alv=2.5e10 !latent heat vaporization
1449 alf=3.336e9
1450 als=2.8336e10
1451 avc=alv/cp
1452 afc=alf/cp
1453 asc=als/cp
1454 rw=4.615e6
1455 cw=4.187e7
1456 ci=2.093e7
1457 !*** DEFINE THE COEFFICIENTS USED IN TERMINAL VELOCITY
1458 !*** DEFINE THE DENSITY AND SIZE DISTRIBUTION OF PRECIPITATION
1459 !********** HAIL OR GRAUPEL PARAMETERS **********
1461 ! tnw ! rain intercept (1/cm**4)
1462 ! roqr ! rain density (g/cm**3)
1463 ! tns ! snow intercept (1/cm**4)
1464 ! roqs ! snow density (g/cm**3)
1465 ! tng ! graupel intercept (1/cm**4)
1466 ! roqg ! graupel density (g/cm**3)
1468 !********** HAIL PARAMETERS **********
1469 ! tnh = 0.002 !hail medium size
1470 tnh = 0.01
1471 ! tnh = 0.0002 !hail large size
1472 ! tnh = 0.02 !hail small size
1473 roqh = 0.9 !hail
1474 cd1=6.e-1 !hail
1475 cd2=4.*grvt/(3.*cd1) !hail
1476 ah=sqrt(cd2*roqh) !hail
1477 bh=0.5 !hail
1480 !!!!!!!!check /08/08/12!!!!!!!!!!!!
1481 roqg=.3 !bulk density
1482 ! ag=341.84 !for bulk density of 0.4
1483 ag=330.22 !for bulk density of 0.3
1484 ! ag=314.54 !for bulk density of 0.2
1485 bg=.36
1486 ag2=544.83 !for bulk density of 0.5
1487 bg2=.54
1488 tng=.04
1490 roqg2=0.5
1491 qrog2=2.0 !g/m**3
1492 ! qrog2=9.0 !g/m**3 !fixed size test
1493 qrog2=qrog2*1.e-6
1494 ag=330.22
1495 ag2=544.83
1496 bg=0.36
1497 bg2=0.54
1499 !********** SNOW PARAMETERS **********
1500 ! 6/15/02
1501 ! TNS=1.
1502 ! TNS=.08 ! if ice913=1, tao's
1503 tns=.16 ! if ice913=0, tao's snow intercept
1504 roqs=.1 ! snow density g/cm
1505 as=78.63154 ! coefficient a in snow terminal velocity
1506 bs=.11 ! coefficient b in snow terminal velocity
1507 roqs=.05 ! snow density (g/cm**3)
1508 tns=0.1 ! snow intercept (1/cm**4)
1509 as=151.01
1510 bs=0.24
1512 !********** RAIN PARAMETERS **********
1513 aw=2115.
1514 bw=.8
1515 roqr=1. !not defined in Steve's
1516 tnw=.08 !not defined in Steve's
1517 !*****************************************************************
1518 bgh=.5*bg
1519 bgh2=.5*bg2
1520 bsh=.5*bs
1521 bwh=.5*bw
1522 bhh=.5*bh
1523 bgq=.25*bg
1524 bgq2=.25*bg2
1525 bsq=.25*bs
1526 bwq=.25*bw
1527 bhq=.25*bh !hail
1529 !*** define the snow and graupel size mapping parameters for improve=3
1530 tslopes=0.11177209 ! increase tns by 50 from 0 to -35C !sgmap000
1531 tslopeg=0.05756866 ! increase tng by 7.5 from 0 to -35C
1533 draimax=0.0500 !maximum rain diameter (cm)
1534 draimax=draimax**4.*roqr*cpi
1535 dsnomin=0.0100 !minimum snow diameter (cm) !fin16
1536 dsnomin4=dsnomin**4.*0.900*cpi !density at 100 microns
1538 dgrpmin=0.0135 !minimum graupel diameter (cm)
1539 dgrpmin4=dgrpmin**4.*roqg*cpi
1541 xs=0.97
1542 sno11=0.70 !fin15
1543 sno00=-0.24 !fin17
1544 dsno11=3.25 !fin8
1545 dsno00=0.40 !fin12
1546 sexp11=1.5 !sgmapfin17
1547 sexp00=0.9 !fin15
1548 stt=-35. !iterate
1549 stexp=0.42 !fin18
1550 slim=1.00 !fin5
1551 sbase=0.04000 !sgmap00s
1553 xg=0.98
1554 grp11=0.30 !sgmap00x
1555 grp00=0.30 !sgmap00x
1556 dgrp11=2.70 !sgmap001
1557 dgrp00=2.60
1558 gexp11=0.20 !sgmap004
1559 gexp00=0.20 !sgmap004
1560 gtt=-35. !sgmap006 made consistent
1561 gtexp=0.40
1562 glim=0.90
1563 gbase=0.0130 !sgmap00u
1565 hai00=3.25 !g/m**3 fin20
1566 hai11=0.50 !g/m**3 fin20
1567 htt0=-5.00 !deg C fin20
1568 htt1=-50.00 !deg C fin20
1569 haixp=2.7 !exponent fin20
1571 !**********GAMMA FUNCTION CALCULATIONS*************
1572 ga3=2.
1573 ga4=6.
1574 ga5=24.
1575 ga6=120.
1576 ga7=720.
1577 ga8=5040.
1578 ga9=40320.
1579 !
1580 ga3b = gammagce(3.+bw)
1581 ga4b = gammagce(4.+bw)
1582 ga6b = gammagce(6.+bw)
1583 ga5bh = gammagce((5.+bw)/2.)
1584 ga3g = gammagce(3.+bg)
1585 ga4g = gammagce(4.+bg)
1586 ga5gh = gammagce((5.+bg)/2.)
1587 ga3d = gammagce(3.+bs)
1588 ga4d = gammagce(4.+bs)
1589 ga5dh = gammagce((5.+bs)/2.)
1591 ga4g=11.63177
1592 ga3g=3.3233625
1593 ga5gh=1.608355
1594 if(bg.eq.0.37) ga4g=9.730877
1595 if(bg.eq.0.37) ga3g=2.887512
1596 if(bg.eq.0.37) ga5gh=1.526425
1597 if(bg.eq.0.36) ga4g=9.599978
1598 if(bg.eq.0.36) ga3g=2.857136
1599 if(bg.eq.0.36) ga5gh=1.520402
1600 if(bg2.eq.0.54) ga4g2=12.298653
1601 if(bg2.eq.0.54) ga3g2=3.474196
1602 if(bg2.eq.0.54) ga5gh2=1.635061
1603 ga3d=2.54925
1604 ga4d=8.285063
1605 ga5dh=1.456943
1606 if(bs.eq.0.57) ga3d=3.59304
1607 if(bs.eq.0.57) ga4d=12.82715
1608 if(bs.eq.0.57) ga5dh=1.655588
1609 if(bs.eq.0.24) ga3d=2.523508
1610 if(bs.eq.0.24) ga4d=8.176166
1611 if(bs.eq.0.24) ga5dh=1.451396
1612 if(bs.eq.0.11) ga3d=2.218906
1613 if(bs.eq.0.11) ga4d=6.900796
1614 if(bs.eq.0.11) ga5dh=1.382792
1616 ga6d=144.93124
1617 if(bs.eq.0.24) ga6d=181.654791
1618 ga3h=gammagce(3.+bh) !hail
1619 ga4h=gammagce(4.+bh) !hail
1620 ga5hh=gammagce((5.+bh)/2.) !hail
1621 !
1622 !CCCCC LIN ET AL., 1983 OR LORD ET AL., 1984 CCCCCCCCCCCCCCCCC
1623 ac1=aw
1624 ac2=ag ! Steve only defines ac1 and ac2
1625 ac3=as ! need to talk about these 3 parameters.
1627 bc1=bw
1628 cc1=as
1629 dc1=bs
1631 zrc=(cpi*roqr*tnw)**0.25
1632 zsc=(cpi*roqs*tns)**0.25
1633 zgc=(cpi*roqg*tng)**0.25
1634 zgc2=(cpi*roqg2*tng)**0.25
1635 zhc=(cpi*roqh*tnh)**0.25 !hail
1637 vrc=aw*ga4b/(6.*zrc**bw)
1639 vrc0=-26.7
1640 vrc1=20600./zrc
1641 vrc2=-204500./(zrc*zrc)
1642 vrc3=906000./(zrc*zrc*zrc)
1644 vsc=as*ga4d/(6.*zsc**bs)
1645 vgc=ag*ga4g/(6.*zgc**bg)
1646 vgc2=ag2*ga4g2/(6.*zgc2**bg2)
1647 vhc=ah*ga4h/(6.*zhc**bh) !hail
1648 ! ****************************
1649 rn1=9.4e-15
1650 rn2=1.e-3
1651 bnd2=2.0e-3 ! if ice913=0 6/15/02 tao's
1653 esi=0.70
1654 rn3=.25*cpi*tns*as*esi*ga3d
1656 esc=0.45
1657 rn4=.25*cpi*esc*tns*as*ga3d
1659 eri=.1 ! 6/17/02 tao's ice913=0 (not 1)
1660 rn5=.25*cpi*eri*tnw
1661 rn50=-.267e2*ga3
1662 rn51=5.15e3*ga4
1663 rn52=-1.0225e4*ga5
1664 rn53=7.55e3*ga6
1666 ami=1./(24.*6.e-9) ! 6/15/02 tao's
1667 rn6=cpi2*eri*tnw*roqr*ami
1668 rn60=-.267e2*ga6
1669 rn61=5.15e3*ga7
1670 rn62=-1.0225e4*ga8
1671 rn63=7.55e3*ga9
1673 ESR=1. ! also if ice913=1 for tao's
1674 rn7=cpi2*esr*tnw*tns*roqs
1675 rn8=cpi2*esr*tnw*tns*roqr
1677 egs=.1
1678 rn9=cpi2*egs*tns*tng*roqs
1680 rn10=4.*tns
1681 rn101=.65
1682 rn102=.44*sqrt(as/dva)*ga5dh
1683 rn10a=alv*als*amw/(tca*ars)
1684 rn10b=alv/tca
1685 rn10c=ars/(dwv*amw)
1687 rn11=2.*cpi*tns*tca/alf
1688 rn11a=cw/alf
1690 ami50=4.8e-7*(dsnomin*1.e4/2./50.)**3
1691 ami100=1.51e-7 ! improve < 3
1692 ami40=2.46e-7*.875**3 !fin24 35 microns
1694 eiw=1.
1695 ui50=100. ! 6/15/02 tao's
1696 ri50=dsnomin/2.
1698 cmn=1.05e-15
1699 rn12=cpi*eiw*ui50*ri50*ri50
1700 do k=1,31
1701 y1=1.-aa2(k)
1702 rn13(k)=aa1(k)*y1/(ami50**y1-ami40**y1)
1703 rn12a(k)=rn13(k)/ami50
1704 rn12b(k)=aa1(k)*ami50**aa2(k)
1705 rn25a(k)=aa1(k)*cmn**aa2(k)
1707 BergCon1(k)=6.*aa1(k)*ami50**(aa2(k)-1.)
1708 BergCon2(k)=-2.*aa1(k)*ami50**aa2(k)*1.2
1710 BergCon3(k)=6.*aa2(k)/((aa2(k)+1.)*(aa2(k)+2.)) &
1711 *aa1(k)*ami50**(aa2(k)-1.)
1712 BergCon4(k)=2.*(1.-aa2(k))/((aa2(k)+1.)*(aa2(k)+2.)) &
1713 *aa1(k)*ami50**aa2(k)*1.2
1714 enddo
1715 !
1716 egc=0.65
1717 rn14=.25*cpi*egc*ag*tng*ga3g
1718 rn142=.25*cpi*egc*ag2*tng*ga3g2 !high dens graupel
1720 egi=.1
1721 rn15=.25*cpi*egi*tng*ag*ga3g
1722 rn15a=.25*cpi*egi*tng*ag*ga3g
1723 rn15a2=.25*cpi*egi*tng*ag2*ga3g2 !high dens graupel
1724 rn152=.25*cpi*egi*tng*ag2*ga3g2 !high dens graupel
1726 egr=1.
1727 rn16=cpi2*egr*tng*tnw*roqr
1728 rn17=2.*cpi*tng
1729 rn17a2=.31*ga5gh2*sqrt(ag2/dva) !high dens graupel
1730 rn17b=cw-ci
1731 rn17c=cw
1732 rn171=2.*cpi*tng*alv*dwv
1733 rn172=2.*cpi*tng*tca
1734 rn17a=.31*ga5gh*sqrt(ag/dva)
1735 !
1736 apri=.66
1737 bpri=1.e-4
1738 bpri=0.5*bpri ! 6/17/02 tao's
1739 rn18=20.*cpi2*bpri*tnw*roqr
1740 rn18a=apri
1741 rn191=.78 !Are this same as rnn191 (listed below)?
1742 rn192=.31*ga5gh*sqrt(ag/dva) !Are this same as rnn192 (listed below)?
1743 rn192_2=.31*ga5gh2*sqrt(ag2/dva) !high dens graupel
1744 rn19b=cw/alf
1745 rn19=2.*cpi*tng*tca/alf
1746 rn19a=cw/alf
1747 !
1748 rn20=2.*cpi*tng
1749 rn20a=als*als*amw/(tca*ars)
1750 rn20b=als/tca
1751 rn30a=alv*alv*amw/(tca*ars)
1752 rn30=2.*cpi*tng
1754 bnd3=2.e-3
1755 rn21=1.e-3
1756 bnd21=1.5e-3
1757 !
1758 erc=1.
1759 rn22=.25*cpi*erc*tnw
1760 !
1761 rn23=2.*cpi*tnw
1762 rn23a=.31*ga5bh*sqrt(ac1)
1763 rn23b=alv*alv/rw
1764 rn231=.78
1765 rn232=.31*ga3*sqrt(3.e3/dva)
1766 !
1767 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1768 !cc
1769 !cc "c0" in routine "consat" (2d), "consatrh" (3d)
1770 !cc if ( itaobraun.eq.1 ) --> betah=0.5*beta=-.46*0.5=-0.23; cn0=1.e-6
1771 !cc if ( itaobraun.eq.0 ) --> betah=0.5*beta=-.6*0.5=-0.30; cn0=1.e-8
1773 itaobraun=0
1775 cn0=1.e-8
1776 beta=-.6
1777 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1778 !
1779 rn25=cn0
1780 rn30b=alv/tca
1781 rn30c=ars/(dwv*amw)
1782 rn31=1.e-17
1784 rn32=4.*51.545e-4
1785 !!
1786 rn33=4.*tns
1787 rn331=.65
1788 rn332=.44*sqrt(as/dva)*ga5dh
1789 !
1790 amc=1./(24.*4.e-9)
1791 rn34=cpi2*esc*amc*as*roqs*tns*ga6d
1792 rn35=alv*alv/(cp*rw)
1794 gn17=2.*cpi*tng !4ice revised
1795 gn17a=.31*ga5gh*sqrt(ag) !4ice revised
1796 gn17a2=.31*ga5gh2*sqrt(ag2) !4ice revised
1798 ehs=1.
1799 hn9=cpi2*ehs*tns*tnh*roqs
1800 ehg=0.3
1801 hn10=cpi2*ehg*tng*tnh*roqg
1802 ehw=1.
1803 hn14=.25*cpi*ehw*tnh*ga3h*ah
1804 ehi=1.
1805 hn15a=.25*cpi*ehi*tnh*ga3h*ah
1806 ehr=1.
1807 hn16=cpi2*ehr*tnh*tnw*roqr
1808 hn17=2.*cpi*tnh
1809 hn17a=.31*ga5hh*sqrt(ah)
1810 hn19=2.*cpi*tnh/alf
1811 hn19a=.31*ga5hh*sqrt(ah)
1812 hn10a=als*als/rw
1813 hn20=2.*cpi*tnh
1814 hn20b=.31*ga5hh*sqrt(ah)
1816 ! JDC add schmidt number term
1817 sc13 = 0.8420526
1818 rn102 = rn102 * sc13
1819 rn17a = rn17a * sc13
1820 rn17a2 = rn17a2 * sc13
1821 rn192 = rn192 * sc13
1822 rn192_2 = rn192_2 * sc13
1823 rn232 = rn232 * sc13
1824 rn332 = rn332 * sc13
1826 END SUBROUTINE consat_s
1828 !JJS
1829 !JJS REAL FUNCTION GAMMA(X)
1830 !JJS Y=GAMMLN(X)
1831 !JJS GAMMA=EXP(Y)
1832 !JJS RETURN
1833 !JJS END
1834 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1835 !JJS real function GAMMLN (xx)
1836 real function gammagce (xx)
1837 !**********************************************************************
1838 implicit none
1840 real*8 cof(6),stp,half,one,fpf,x,tmp,ser
1841 data cof,stp / 76.18009173,-86.50532033,24.01409822, &
1842 -1.231739516,.120858003e-2,-.536382e-5, 2.50662827465 /
1843 data half,one,fpf / .5, 1., 5.5 /
1844 !
1845 real xx
1846 real gammln
1847 integer j
1849 x=xx-one
1850 tmp=x+fpf
1851 tmp=(x+half)*log(tmp)-tmp
1852 ser=one
1853 do j=1,6
1854 x=x+one
1855 ser=ser+cof(j)/x
1856 enddo !j
1857 gammln=tmp+log(stp*ser)
1858 !JJS
1859 gammagce=exp(gammln)
1860 !JJS
1862 END FUNCTION gammagce
1864 !DIR$ ATTRIBUTES FORCEINLINE :: sgmap
1865 SUBROUTINE sgmap(isg,qcs,qcg,qcg2,qch,qch2,r00,tairc,ftnsg)
1866 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1867 ! compute base snow/graupel intercept scaling factor - ftnsg
1868 ! isg -- flag
1869 ! 1 ftnsg=snow intercept scaling factor
1870 ! 2 ftnsg=graupel intercept scaling factor
1871 ! 3 ftnsg=snow density scaling factor
1872 ! tns is the actual intercept
1873 !
1874 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1875 IMPLICIT NONE
1877 ! common/size/ tnw,tns,tng,roqs,roqg,roqr !defined in the beginning of the module
1878 integer, intent(in) :: isg
1879 ! integer, intent(in) :: i, j, k
1880 ! INTEGER, INTENT(IN) :: ims,ime, jms,jme, kms,kme, &
1881 ! its,ite, jts,jte, kts,kte
1882 real, intent(in) :: r00, tairc
1883 ! real, DIMENSION( ims:ime, jms:jme, kms:kme ), intent(in) :: qcs,qcg, qch
1884 ! real, DIMENSION( its:ite, jts:jte, kts:kte ), intent(in) :: qcs,qcg, qch
1885 real, intent(in) :: qcs,qcg,qcg2,qch,qch2
1886 real, intent(out) :: ftnsg
1889 ! LOCAL variables
1891 !! for snomap
1892 ! real :: xs,sno11,sno00,dsno11,dsno00,sexp11,sexp00,stt, &
1893 ! stexp,sbase,tslopes,dsnomin,slim
1894 !! for grpmap
1895 ! real :: xg,grp11,grp00,dgrp11,dgrp00,gexp11,gexp00,gtt, &
1896 ! gtexp,gbase,tslopeg,dgrpmin, glim
1898 real :: taird, qsg, qsg1, xx, fexp !, cpi, cmin Di
1899 real :: ftnsT, sno1, dsno1, sexp1, ftnsQ, tnsmax, densno
1900 real :: ftngT, grp1, dgrp1, gexp1, ftngQ, tngmax
1901 real :: hx, gx, hgx
1902 real :: kk
1903 real :: hai2, dhai1
1904 ! CPI=4.*ATAN(1.)
1905 ! cmin=1.e-20 !4ice
1907 ftnsg=1.
1909 ! if (isg.eq.1.or.isg.eq.3) qsg=qcs(i,k,j)
1910 ! if (isg.eq.2.) qsg = qcg(i,k,j)
1911 if (isg.eq.1.or.isg.eq.3) qsg=qcs
1912 !NUWRF EMK...Fix options 2 and 4.
1913 if (isg.eq.2) qsg = qcg
1914 if (isg.eq.4) qsg = qch
1916 if (qsg .gt. cmin) then
1918 qsg1=qsg*r00*1.e6
1920 ! if (isg.eq.1) then !snow
1921 if(isg.eq.1.or.isg.eq.3)then !snow 4ice
1922 taird=min(0.,max(stt,tairc)+0.0)
1923 ftnsT=exp(-1.*tslopes*taird)
1924 sno1=sno11
1925 dsno1=dsno11
1926 sexp1=sexp11
1928 if (taird.gt.stt) then
1929 sno1=sno00-(sno00-sno11)*(taird/stt)**stexp
1930 dsno1=dsno00-(dsno00-dsno11)*(taird/stt)**stexp
1931 sexp1=sexp00-(sexp00-sexp11)*(taird/stt)**stexp
1932 endif !taird
1934 hx=0.
1935 gx=0.
1936 hgx=1.
1938 if(qch2*1.e6*r00.gt.0.008) &
1939 hx=max(1.,qch2*1.e6*r00*125.)
1940 if(qcg2*1.e6*r00.gt.0.040) &
1941 gx=max(1.,qcg2*1.e6*r00*25.)
1942 hgx=hx+gx
1943 hgx=max(1.,hgx)
1945 xx=xs-xs*min(slim,max(0.,(qsg1-sno1)/dsno1)**sexp1)
1946 ftnsT=ftnsT**xx
1947 fexp=xx
1948 ftnsQ=1.0
1949 ftnsQ=(qsg1/sbase)**fexp
1950 ftnsg=ftnsT*ftnsQ*hgx
1951 tnsmax=r00*qsg/dsnomin4
1953 if (ftnsg*tns.gt.tnsmax) ftnsg=tnsmax/tns !4ice
1954 ! densno=0.004168*(ftnsg*tns/qsg/r00)**0.3115 !4ice Heymsfield et al. 2004
1955 densno=0.001996*(ftnsg*tns/qsg/r00)**0.2995 !Brandes et al. 2007
1956 if(isg.eq.3) ftnsg=min(densno,0.9)/roqs !4ice
1958 else if (isg.eq.2) then !graupel
1959 taird=min(0.,max(gtt,tairc)+0.0)
1960 ftngT=exp(-1.*tslopeg*taird)
1961 grp1=grp11
1962 dgrp1=dgrp11
1963 gexp1=gexp11
1965 if (taird.gt.gtt) then
1966 grp1=grp00-(grp00-grp11)*(taird/gtt)**gtexp
1967 dgrp1=dgrp00-(dgrp00-dgrp11)*(taird/gtt)**gtexp
1968 gexp1=gexp00-(gexp00-gexp11)*(taird/gtt)**gtexp
1969 endif !taird
1971 xx=xg-xg*min(glim,max(0.0,(qsg1-grp1)/dgrp1)**gexp1)
1972 ftngT=ftngT**xx
1973 fexp=xx
1974 ftngQ=1.0
1975 ftngQ=(qsg1/gbase)**fexp
1976 ftnsg=ftngT*ftngQ
1977 tngmax=r00*qsg/dgrpmin4
1978 if (ftnsg*tng.gt.tngmax) ftnsg=tngmax/tng
1980 elseif(isg.eq.4)then !hail
1981 hai2=hai00
1982 if(tairc.le.htt0.and.tairc.ge.htt1) & !fin18
1983 hai2=hai11+(hai00-hai11)*((tairc-htt1)/(htt0-htt1))**haixp !fin18
1984 if(tairc.le.htt1) hai2=hai11 !fin18
1985 dhai1=hai2 !x2
1986 if(qsg1.ge.hai2) &
1987 ftnsg=1.0-0.80*min(max((qsg1-hai2)/dhai1,0.),1.) ! fin16
1989 endif !isg
1991 endif !qsg
1993 end subroutine sgmap
1995 ! compute fall speed of cloud rain and ice
1996 SUBROUTINE vqrqi(isg,r00,fv,qri,ql,tair,ww1)
1997 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1998 ! compute fall speed of cloud rain and ice
1999 ! isg=1, for rain
2000 ! isg=2, for ice
2001 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2002 implicit none
2004 integer, intent(in) :: isg
2005 real, intent(in) :: r00, fv,qri,ql,tair
2006 real, intent(inout) :: ww1
2008 ! LOCAL variables
2009 integer :: ic
2010 real :: y1,vr,vs,vg
2011 real :: const_vt, const_d, const_m !cpi, cmin Di
2012 real :: bb1, bb2,ice_fall
2013 real :: bin_factor, ftnw, ftnwmin
2014 real, dimension(7) :: aice, vice
2015 data aice/1.e-6, 1.e-5, 1.e-4, 1.e-3, 0.01, 0.1, 1./
2016 data vice/5,15,30,35,40,45,50/
2018 ! cmin=1.e-40
2019 ! CPI=4.*ATAN(1.)
2020 ice_fall=0.
2021 const_vt=1.49e4
2022 const_d=11.9
2023 const_m=1./5.38e7
2025 y1=r00*qri
2026 ww1=0.
2028 if (y1 .gt. cmin) then
2030 if (isg.eq.1) then ! rain
2032 ftnw=1.
2033 if(ql.lt.cmin .and. tair .gt. t0)then
2034 bin_factor=0.11*(1000.*qri)**(-1.27) + 0.98
2035 bin_factor=min(bin_factor, 1.30)
2036 ftnw=1./bin_factor**3.35
2037 ftnwmin=r00*qri/draimax
2038 if(qri.le.0.001) ftnw=max(ftnw,ftnwmin/tnw)
2039 endif
2041 vs=sqrt( y1 )
2042 vg=sqrt( vs)
2043 ! vr=vrc0+vrc1*vg+vrc2*vs+vrc3*vg*vs
2044 vr=vrc0+vrc1*vg/ftnw**0.25+vrc2*vs/ftnw**0.50+vrc3*vg*vs/ftnw**0.75
2045 ww1=max(fv*vr, 0.e0)
2048 else if (isg.eq.2) then ! cloud ice
2050 y1=1.e6*r00*qri ! to g/m**3
2052 if (y1 .gt. 1.e-6) then
2053 y1=y1*1.e-3
2054 bb1=const_m*y1**0.25
2055 bb2=const_d*bb1**0.5
2056 ww1=max(const_vt*bb2**1.31, 0.0)
2057 ww1=ww1*100. !cm/s
2058 if (ww1 .gt. 50.) ww1=50. ! SLang
2059 endif !y1
2060 endif !isg
2061 endif !y1
2063 end subroutine vqrqi
2065 SUBROUTINE saticel_s (dt, dx, itaobraun, istatmin, &
2066 new_ice_sat, id, improve, &
2067 ptwrf, qvwrf, qlwrf, qrwrf, &
2068 qiwrf, qswrf, qgwrf, qhwrf, &
2069 rho_mks, pi_mks, p0_mks, w_mks, &
2070 itimestep, xland, &
2071 refl_10cm, diagflag, do_radar_ref, & ! GT added for reflectivity calcs
2072 ids,ide, jds,jde, kds,kde, &
2073 ims,ime, jms,jme, kms,kme, &
2074 its,ite, jts,jte, kts,kte, &
2075 !NUWRF BEGIN
2076 re_cloud_gsfc, re_rain_gsfc, re_ice_gsfc, &
2077 re_snow_gsfc, re_graupel_gsfc, re_hail_gsfc, & ! cloud effective radius
2078 physc, physe, physd, physs, physm, physf, &
2079 acphysc, acphyse, acphysd, acphyss, acphysm, acphysf &
2080 #if ( WRF_CHEM == 1)
2081 !JJS 20110525 vvvvv
2082 ,aero, icn_diag, nc_diag, gid, &
2083 !JJS 20110525 ^^^^^
2084 ! EMK
2085 chem_opt, &
2086 gsfcgce_gocart_coupling &
2087 #endif
2088 ,ii,j,irestrict)
2089 !NUWRF END
2091 !-----------------------------------------------------------------------
2092 ! USE module_dm
2093 IMPLICIT NONE
2094 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2095 ! c
2096 ! History: c
2097 ! c
2098 ! Coded by Tao (1989-2003), modified by S. Lang (2006/07) c
2099 ! c
2100 ! Implemented into WRF by Jainn Shi 2006/2007 c
2101 ! Improved by S. Lang (2008-2009) c
2102 ! Implemented by Tao, Jainn Shi and tested by Jainn Shi c
2103 ! Modified by Tao, Jul. 2010 c
2104 ! Modified by Tao, Aug. 2010 c
2105 ! Added 4ICE code (from S. Lang) into NU-WRF by Di Wu in 2014 c
2106 ! Added aerosol coupling by Jainn Shi, Mar. 2014 c
2107 ! Added cloud droplet eff. radius code by Jainn Shi, Dec. 2014 c
2108 ! Code optimization by J. Mielikainen (2016) c
2109 ! c
2110 ! References: c
2111 ! c
2112 ! Tao, W.-K., and J. Simpson, 1989: Modeling study of a tropical c
2113 ! squall-type convective line. J. Atmos. Sci., 46, 177-202. c
2114 ! c
2115 ! Tao, W.-K., J. Simpson and M. McCumber, 1989: An ice-water c
2116 ! saturation adjustment. Mon. Wea. Rev., 117, 231-235. c
2117 ! c
2118 ! c
2119 ! Tao, W.-K., and J. Simpson, 1993: The Goddard Cumulus Ensemble c
2120 ! Model. Part I: Model description. Terrestrial, Atmospheric and c
2121 ! Oceanic Sciences, 4, 35-72. c
2122 ! c
2123 ! Tao, W.-K., J. Simpson, D. Baker, S. Braun, M.-D. Chou, B. c
2124 ! Ferrier,D. Johnson, A. Khain, S. Lang, B. Lynn, C.-L. Shie, c
2125 ! D. Starr, C.-H. Sui, Y. Wang and P. Wetzel, 2003: Microphysics, c
2126 ! radiation and surface processes in the Goddard Cumulus Ensemble c
2127 ! (GCE) model, A Special Issue on Non-hydrostatic Mesoscale c
2128 ! Modeling, Meteorology and Atmospheric Physics, 82, 97-137. c
2129 ! c
2130 ! Lang, S., W.-K. Tao, R. Cifelli, W. Olson, J. Halverson, S. c
2131 ! Rutledge, and J. Simpson, 2007: Improving simulations of c
2132 ! convective system from TRMM LBA: Easterly and Westerly regimes. c
2133 ! J. Atmos. Sci., 64, 1141-1164. c
2134 ! c
2135 ! Tao, W.-K., J. J. Shi, S. Lang, C. Peters-Lidard, A. Hou, S. c
2136 ! Braun, and J. Simpson, 2007: New, improved bulk-microphysical c
2137 ! schemes for studying precipitation processes in WRF. Part I: c
2138 ! Comparisons with other schemes. c
2139 ! c
2140 ! Lang, S., W.-K. Tao, X. Zeng, and Y. Li, 2011: Reducing the Biases in c
2141 ! Simulated Radar Reflectivities from a Bulk Microphysics Scheme: c
2142 ! Tropical Convective Systems. J. Atmos. Sci., 68, 2306-2320. c
2143 ! c
2144 ! Shi, J. J., T. Matsui, W.-K. Tao, C. Peters-Lidard, M. Chin, Q. Tan, c
2145 ! K. Pickering, N. Guy, S. Lang, and E. Kemp., 2014: Implementation of c
2146 ! an Aerosol-Cloud Microphysics-Radiation Coupling into the NASA c
2147 ! Unified WRF: Simulation Results for the 6-7 August 2006 AMMA Special c
2148 ! Observing Period. Quart. J. Roy. Meteor. Soc., 140, 2158-2175, c
2149 ! doi:10.1002/qj.2286. c
2150 ! c
2151 ! Stephen E. Lang, Wei-Kuo Tao, Jiun-Dar Chern, Di Wu, c
2152 ! and Xiaowen Li, 2014: Benefits of a Fourth Ice Class c
2153 ! in the Simulated Radar Reflectivities of Convective Systems c
2154 ! Using a Bulk Microphysics Scheme. J. Atmos. Sci., 71, 3583-3612. c
2155 ! doi: http://dx.doi.org/10.1175/JAS-D-13-0330.1 c
2156 ! c
2157 ! Lang, S., W.-K. Tao, J.-D. Chern, D. Wu, and X. Li, 2014: c
2158 ! Benefits of a 4th ice class in the simulated radar reflectivities c
2159 ! of convective systems using a bulk microphysics scheme. J. Atmos. c
2160 ! Sci., 71, 3583-3612. doi: http://dx.doi.org/10.1175/JAS-D-13-0330.1 c
2161 ! c
2162 ! Tao, W.-K., D. Wu, S. Lang, J.-D. Chern, C. Peters-Lidard, c
2163 ! A. Fridlind, and T. Matsui (2016), High-resolution NU-WRF c
2164 ! simulations of a deep convective-precipitation system during c
2165 ! MC3E: Further improvements and comparisons between Goddard c
2166 ! microphysics schemes and observations, J. Geophys. Res. Atmos., c
2167 ! 121, 1278-1305, doi:10.1002/2015JD023986. c
2168 ! c
2169 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2170 !
2171 ! COMPUTE ICE PHASE MICROPHYSICS AND SATURATION PROCESSES
2172 !
2174 !cc using scott braun's way for pint, pidep computations
2175 INTEGER, INTENT(IN ) :: itaobraun, improve, new_ice_sat
2176 integer, intent(in) :: id
2177 integer, intent(in) :: itimestep,istatmin
2178 real, intent(in) :: dt ! timestep (second)
2179 real, intent(in) :: dx ! grid resolution (meters)
2180 real :: thresh_evap
2181 !cc
2183 !JJS 20090623 vvvvv
2184 integer, intent(in) :: ids,ide,jds,jde,kds,kde
2185 integer, intent(in) :: ims,ime,jms,jme,kms,kme
2186 integer, intent(in) :: its,ite,jts,jte,kts,kte
2187 integer, intent(in) :: ii,j,irestrict ! global i-index inside local i-loops: ii+i-1
2188 integer i, k, kp
2190 real, dimension(CHUNK) :: afcp, alvr, ascp, avcp, rp0, pi0, pir, &
2191 pr0, r00, rrs, rrq, fv0, fvs, cp409, &
2192 rr0, zrr, zsr, zgr, zhr, &
2193 cp580, cs580, cv409, vscf, vgcf, vgcf2, &
2194 vhcr, dwvp, r3f, r4f, r5f, r6f, &
2195 r12r, r14f ,r14f2, r15af, r15af2, r15f, r18r, &
2196 r22f, r25rt, r32rt, r331r, r332rf, &
2197 r34f
2199 real :: bg3, bg3_2, bgh5, bgh5_2, bs3 ,bs6, bsh5, bw3 ,bw6 ,bwh5, &
2200 cmin, cmin1, cmin2, d2t, del, f2 ,f3, ft, qb0, r25a, r_nci, &
2201 sccc, sddd, seee, sfff, smmm, ssss, tb0, temp, ucog ,ucog2, &
2202 ucor ,ucos, ucoh, uwet, rdt, bnd1, c_nci, &
2203 r10t, r20t, r23t
2205 real :: a_1, a_2, a_3, a_4
2206 real :: a_11, a_22, a_33, a_44
2207 real :: zdry, zwet, zwet0
2208 real :: vap_frac
2209 !JJS 20090623 ^^^^^
2211 real, dimension (CHUNK, kts:kte) :: fv
2212 real, dimension (CHUNK, kts:kte) :: dpt, dqv
2213 real, dimension (CHUNK, kts:kte) :: qcl, qrn, &
2214 qci, qcs, qcg, qch
2215 real, dimension (CHUNK, kts:kte) :: qsz, qgz,qhz
2216 !JJS
2217 real, dimension (CHUNK, kms:kme), INTENT(INOUT) &
2218 :: ptwrf, qvwrf, &
2219 qlwrf, qrwrf, &
2220 qiwrf, qswrf, &
2221 qgwrf, qhwrf
2223 !JJS in MKS
2224 real, dimension (CHUNK, kms:kme), INTENT(IN ) &
2225 :: rho_mks, &
2226 pi_mks, &
2227 p0_mks, &
2228 w_mks
2229 !JJS COMMON /BADV/
2230 real, dimension (CHUNK) :: &
2231 vg, zg, &
2232 ps, pg, &
2233 prn, psn, &
2234 pwacs, wgacr, &
2235 pidep, pint, &
2236 qsi, ssi, &
2237 esi, esw, &
2238 qsw, pr, &
2239 ssw, pihom, &
2240 pidw, pimlt, &
2241 psaut, qracs, &
2242 psaci, psacw, &
2243 qsacw, praci, &
2244 pmlts, pmltg, &
2245 asss
2247 !JJS COMMON/BSAT/
2248 real, dimension (CHUNK) :: &
2249 praut, pracw, &
2250 psfw, psfi, &
2251 dgacs, dgacw, &
2252 dgaci, dgacr, &
2253 pgacs, wgacs, &
2254 qgacw, wgaci, &
2255 qgacr, pgwet, &
2256 pgaut, pracs, &
2257 psacr, qsacr, &
2258 pgfr, psmlt, &
2259 pgmlt, psdep, &
2260 pgdep, piacr, &
2261 egs
2263 !JJS COMMON/BSAT1/
2264 real, dimension (CHUNK) :: &
2265 pt, qv, &
2266 qc, qr, &
2267 qi, qs, &
2268 qg, qh, &
2269 tair, &
2270 tairc, rtair, &
2271 dep, dd, &
2272 dd1, qvs, &
2273 dm, rq, &
2274 rsub1, col, &
2275 cnd, ern, &
2276 dlt1, dlt2, &
2277 dlt3, dlt4, &
2278 zr, vr, &
2279 zs, vs
2282 !JJS COMMON/BSAT2H/
2283 real, dimension (CHUNK) :: &
2284 phfr,phmlt, & !4ice
2285 dhacw,qhacw,dhacr,qhacr,whacr, & !4ice
2286 dhaci,whaci,dhacs,phacs,whacs, & !4ice
2287 dhacg,whacg,phwet,phdep,phsub, & !4ice
2288 pvaph,primh,scv,dwv,tca !4ice
2290 !JJS COMMON/B5/
2291 real, dimension (CHUNK,kts:kte) :: rho !only in satice in cgs
2293 !JJS COMMON/B6/
2294 real, dimension (CHUNK, kts:kte) :: p0, pi, f0, ww1
2295 real, dimension (CHUNK, kts:kte) :: &
2296 fd, fe, &
2297 st, sv, &
2298 sq, sc, &
2299 se, sqa
2301 !JJS COMMON/BI/ IT(its:ite,jts:jte), ICS(its:ite,jts:jte,4)
2302 integer, dimension (CHUNK) :: it
2303 integer, dimension (CHUNK, 4) :: ics
2305 integer :: i24h
2306 real :: r2is, r2ig, r2ih
2309 !JJS COMMON/MICRO/
2310 real, dimension (CHUNK, kms:kme), INTENT(INOUT) :: &
2311 physc, physe, physd, &
2312 physs, physm, physf, &
2313 acphysc, acphyse, acphysd, &
2314 acphyss, acphysm, acphysf
2316 ! EMK NUWRF
2317 real, dimension(CHUNK, kts:kte) :: dbz
2319 !JJS 9/30/2009 for Steve's new improvement
2321 integer :: ihalmos
2322 real :: xnsplnt, xmsplnt
2323 real :: hmtemp1, hmtemp2, hmtemp3, hmtemp4
2324 real :: ftnw, ftnwmin
2325 real :: xssi, fssi, rssi, xsubi, wssi
2326 real :: dmicrons, dmicrong, dvair, alpha
2327 real, dimension (CHUNK) :: tairN, tairI, &
2328 ftns, ftng, &
2329 ftns0, ftng0, &
2330 ftnh, ftnh0, &
2331 pihms, pihmg, &
2332 pimm, pcfr, &
2333 pssub, pgsub, &
2334 fros, fros0, &
2335 vi, zh, &
2336 vh, pihmh, &
2337 dda0, pvapg, &
2338 pracg,qracg, &
2339 qrimh,pg2h !4ice revised
2341 real, dimension (CHUNK) :: y1, y2, y3, y4, &
2342 y5, y6, y7, y8
2343 ! for Xiping's new dbz code
2344 real :: hfact, sfact, yy1
2345 real :: xncld, esat, rv, rlapse_m
2346 real :: delT, bhi !Di deleted cpi
2347 real :: rc, ra, cna
2348 real :: xccld, xknud, cunnf, diffar
2349 real :: qgz2, qhz2
2351 real :: r11t, r19t, r19at, r30t, r33t
2352 real, dimension(CHUNK) :: r7rf, r8rf, r9rf, r16rf
2353 real, dimension(CHUNK) :: r101r, r102rf, r191r, r192rf, r192rf2
2354 real, dimension(CHUNK) :: r231r, r232rf
2355 real, dimension(CHUNK) :: h9r, h10r, h14r, h15ar, h16r, h17r, h17aq, &
2356 h19aq, h19rt, h10ar, h20t, h20bq
2357 real, dimension(CHUNK) :: bin_factor, rim_frac
2358 real :: term1, term2, fdwv, dwv0 !JDC water vapor diffusivity correction term
2359 integer :: iter
2361 !NUWRF BEGIN
2363 #if ( WRF_CHEM == 1)
2364 ! JJS 20110525 vvvvv
2365 ! for inline Gocart coupling
2366 INTEGER, INTENT(IN ) :: gid
2367 INTEGER, PARAMETER :: num_go = 14 ! number of the gocart aerosol species
2368 REAL, DIMENSION( CHUNK, kms:kme, num_go), intent(in) :: aero
2369 REAL, DIMENSION( CHUNK, kms:kme ), intent(out) :: icn_diag !IN concentration [#/Litre]
2370 REAL, DIMENSION( CHUNK, kms:kme ), intent(out) :: nc_diag !cloud concentration [#/cm3]
2371 integer,intent(in) :: chem_opt ! EMK
2372 integer,intent(in) :: gsfcgce_gocart_coupling ! EMK
2374 ! Local Variables
2376 real :: e_sat, e_dry !saturated and dry air water vapor [hPa, mb]
2377 real :: rh_rad ! relative humidity [%]
2378 real :: super_sat !super saturation [%]
2379 real :: ccn_out(CHUNK) ! CCN conc [#/cm3] ! EMK TEST
2380 real :: icn_out(CHUNK) ! IN conc [#/Litter] ! EMK TEST
2381 real :: P_liu_daum ! autoconversion rate [g/cm3 s-1] !
2382 real :: re_liu_daum ! effective radius of cloud [micron] !
2383 real,parameter :: min_icn = 0.01 !minimum # conc of IN [#/Litre]
2385 ! JJS 20110525 ^^^^^
2386 #endif
2388 !JJS 20140226 variables for the calculation of effective radius of cloud species
2389 real, dimension (CHUNK, kms:kme) , INTENT(INOUT ) &
2390 :: re_cloud_gsfc, re_rain_gsfc, &
2391 re_ice_gsfc, re_snow_gsfc, &
2392 re_graupel_gsfc, re_hail_gsfc
2393 REAL , DIMENSION( CHUNK ) , INTENT(IN) :: XLAND
2394 real, parameter :: roqi = 0.9179 ! ice density
2395 real, parameter :: ccn_over_land = 1500 ! [#/cm3] climatological value
2396 real, parameter :: ccn_over_water = 150 ! [#/cm3] climatological value
2397 real :: L_cloud ! cloud water [g/cm3] !
2398 real :: I_cloud ! cloud water [g/cm3] !
2399 real :: mu, ccn_ref, lambda
2400 real :: gamfac1, gamfac3
2401 !JJS 20140226 ^^^^^
2402 !NUWRF END
2404 !+---+-----------------------------------------------------------------+
2405 REAL, DIMENSION(CHUNK, kms:kme), INTENT(INOUT):: refl_10cm ! GT
2407 LOGICAL, OPTIONAL, INTENT(IN) :: diagflag
2408 INTEGER, OPTIONAL, INTENT(IN) :: do_radar_ref
2409 !+---+-----------------------------------------------------------------+
2411 ! JDC dwv0 is water vapor diffusivity at STP
2412 dwv0 = 0.226
2414 !JJS20090623 save
2416 if (itimestep.eq.1) then
2417 do k = kts, kte
2418 !dir$ vector aligned
2419 DO i=1,irestrict
2420 physc(i,k)=0.
2421 physe(i,k)=0.
2422 physd(i,k)=0.
2423 physs(i,k)=0.
2424 physf(i,k)=0.
2425 physm(i,k)=0.
2426 acphysc(i,k)=0.
2427 acphyse(i,k)=0.
2428 acphysd(i,k)=0.
2429 acphyss(i,k)=0.
2430 acphysf(i,k)=0.
2431 acphysm(i,k)=0.
2432 ENDDO
2433 enddo !k
2434 if ( wrf_dm_on_monitor() .and. i.eq.its .and. j.eq.jts ) then
2435 write(6, *) ' latent heating variables have been initialized to 0. at timestep = ', itimestep
2436 endif
2437 endif
2439 !JJS convert from mks to cgs, and move from WRF grid to GCE grid
2440 do k=kts,kte
2441 !dir$ vector aligned
2442 DO i=1,irestrict
2443 rho(i,k)=rho_mks(i,k)*0.001
2444 p0(i,k)=p0_mks(i,k)*10.0
2445 pi(i,k)=pi_mks(i,k)
2446 ww1(i,k)=w_mks(i,k)*100.
2447 dpt(i,k)=ptwrf(i,k)
2448 dqv(i,k)=qvwrf(i,k)
2449 qcl(i,k)=qlwrf(i,k)
2450 qrn(i,k)=qrwrf(i,k)
2451 qci(i,k)=qiwrf(i,k)
2452 qcs(i,k)=qswrf(i,k)
2453 qcg(i,k)=qgwrf(i,k)
2454 qch(i,k)=qhwrf(i,k)
2455 ENDDO
2456 enddo !k
2458 do k=kts,kte
2459 !dir$ vector aligned
2460 DO i=1,irestrict
2461 fv(i,k)=sqrt(rho(i,1)/rho(i,k))
2462 ENDDO
2463 enddo !k
2464 !JJS
2466 !
2467 ! ****** THREE CLASSES OF ICE-PHASE (LIN ET AL, 1983) *********
2469 d2t=dt
2471 r2ig=1.
2472 r2is=1.
2473 r2ih=1.
2475 !C TAO 2007 END
2477 cmin=1.e-20
2478 !JJS 10/7/2008
2479 cmin1=1.e-20
2480 cmin2=1.e-35
2482 xssi=0.16 ! maximum allowable ice supersaturation from SEL
2483 rssi=0.05 ! maximum allowable residual ice supersaturation SEL
2484 xsubi=0.70 ! ice RH below which qi must sublimate
2486 ! HALLET-MOSSOP RIME SPLINTERING parameters
2487 ihalmos=1
2488 xnsplnt=370. ! peak # splinters per milligram of rime
2489 xmsplnt=4.4e-8 ! mass of a splinter (from Ferrier 1994)
2490 hmtemp1=-2.
2491 hmtemp2=-4.
2492 hmtemp3=-6.
2493 hmtemp4=-8.
2495 it(:)=1
2496 !!!!!!!!08/07/12
2497 f2=rd1*d2t
2498 f3=rd2*d2t
2500 ft=dt/d2t
2502 bw3=bw+3.
2503 bs3=bs+3.
2504 bg3=bg+3.
2505 bh3=bh+3
2506 bsh5=2.5+bsh
2507 bgh5=2.5+bgh
2508 bhh5=2.5+bhh
2509 bwh5=2.5+bwh
2510 bw6=bw+6.
2511 bs6=bs+6.
2512 ! betah=.5*beta
2514 r10t=rn10*d2t
2515 r25a=rn25
2517 rdt=1./d2t
2518 r11t=rn11*d2t
2519 r19t=rn19*d2t
2520 r19at=rn19a*d2t
2521 r20t=rn20*d2t
2522 r23t=rn23*d2t
2523 r30t=rn30*d2t
2524 r33t=rn33*d2t
2525 bg3_2=bg2+3
2526 bgh5_2=2.5+bgh2
2528 !JJS 20150831
2529 ! Calculate the threshold for reducing spurious evaporation
2530 ! a function of dx (grid resolution in km)
2532 thresh_evap = -39.974 * exp(-1.194 * dx/1000.)
2534 if ( wrf_dm_on_monitor() .and. itimestep.eq.1 .and. &
2535 i.eq.its .and. j.eq.jts ) then
2536 print *,'GSFCGCE 4ice scheme inside satice improve=',improve
2537 print *,'dx, thresh_evap = ', dx, thresh_evap
2538 print *,'no reduce suprious evaporation adjustment'
2539 endif
2541 Rc=1.e-3 ! cloud droplet radius 10 microns
2542 Ra=1.e-5 ! aerosol radius 0.1 microns
2543 Cna=500. ! contact nuclei conc per cc
2544 Bhi=1.01e-2 ! pollen (Deihl et al. 2006)
2546 !C ******************************************************************
2548 do 1000 k=kts,kte
2549 kp=k+1
2550 tb0=0.
2551 qb0=0.
2553 !dir$ vector aligned
2554 DO i=1,irestrict
2556 ! EMK BUG FIX...Initialize variable
2557 #if ( WRF_CHEM == 1)
2558 ccn_out(i) = 0.e0
2559 icn_out(i) = 0.e0
2560 if ( (chem_opt == 300 .or. chem_opt == 301 .or. &
2561 chem_opt == 302 .or. chem_opt == 303) .and. &
2562 (gsfcgce_gocart_coupling == 1) ) then
2563 icn_diag(i,k) = icn_out(i) ! #/Litre
2564 nc_diag(i,k) = ccn_out(i) ! #/cm3
2565 else
2566 icn_diag(i,k) = 0.
2567 nc_diag(i,k) = 0.
2568 end if
2569 #endif
2570 ! EMK END
2572 rp0(i)=3.799052e3/p0(i,k)
2573 pi0(i)=pi(i,k)
2574 pir(i)=1./(pi(i,k))
2575 pr0(i)=1./p0(i,k)
2576 r00(i)=rho(i,k)
2577 rr0(i)=1./rho(i,k)
2578 rrs(i)=sqrt(rr0(i))
2579 rrq(i)=sqrt(rrs(i))
2580 f0(i,k)=al/cp/pi(i,k)
2581 fv0(i)=fv(i,k)
2582 fvs(i)=sqrt(fv(i,k))
2583 zrr(i)=1.e5*zrc*rrq(i)
2584 zsr(i)=1.e5*zsc*rrq(i)
2585 zgr(i)=1.e5*zgc*rrq(i)
2586 zhr(i)=1.e5*zhc*rrq(i)
2587 cp409(i)=c409*pi0(i)
2588 cv409(i)=c409*avc !
2589 cp580(i)=c580*pi0(i)
2590 cs580(i)=c580*asc !
2591 alvr(i)=r00(i)*alv
2592 afcp(i)=afc*pir(i)
2593 avcp(i)=avc*pir(i)
2594 ascp(i)=asc*pir(i)
2595 vscf(i)=vsc*fv0(i)
2596 vgcf(i)=vgc*fv0(i)
2597 vgcf2(i)=vgc2*fv0(i)
2598 vhcr(i)=vhc*rrs(i)
2599 dwvp(i)=c879*pr0(i)
2601 r3f(i)=rn3*fv0(i)
2602 r4f(i)=rn4*fv0(i)
2603 r5f(i)=rn5*fv0(i)
2604 r6f(i)=rn6*fv0(i)
2606 r12r(i)=rn12*r00(i)
2607 r14f(i)=rn14*fv0(i)
2608 r14f2(i)=rn142*fv0(i)
2609 r15f(i)=rn15*fv0(i)
2610 r15af(i)=rn15a*fv0(i) !4ice revised
2611 r15af2(i)=rn15a2*fv0(i) !4ice revised
2612 r18r(i)=rn18*rr0(i)
2613 r22f(i)=rn22*fv0(i)
2614 r25rt(i)=rn25*rr0(i)*d2t
2615 r32rt(i)=rn32*d2t*rrs(i)
2618 !JJS added 10/1/2009
2619 r7rf(i)=rn7*rr0(i)*fv0(i)
2620 r8rf(i)=rn8*rr0(i)*fv0(i)
2621 r9rf(i)=rn9*rr0(i)*fv0(i)
2622 r16rf(i)=rn16*rr0(i)*fv0(i)
2623 r101r(i)=rn101*rr0(i)
2624 r102rf(i)=rn102*rrs(i)*fvs(i)
2625 r191r(i)=rn191*rr0(i)
2626 r192rf(i)=rn192*rrs(i)*fvs(i)
2627 r192rf2(i)=rn192_2*rrs(i)*fvs(i)
2628 r331r(i)=rn331*rr0(i)
2629 r332rf(i)=rn332*rrs(i)*fvs(i)
2630 r34f(i)=rn34*fv0(i)
2632 r231r(i)=rn231*rr0(i)
2633 r232rf(i)=rn232*rrs(i)*fvs(i)
2635 h9r(i)=hn9*rr0(i)
2636 h10r(i)=hn10*rr0(i)
2637 h14r(i)=hn14*rrs(i)
2638 h15ar(i)=hn15a*rrs(i)
2639 h16r(i)=hn16*rr0(i)
2640 h17r(i)=hn17*rr0(i)
2641 h17aq(i)=hn17a*rrq(i)
2642 h19aq(i)=hn19a*rrq(i)
2643 h19rt(i)=hn19*rr0(i)*d2t
2644 h10ar(i)=hn10a*r00(i)
2645 h20t(i)=hn20*d2t !
2646 h20bq(i)=hn20b*rrq(i)
2648 pt(i)=dpt(i,k)
2649 qv(i)=dqv(i,k)
2650 qc(i)=qcl(i,k)
2651 qr(i)=qrn(i,k)
2652 qi(i)=qci(i,k)
2653 qs(i)=qcs(i,k)
2654 qg(i)=qcg(i,k)
2655 qh(i)=qch(i,k)
2656 if (qc(i) .le. cmin) qc(i)=0.0
2657 if (qr(i) .le. cmin) qr(i)=0.0
2658 if (qi(i) .le. cmin) qi(i)=0.0
2659 if (qs(i) .le. cmin) qs(i)=0.0
2660 if (qg(i) .le. cmin) qg(i)=0.0
2661 if (qh(i) .le. cmin) qh(i)=0.0
2662 tair(i)=(pt(i)+tb0)*pi0(i)
2663 tairc(i)=tair(i)-t0
2664 zr(i)=zrr(i)
2665 zs(i)=zsr(i)
2666 zg(i)=zgr(i)
2667 zh(i)=zhr(i)
2668 vr(i)=0.0
2669 vs(i)=0.0
2670 vg(i)=0.0
2671 vi(i)=0.0
2672 vh(i)=0.0
2674 ftns(i)=1.
2675 ftng(i)=1.
2676 ftns0(i)=1.
2677 ftng0(i)=1.
2678 fros(i)=1.
2679 ftnh(i)=1.
2680 ftnh0(i)=1.
2682 cnd(i)=0.0
2683 dep(i)=0.
2684 ern(i)=0.0
2685 pint(i)=0.0
2686 pidep(i)=0.0
2688 psdep(i)=0.
2689 pgdep(i)=0.
2690 dd1(i)=0.
2691 dd(i)=0.
2692 pgsub(i)=0.
2693 psmlt(i)=0.
2694 pgmlt(i)=0.
2695 pimlt(i)=0.
2696 psacw(i)=0.
2697 piacr(i)=0.
2699 pssub(i)=0.0
2700 pgsub(i)=0.0
2702 psfw(i)=0.0
2703 psfi(i)=0.0
2704 pidep(i)=0.0
2706 pgfr(i)=0.
2707 psacr(i)=0.
2708 wgacr(i)=0.
2709 pihom(i)=0.
2710 pidw(i)=0.0
2712 psaut(i)=0.0
2713 psaci(i)=0.0
2714 praci(i)=0.0
2715 pwacs(i)=0.0
2716 qsacw(i)=0.0
2718 pracs(i)=0.0
2719 qracs(i)=0.0
2720 qsacr(i)=0.0
2721 pgaut(i)=0.0
2723 praut(i)=0.0
2724 pracw(i)=0.0
2725 pgfr(i)=0.0
2727 qracs(i)=0.0
2729 pgacs(i)=0.0
2730 qgacw(i)=0.0
2731 dgaci(i)=0.0
2732 dgacs(i)=0.0
2733 wgacs(i)=0.0
2734 wgaci(i)=0.0
2735 dgacw(i)=0.0
2736 dgacr(i)=0.
2737 pgwet(i)=0.0
2739 qgacr(i)=0.0
2741 pihom(i)=0.0
2742 pimlt(i)=0.0
2743 pidw(i)=0.0
2744 pimm(i)=0.0
2745 pcfr(i)=0.0
2747 pihms(i)=0.0
2748 pihmg(i)=0.0
2749 ftns(i)=1.
2750 ftng(i)=1.
2751 pmlts(i)=0.0
2752 pmltg(i)=0.0
2754 phfr(i)=0.0
2755 phmlt(i)=0.0
2756 dhacw(i)=0.0
2757 qhacw(i)=0.0
2758 dhacr(i)=0.0
2759 qhacr(i)=0.0
2760 whacr(i)=0.0
2761 dhaci(i)=0.0
2762 whaci(i)=0.0
2763 dhacs(i)=0.0
2764 phacs(i)=0.0
2765 whacs(i)=0.0
2766 dhacg(i)=0.0
2767 whacg(i)=0.0
2768 phwet(i)=0.0
2769 phdep(i)=0.0
2770 phsub(i)=0.0
2771 pvapg(i)=0.0
2772 pvaph(i)=0.0
2773 primh(i)=0.0
2775 dlt4(i)=0.0
2776 dlt3(i)=0.0
2777 dlt2(i)=0.0
2779 ! ******************************************************************
2780 ! *** Y1 : DYNAMIC VISCOSITY OF AIR (U)
2781 ! *** DWV : DIFFUSIVITY OF WATER VAPOR IN AIR (PI)
2782 ! *** TCA : THERMAL CONDUCTIVITY OF AIR (KA)
2783 ! *** Y2 : KINETIC VISCOSITY (V)
2785 y1(i)=c149*tair(i)**1.5/(tair(i)+120.)
2786 dwv(i)=dwvp(i)*tair(i)**1.81
2787 tca(i)=c141*y1(i)
2788 scv(i)=1./((rr0(i)*y1(i))**.1666667*dwv(i)**.3333333)
2790 ENDDO
2791 !dir$ vector aligned
2792 DO i=1,irestrict
2794 !JJS for calculating processes related to both ice and warm rain
2796 ! *** COMPUTE ZR,ZS,ZG,VR,VS,VG *****************************
2798 if (qr(i) .gt. cmin) then
2799 dd(i)=r00(i)*qr(i)
2800 y1(i)=sqrt(dd(i))
2801 y2(i)=sqrt(y1(i))
2802 zr(i)=zrc/y2(i)
2803 endif
2805 call vqrqi(1,r00(i),fv0(i),qr(i),qc(i),tair(i),vr(i))
2806 call vqrqi(2,r00(i),fv0(i),qi(i),qc(i),tair(i),vi(i))
2808 ftns(i)=1.
2809 ftns0(i)=1.
2811 qhz2=qhwrf(i,k)
2812 qgz2=qgwrf(i,k)
2813 if (k .lt. kte-2 .and. tairc(i) .ge. -5) then
2814 qhz2=qhwrf(i,k+1)
2815 qgz2=qgwrf(i,k+1)
2816 endif
2817 call sgmap(1,qs(i),qg(i),qgz2,qh(i),qhz2,r00(i),tairc(i),ftns0(i))
2818 call sgmap(3,qs(i),qg(i),qgz2,qh(i),qhz2,r00(i),tairc(i),fros0(i))
2820 if (qs(i) .gt. cmin) then
2821 dd(i)=r00(i)*qs(i)
2822 y1(i)=dd(i)**.25
2823 ftns(i)=1.
2824 ftns(i)=ftns0(i)**0.25
2825 fros(i)=1 !improve4
2826 fros(i)=fros0(i)**0.25 !improve4
2827 ZS(i)=ZSC/Y1(i)*ftns(i)*fros(i) !improve4
2828 ftns(i)=ftns0(i)**bsq
2829 fros(i)=fros0(i)**bsq !improve4
2830 VS(i)=MAX(vscf(i)*DD(i)**BSQ/ftns(i)/fros(i), 0.)
2831 endif
2833 ftng(i)=1.
2834 ftng0(i)=1.
2835 call sgmap(2,qs(i),qg(i),qgz2,qh(i),qhz2,r00(i),tairc(i),ftng0(i))
2837 if (qg(i) .gt. cmin) then
2838 dd(i)=r00(i)*qg(i)
2839 y1(i)=dd(i)**.25
2840 ftng(i)=1.
2841 ftng(i)=ftng0(i)**0.25
2843 zg(i)=zgc/y1(i)*ftng(i)
2844 if(dd(i).gt.qrog2) zg(i)=zgc2/y1(i)*ftng(i)
2846 ftng(i)=ftng0(i)**bgq
2847 vg(i)=max(vgcf(i)*dd(i)**bgq/ftng(i), 0.0)
2848 if(dd(i).gt.qrog2)then
2849 ftng(i)=ftng0(i)**bgq2
2850 vg(i)=max(vgcf2(i)*dd(i)**bgq2/ftng(i), 0.e0)
2851 endif !improve4
2852 endif !qg
2854 call sgmap(4,qs(i),qg(i),qgz2,qh(i),qhz2,r00(i),tairc(i),ftnh0(i))
2856 if (qh(i) .gt. cmin) then
2857 dd(i)=r00(i)*qh(i)
2858 y1(i)=dd(i)**.25
2859 ftnh(i)=ftnh0(i)**0.25
2860 zh(i)=zhc/y1(i)*ftnh(i)
2861 ftnh(i)=ftnh0(i)**bhq
2862 vh(i)=max(vhcr(i)*dd(i)**bhq/ftnh(i), 0.0)
2863 endif
2865 if (qr(i) .le. cmin1) vr(i)=0.0
2866 if (qs(i) .le. cmin1) vs(i)=0.0
2867 if (qg(i) .le. cmin1) vg(i)=0.0
2868 if (qi(i) .le. cmin1) vi(i)=0.0
2869 if (qh(i) .le. cmin1) vh(i)=0.0
2871 !! ******************************************************************
2872 !! *** Y1 : DYNAMIC VISCOSITY OF AIR (U)
2873 !! *** DWV : DIFFUSIVITY OF WATER VAPOR IN AIR (PI)
2874 !! *** TCA : THERMAL CONDUCTIVITY OF AIR (KA)
2875 !! *** Y2 : KINETIC VISCOSITY (V)
2876 !
2878 !* 1 * PSAUT : AUTOCONVERSION OF QI TO QS ***1**
2879 !* 3 * PSACI : ACCRETION OF QI TO QS ***3**
2880 !* 4 * PSACW : ACCRETION OF QC BY QS (RIMING) (QSACW FOR PSMLT) ***4**
2881 !* 5 * PRACI : ACCRETION OF QI BY QR ***5**
2882 !* 6 * PIACR : ACCRETION OF QR OR QG BY QI ***6**
2883 !* 34 * pwacs : collection of qs by qc **34**
2885 pihms(i)=0.0
2886 pihmg(i)=0.0
2887 psaut(i)=0.0
2888 psaci(i)=0.0
2889 praci(i)=0.0
2890 piacr(i)=0.0
2891 psacw(i)=0.0
2892 pwacs(i)=0.0
2893 qsacw(i)=0.0
2894 ftns(i)=1.
2895 ftng(i)=1.
2896 ftnh(i)=1.
2897 ftns(i)=ftns0(i)
2898 ftng(i)=ftng0(i)
2900 if (tair(i).lt.t0) then
2902 rn1=1./300.
2903 bnd1=6.e-5
2904 esi(i)=0.25
2905 psaut(i)=r2is*max(rn1*esi(i)*(qi(i)-bnd1*fv0(i)*fv0(i)) ,0.0)
2906 ftns(i)=ftns0(i)
2907 ftng(i)=ftng0(i)
2908 fros(i)=fros0(i)
2909 ftnh(i)=ftnh0(i)
2910 esi(i)=1.0 !esi constant in consatrh via r3f/rn3; use esi( ) to make f(T)/f(q)
2911 dmicrons=(r00(i)*qs(i)/(roqs*fros(i))/cpi/(tns*ftns(i)))**.25*1.e4
2912 esi(i)=min(1.,(dmicrons/375.)**3.) ! f(dmicrons)
2914 y1(i)=1.0
2915 if (vs(i).gt.0.) y1(i)=abs((vs(i)-vi(i)) &
2916 /vs(i))
2917 psaci(i)=r2is*y1(i)*r3f(i)*qi(i)/zs(i)**bs3*ftns(i)*esi(i)
2918 if(qs(i).le.cmin) psaci(i)=0.
2919 psacw(i)=r2is*r4f(i)*qc(i)/zs(i)**bs3*ftns(i)
2920 if(qs(i).le.cmin) psacw(i)=0.
2921 if (ihalmos.eq.1)then
2922 y2(i)=0.
2923 if((tairc(i).le.hmtemp1).and.(tairc(i).ge.hmtemp4)) &
2924 y2(i)=0.5
2925 if((tairc(i).le.hmtemp2).and.(tairc(i).ge.hmtemp3)) &
2926 y2(i)=1.
2927 pihms(i)=r2ih*psacw(i)*y2(i)*xnsplnt*1000.*xmsplnt
2928 psacw(i)=psacw(i)-pihms(i)
2929 endif
2930 pwacs(i)=r2is*r34f(i)*qc(i)/zs(i)**bs6*ftns(i)*fros(i) !improve4
2931 if(qs(i).le.cmin) pwacs(i)=0.
2932 y1(i)=1./zr(i)
2933 y2(i)=y1(i)*y1(i)
2934 y3(i)=y1(i)*y2(i)
2935 y5(i)=1.0
2936 if (vr(i).gt.0.) y5(i)=abs((vr(i)-vi(i)) &
2937 /vr(i))
2938 dd(i)=y5(i)*r5f(i)*qi(i)*y3(i)*(rn50+rn51*y1(i) &
2939 +rn52*y2(i)+rn53*y3(i))
2940 praci(i)=max(dd(i),0.0)
2941 if (qr(i) .le. cmin) praci(i)=0.
2942 y4(i)=y3(i)*y3(i)
2943 dd1(i)=y5(i)*r6f(i)*qi(i)*y4(i)*(rn60+rn61*y1(i) &
2944 +rn62*y2(i)+rn63*y3(i))
2946 piacr(i)=max(dd1(i),0.0)
2947 if (qr(i) .le. cmin) piacr(i)=0.
2948 else
2949 qsacw(i)=r2is*r4f(i)*qc(i)/zs(i)**bs3*ftns(i)
2950 if (qs(i) .le. cmin) psacw(i)=0.
2951 endif !tairc
2954 !23456789012345678901234567890123456789012345678901234567890123456789012
2955 !* 21 * PRAUT AUTOCONVERSION OF QC TO QR **21**
2956 !* 22 * PRACW : ACCRETION OF QC BY QR **22**
2958 #if (WRF_CHEM == 1)
2959 !JJS 20110602 vvvvv
2960 ! EMK...Only execute when GOCART and coupling are selected
2961 if ( (chem_opt == 300 .or. chem_opt == 301 .or. &
2962 chem_opt == 302 .or. chem_opt == 303) .and. &
2963 (gsfcgce_gocart_coupling == 1) ) then
2964 ! sat vapor pressure [hPa,mb]
2965 e_sat = 6.11 * exp( 5423. *( 1.0/273.15 - 1./tair(i) ) )
2966 ! dry air vapor pressure [hPa, mb]
2967 e_dry = qv(i) / ( qv(i) + 0.622 ) * p0(i,k) * 1.e-3 ! p0 in [g*cm/s2/cm2]
2968 rh_rad = max(1.e-6, e_dry/e_sat*100.) ! relative humidity [%]
2969 super_sat = max(0.001, rh_rad - 100.e0) !super saturation [%]
2970 !
2971 ! convert gocart aerosol mass conc to CN
2972 !
2973 call mass2ccn(tair(i),super_sat,aero(i,k,:),ccn_out(i) )
2974 ! nc_cgs(i,k) = max(100., ccn_out) !diagnostic cloud droplet conc [#/cm3]
2975 ccn_out(i) = max(100., ccn_out(i))
2977 ! rho_dryair = p0(i,k) / ( tair(i) * 2.87 * 1.0e6) ! dry air density [g/cm3]
2978 L_cloud = qc(i) * rho(i,k) ! cloud water [g/cm3]
2979 ! g/g g/cm3
2980 ! call auto_conversion( L_cloud, nc_cgs(i,k), P_liu_daum, re_liu_daum )
2981 call auto_conversion( L_cloud, ccn_out(i), P_liu_daum, re_liu_daum )
2983 praut(i) = P_liu_daum / rho(i,k) !autoconversion rate [g/g s-1]
2984 else
2985 praut(i)=max(rn21*(qc(i)-bnd21),0.0)
2986 end if ! if (gsfcgce_gocart_coupling == 1)
2987 #else
2988 praut(i)=max(rn21*(qc(i)-bnd21),0.0)
2989 !JJS 20110602 ^^^^^
2990 #endif
2992 y1(i)=1./zr(i)
2993 y2(i)=y1(i)*y1(i)
2994 y3(i)=y1(i)*y2(i)
2995 y4(i)=r22f(i)*qc(i)*y3(i)*(rn50+rn51*y1(i)+ &
2996 rn52*y2(i)+rn53*y3(i))
2997 pracw(i)=max(y4(i), 0.0)
2998 if(qr(i) .le. cmin) pracw(i)=0.
3000 !* 12 * PSFW : BERGERON PROCESSES FOR QS (KOENING, 1971) **12**
3001 !* 13 * PSFI : BERGERON PROCESSES FOR QS **13**
3003 pidep(i)=0.0
3004 psfw(i)=0.0
3005 psfi(i)=0.0
3007 ENDDO
3008 !dir$ vector aligned
3009 DO i=1,irestrict
3011 if (tair(i) .lt. t0) then
3012 y1(i)=max( min(tairc(i), -1.), -31.)
3013 it(i)=int(abs(y1(i)))
3014 y1(i)=rn12a(it(i))
3015 y2(i)=rn12b(it(i))
3016 y3(i)=rn13(it(i))
3017 psfw(i)=r2is*max(d2t*y1(i)*(y2(i)+r12r(i)*qc(i))* &
3018 qi(i), 0.0)
3019 psfi(i)=r2is*y3(i)*qi(i)
3020 !
3021 y4(i)=1./(tair(i)-c358)
3022 y5(i)=1./(tair(i)-c76)
3023 qsw(i)=rp0(i)*exp(c172-c409*y4(i))
3024 qsi(i)=rp0(i)*exp(c218-c580*y5(i))
3025 hfact=(qv(i)+qb0-qsi(i))/(qsw(i)-qsi(i)+cmin1)
3026 ! add cmin1 to prevent overflow of hfact
3027 if(hfact.gt.1.) hfact=1.
3028 sfact=1
3029 SSI(i)=(qv(i)+qb0)/qsi(i)-1.
3031 fssi=min(xssi,max(rssi,xssi*(tairc(i)+44.)/(44.0-38.0))) !improve4 max ssi f(tair)
3032 fssi=rssi
3033 wssi=ww1(i,k)/100.-2.0
3034 if(wssi.gt.0.) fssi=rssi+min(xssi,wssi*0.01)
3035 fssi=min(ssi(i),fssi)
3037 ! STEVE : PLEASE CHECK
3038 if (tairc(i).le.-5.) then !meyers
3039 r_nci=max(1.e-3*exp(-.639+12.96*fssi),0.528e-3) !meyers et al.
3040 else
3041 r_nci=min(1.e-3*exp(-.639+12.96*fssi),0.528e-3) !meyers et al.
3042 endif !tairc
3044 if (r_nci.gt.15.) r_nci=15. !cap at 15000/liter
3046 ! Cooper curve
3047 if( tairc(i) .lt. -40.0 ) then
3048 c_nci = 5.0e-6*exp(0.304*40.0)
3049 else
3050 c_nci = 5.0e-6*exp(0.304*abs(tairc(i)))
3051 end if
3052 r_nci = c_nci
3053 ! JDC cap r_nci to the amount corresponding to crystal size of ami50
3054 r_nci = max(r_nci,r00(i)*qi(i)/ami50)
3056 #if (WRF_CHEM == 1)
3057 ! EMK...Only execute when GOCART and coupling is turned on.
3058 if ( (chem_opt == 300 .or. chem_opt == 301 .or. &
3059 chem_opt == 302 .or. chem_opt == 303) .and. &
3060 (gsfcgce_gocart_coupling == 1) ) then
3061 !JJS 20110602 vvvvv
3062 ! Conversion rate of cloud water to ice in the Bergeron porcess based on Meyer + DeMott formulae
3063 !
3064 ! convert gocart aerosol mass conc to IN
3065 ! p0 need to be converted from g*cm/s2/cm2 to mb (hPa)
3066 ! call mass2icn(p0(i,k)*0.001,tair(i),aero(i,k,:), icn_out)
3067 ! call mass2icn(p0(i,k)*0.001,tair(i),aero(i,k,:), icn_out,i,k)
3068 ! EMK...mass2icn requires i,j,k indices, but only prints
3069 ! the values when an error occurs. In no case are the i,j,k
3070 ! used to access a value in an array.
3071 ! We will simply pass a bogus value of j=1 in the call
3072 ! below.
3073 call mass2icn(p0(i,k)*0.001,tair(i),aero(i,k,:), icn_out(i),&
3074 i,1,k)
3076 icn_out(i) = min(1.e3, max(0.01e0 , icn_out(i)) )
3077 r_nci = icn_out(i) * 1.e-3 !DeMotto's formuale
3078 !JJS 20110602 ^^^^^
3079 end if ! if (gsfcgce_gocart_coupling == 1)
3080 #endif
3082 dd(i)=min((r00(i)*qi(i)/r_nci), ami40) !mean cloud ice mass
3083 yy1=1.-aa2(it(i))
3084 sfact=(AMI50**YY1-AMI40**YY1)/(AMI50**YY1-dd(i)**YY1)
3086 !JDC water vapor diffusivity correction term
3087 esi(i) = qsi(i)/rp0(i)*c610
3088 y3(i) = 1./tair(i)
3089 ! term1 = y3(i)*(rn10a*y3(i)-rn10b)
3090 term1 = y3(i)*(rn20a*y3(i)-rn20b)
3091 term2 = rn10c*tair(i)/esi(i)
3092 dd(i) = term1+term2
3093 fdwv = dd(i)/(term1+term2*dwv0/dwv(i))
3094 psfw(i) = max( d2t*y1(i)*fdwv*(y2(i)*fdwv &
3095 + r12r(i)*qc(i))*qi(i),0.0 )
3097 psfw(i)=0.0
3099 if (hfact.gt.0.) then
3100 psfi(i)=r2is*psfi(i)*hfact*sfact*fdwv
3101 else
3102 psfi(i)=0.
3103 endif !hfact
3105 if(qi(i).le.1.e-5*fv0(i)*fv0(i)) psfi(i)=0.0
3107 endif !tair
3109 ENDDO
3110 !dir$ vector aligned
3111 DO i=1,irestrict
3113 !TTT***** QG=QG+MIN(PGDRY,PGWET)
3114 !* 9 * PGACS : ACCRETION OF QS BY QG (DGACS,WGACS: DRY AND WET) ***9**
3115 !* 14 * DGACW : ACCRETION OF QC BY QG (QGACW FOR PGMLT) **14**
3116 !* 16 * DGACR : ACCRETION OF QR TO QG (QGACR FOR PGMLT) **16**
3117 !*******PGDRY : DGACW+DGACI+DGACR+DGACS ******
3118 !* 15 * DGACI : ACCRETION OF QI BY QG (WGACI FOR WET GROWTH) **15**
3119 !* 17 * PGWET : WET GROWTH OF QG **17**
3120 !* Steve turned off PGWET, set PGWET = 0.
3121 !* Steve turned off wgaci, set wgaci = 0.
3123 y1(i)=abs( vg(i)-vs(i) )
3124 y2(i)=zs(i)*zg(i)
3125 y3(i)=5./y2(i)
3126 y4(i)=.08*y3(i)*y3(i)
3127 y5(i)=.05*y3(i)*y4(i)
3128 y2(i)=y1(i)*(y3(i)/zs(i)**5+y4(i)/zs(i)**3 &
3129 +y5(i)/zs(i))
3131 pgacs(i)=r2ig*r2is*r9rf(i)*y2(i)*ftns(i)*ftng(i)*fros(i) !improve4
3132 if(qs(i).le.cmin) pgacs(i)=0.
3133 if(qg(i).le.cmin) pgacs(i)=0.
3134 dgacs(i)=pgacs(i)
3135 dgacs(i)=0.0 !Lang et al. 2007
3136 !crh wgacs(i)=10.*r9rf(i)*y2(i)
3137 wgacs(i)=0.0
3138 wgacs(i)=10.*r9rf(i)*y2(i)*ftns(i)*ftng(i)*fros(i) !4ice revised
3139 ! if(r00(i)*qg(i).gt.qrog2)then !4ice revised
3140 ! wgacs(i)=10.*r9rf(i)*y2(i)*ftns(i)*ftng(i)*fros(i) !4ice revised
3141 ! endif !4ice revised
3142 if(qg(i).le.cmin) wgacs(i)=0. !4ice revised2
3143 if(qs(i).le.cmin) wgacs(i)=0. !4ice revised2
3145 dhacs(i)=0.0
3146 whacs(i)=0.0
3147 y1(i)=abs( vh(i)-vs(i) )
3148 y2(i)=zs(i)*zh(i)
3149 y3(i)=5./y2(i)
3150 y4(i)=.08*y3(i)*y3(i)
3151 y5(i)=.05*y3(i)*y4(i)
3152 dd(i)=Y1(i)*(Y3(i)/ZS(i)**5+Y4(i)/ZS(i)**3 &
3153 +Y5(i)/ZS(i))
3154 whacs(i)=r2ih*r2is*min(h9r(i)*dd(i)*ftnh(i)*ftns(i)*fros(i), &
3155 qs(i)/d2t)
3156 if(qs(i).le.cmin) whacs(i)=0.
3157 if(qh(i).le.cmin) whacs(i)=0.
3159 dhacg(i)=0.0
3160 whacg(i)=0.0
3161 y1(i)=abs( vh(i)-vg(i) )
3162 y2(i)=zg(i)*zh(i)
3163 y3(i)=5./y2(i)
3164 y4(i)=.08*y3(i)*y3(i)
3165 y5(i)=.05*y3(i)*y4(i)
3166 dd(i)=Y1(i)*(Y3(i)/ZG(i)**5+Y4(i)/ZG(i)**3 &
3167 +Y5(i)/ZG(i))
3168 whacg(i)=r2ih*r2ig*min(h10r(i)*dd(i)*ftnh(i)*ftng(i), &
3169 qg(i)/d2t)
3170 if(r00(i)*qg(i).gt.qrog2) &
3171 whacg(i)=whacg(i)/roqg*roqg2*0.5 !reduce ehg for high dens grp
3172 if(qg(i).le.cmin) whacg(i)=0.
3173 if(qh(i).le.cmin) whacg(i)=0.
3175 y1(i)=1./zg(i)**bg3
3176 esi(i)=1.0 !egc constant in consatrh via r14f/rn14; use esi( ) to make f(T)/f(q)
3177 dmicrong=(r00(i)*qg(i)/roqg/cpi/(tng*ftng(i)))**.25*1.e4
3178 if(r00(i)*qg(i).gt.qrog2)then
3179 y1(i)=1./zg(i)**bg3_2
3180 dmicrong=(r00(i)*qg(i)/roqg2/cpi/(tng*ftng(i)))**.25*1.e4
3181 endif
3182 esi(i)=min(1.,(dmicrong/500.)**1.1) ! f(dmicrons)
3184 dgacw(i)= r2ig*esi(i)*r14f(i)*qc(i)*y1(i)*ftng(i)
3185 if(r00(i)*qg(i).gt.qrog2) &
3186 dgacw(i)=r2ig*esi(i)*r14f2(i)*qc(i)*y1(i)*ftng(i)
3187 if(qg(i).le.cmin) dgacw(i)=0.
3189 dhacw(i)=0.0 !4ice
3190 y2(i)=1./zh(i)**bh3 !4ice
3191 dhacw(i)=r2ih*max(h14r(i)*qc(i)*y2(i)*ftnh(i), 0.0) !4ice
3192 if(qh(i).le.cmin) dhacw(i)=0.
3194 if(ihalmos.eq.1)then
3195 y2(i)=0.
3196 if ((tairc(i).le.hmtemp1).and.(tairc(i).ge.hmtemp4)) &
3197 y2(i)=0.5
3198 if ((tairc(i).le.hmtemp2).and.(tairc(i).ge.hmtemp3)) &
3199 y2(i)=1.
3200 pihmg(i)=r2ig*dgacw(i)*y2(i)*xnsplnt*1000.*xmsplnt
3201 dgacw(i)=r2ig*(dgacw(i)-pihmg(i))
3202 pihmh(i)=0.0 !4ice
3203 pihmh(i)=r2ih*dhacw(i)*y2(i)*xnsplnt*1000.*xmsplnt !4ice
3204 dhacw(i)=r2ih*(dhacw(i)-pihmh(i)) !4ice
3205 endif !ihalmos
3207 qgacw(i)=r2ig*dgacw(i)
3208 qhacw(i)=r2ih*dhacw(i) !4ice
3209 y1(i)=1./zg(i)**bg3
3210 y5(i)=1.0
3211 if (vg(i).gt.0.) y5(i)=abs((vg(i)-vi(i))/vg(i))
3212 dgaci(i)= r2ig*y5(i)*r15f(i)*qi(i)*y1(i)*ftng(i)
3213 if(qg(i).le.cmin) dgaci(i)=0.
3214 dgaci(i)=0.0
3215 wgaci(i)=0.0
3216 wgaci(i)=y5(i)*r15af(i)*qi(i)*y1(i)*ftng(i) !4ice revised
3217 if(r00(i)*qg(i).gt.qrog2)then !4ice revised
3218 y1(i)=1./zg(i)**bg3_2
3219 wgaci(i)=y5(i)*r15af2(i)*qi(i)*y1(i)*ftng(i) !4ice revised
3220 endif !4ice revised
3221 if(qg(i).le.cmin) wgaci(i)=0. !4ice revised2
3223 dhaci(i)=0.0 !4ice
3224 whaci(i)=0.0 !4ice
3225 y5(i)=1.0 !4ice
3226 if(vh(i).gt.0.) y5(i)=abs((vh(i)-vi(i))/vh(i)) !4ice
3227 y2(i)=1./zh(i)**bh3 !4ice
3228 whaci(i)=r2ih*min(y5(i)*h15ar(i)*qi(i)*y2(i)*ftnh(i), &
3229 qi(i)/d2t) !4ice
3230 if(qh(i).le.cmin) whaci(i)=0.
3232 y1(i)=abs( vg(i)-vr(i) )
3233 y2(i)=zr(i)*zg(i)
3234 y3(i)=5./y2(i)
3235 y4(i)=.08*y3(i)*y3(i)
3236 y5(i)=.05*y3(i)*y4(i)
3237 dd(i)=r16rf(i)*y1(i)*(y3(i)/zr(i)**5+y4(i)/zr(i)**3 &
3238 +y5(i)/zr(i))*ftng(i)
3239 dgacr(i)=r2ig*max(dd(i),0.0)
3240 if(qg(i).le.cmin) dgacr(i)=0.
3241 if(qr(i).le.cmin) dgacr(i)=0.
3242 qgacr(i)=dgacr(i)
3244 pracg(i)=0.
3245 qracg(i)=0.
3246 y2(i)=zr(i)*zg(i)
3247 y3(i)=5./y2(i)
3248 y4(i)=.08*y3(i)*y3(i)
3249 y5(i)=.05*y3(i)*y4(i)
3250 pracg(i)=r16rf(i)/roqr*roqg*y1(i)*(y3(i)/zg(i)**5+y4(i) &
3251 /zg(i)**3+y5(i)/zg(i))*ftng(i)
3252 if(r00(i)*qg(i).gt.qrog2) pracg(i)=pracg(i)/roqg*roqg2
3253 if(qg(i).le.cmin) pracg(i)=0.
3254 if(qr(i).le.cmin) pracg(i)=0.
3255 qracg(i)=min(d2t*pracg(i), qg(i))
3257 y1(i)=abs( vh(i)-vr(i) )
3258 y2(i)=zr(i)*zh(i)
3259 y3(i)=5./y2(i)
3260 y4(i)=.08*y3(i)*y3(i)
3261 y5(i)=.05*y3(i)*y4(i)
3262 DD(i)=h16r(i)*Y1(i)*ftnh(i)*(Y3(i)/ZR(i)**5 &
3263 +Y4(i)/ZR(i)**3+Y5(i)/ZR(i))
3264 dhacr(i)=r2ih*max(dd(i), 0.0)
3265 if(qh(i).le.cmin) dhacr(i)=0.
3266 if(qr(i).le.cmin) dhacr(i)=0.
3267 qhacr(i)=dhacr(i)
3269 if (tair(i) .ge. t0) then
3270 dgacs(i)=0.0
3271 wgacs(i)=0.0 !4ice revised
3272 whacs(i)=0.0
3273 whacg(i)=0.0
3274 dgacw(i)=0.0
3275 dhacw(i)=0.0 !ice4
3276 dgaci(i)=0.0
3277 wgaci(i)=0.0 !4ice revised
3278 whaci(i)=0.0 !ice4
3279 dgacr(i)=0.0
3280 pracg(i)=0.0
3281 dhacr(i)=0.0
3282 else
3283 pgacs(i)=0.0
3284 qgacw(i)=0.0
3285 qhacw(i)=0.0 !ice4
3286 qgacr(i)=0.0
3287 qracg(i)=0.0
3288 qhacr(i)=0.0
3289 endif
3291 PGWET(i)=0.0
3292 if(tair(i) .lt. t0)then !4ice revised
3294 ! NUWRF corrected by Steve and Roger to prevent division by 0
3295 y1(i)=1./(alf+rn17c*max(tairc(i),-75.0))
3296 y2(i)=r191r(i)/zg(i)**2+r192rf(i)/zg(i)**bgh5 !4ice revised chern
3297 if(r00(i)*qg(i).gt.qrog2) & !4ice revised chern
3298 y2(i)=(r191r(i)/zg(i)**2+r192rf2(i)/zg(i)**bgh5_2) !4ice revised chern
3299 Y4(i)=ALVR(i)*DWV(i)*(RP0(i)-(QV(i)+QB0))-TCA(i)*TAIRC(i) !4ice revised
3300 DD(i)=Y1(i)*(rn20*ftng(i)*Y4(i)*Y2(i)+(WGACI(i) & !4ice revised chern
3301 +WGACS(i))*(ALF+RN17B*TAIRC(i))) !4ice revised
3302 PGWET(i)=max(DD(i), 0.0) !4ice revised
3303 if(qg(i).le.cmin) pgwet(i)=0. !4ice revised
3304 endif !4ice revised
3308 phwet(i)=0.0
3309 if (tair(i) .lt. t0) then
3310 y1(i)=1./(alf+rn17c*tairc(i))
3311 y3(i)=.78/zh(i)**2+h17aq(i)*scv(i)/zh(i)**bhh5
3312 Y4(i)=ALVR(i)*DWV(i)*(RP0(i)-(QV(i)+QB0))-TCA(i)*TAIRC(i)
3313 DD(i)=Y1(i)*(h17r(i)*ftnh(i)*Y4(i)*Y3(i)+(WHACI(i)+WHACS(i) &
3314 +WHACG(i))*(ALF+RN17B*TAIRC(i)))
3315 phwet(i)=r2ih*max(DD(i), 0.0)
3316 if(qh(i).le.cmin) phwet(i)=0.
3317 endif !tair
3319 ENDDO
3320 !dir$ vector aligned
3321 DO i=1,irestrict
3323 !******** HANDLING THE NEGATIVE CLOUD WATER (QC) ******************
3324 y1(i)=qc(i)/d2t
3325 psacw(i)=min(y1(i), psacw(i))
3326 pihms(i)=min(y1(i), pihms(i))
3327 praut(i)=min(y1(i), praut(i))
3328 pracw(i)=min(y1(i), pracw(i))
3329 psfw(i)= min(y1(i), psfw(i))
3330 dgacw(i)=min(y1(i), dgacw(i))
3331 pihmg(i)=min(y1(i), pihmg(i))
3332 dhacw(i)=min(y1(i), dhacw(i)) !4ice
3333 pihmh(i)=min(y1(i), pihmh(i)) !4ice
3334 qsacw(i)=min(y1(i), qsacw(i))
3335 qgacw(i)=min(y1(i), qgacw(i))
3336 qhacw(i)=min(y1(i), qhacw(i)) !4ice
3338 y1(i)=d2t*(psacw(i)+praut(i)+pracw(i)+psfw(i) &
3339 +dgacw(i)+qsacw(i)+qgacw(i)+pihms(i)+pihmg(i) &
3340 +dhacw(i)+qhacw(i)+pihmh(i)) !4ice
3342 qc(i)=qc(i)-y1(i)
3343 !
3344 if (qc(i) .lt. 0.0) then
3345 y2(i)=1.
3346 if (y1(i) .ne. 0.) y2(i)=qc(i)/y1(i)+1.
3347 psacw(i)=psacw(i)*y2(i)
3348 praut(i)=praut(i)*y2(i)
3349 pracw(i)=pracw(i)*y2(i)
3350 psfw(i)=psfw(i)*y2(i)
3351 dgacw(i)=dgacw(i)*y2(i)
3352 dhacw(i)=dhacw(i)*y2(i)
3353 qsacw(i)=qsacw(i)*y2(i)
3354 qgacw(i)=qgacw(i)*y2(i)
3355 qhacw(i)=qhacw(i)*y2(i)
3356 pihms(i)=pihms(i)*y2(i)
3357 pihmg(i)=pihmg(i)*y2(i)
3358 pihmh(i)=pihmh(i)*y2(i) !4ice
3359 qc(i)=0.0
3360 endif
3362 ! !4ice revised
3363 ! convert wet graupel to hail !4ice revised
3364 ! !4ice revised
3365 y1(i)=dgacw(i)+dgacr(i) !4ice revised
3366 if(Y1(i).lt..95*pgwet(i).or.y1(i).eq.0.)THEN !4ice revised
3367 wgaci(i)=0.0 !4ice revised
3368 wgacs(i)=0.0 !4ice revised
3369 endif !4ice revised
3370 pg2h(i)=0.0 !4ice revised
3371 if(y1(i).gt.0..and.pgwet(i).gt.0..and. & !4ice revised
3372 Y1(i).gt.1.0*pgwet(i))THEN !4ice revised
3373 !test pg2h(i)=2.0*y1(i) !4ice revised
3374 pg2h(i)=qg(i)/d2t !4ice revised
3375 ! pg2h(i)=qg(i)*0.647232/d2t !4ice revised Deff
3376 pg2h(i)=min(pg2h(i), qg(i)/d2t) !4ice revised
3377 endif !4ice revised
3379 whacr(i)=phwet(i)-dhacw(i)-whaci(i)-whacs(i)-whacg(i) !4ice
3380 y2(i)=dhacw(i)+dhacr(i) !4ice
3381 if(y2(i).lt.0.95*phwet(i).or.y2(i).eq.0.) THEN
3382 whacr(i)=0.0 !4ice
3383 whaci(i)=0.0 !4ice
3384 whacg(i)=0.0
3385 whacs(i)=0.0
3386 endif
3388 primh(i)=0.0
3389 if(y2(i).gt.0..and.phwet(i).gt.0..and. Y2(i).lt..95*phwet(i)) THEN
3391 rim_frac(i)=2.0*min((1.-y2(i)/phwet(i))**2,1.0)
3392 rim_frac(i)=rim_frac(i)*min((tairc(i)/(t00-t0))**2,1.0)
3393 ! primh(i)=rim_frac(i)*y2(i)
3394 primh(i)=rim_frac(i)*dhacw(i)
3395 primh(i)=min(primh(i), qh(i)/d2t)
3396 endif
3398 ENDDO
3399 !dir$ vector aligned
3400 DO i=1,irestrict
3402 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3403 !******** HANDLING THE NEGATIVE CLOUD ICE (QI) ******************
3404 y1(i)=qi(i)/d2t
3405 psaut(i)=min(y1(i), psaut(i))
3406 psaci(i)=min(y1(i), psaci(i))
3407 praci(i)=min(y1(i), praci(i))
3408 psfi(i)= min(y1(i), psfi(i))
3409 dgaci(i)=min(y1(i), dgaci(i))
3410 wgaci(i)=min(y1(i), wgaci(i))
3411 whaci(i)=min(y1(i), whaci(i)) !4ice
3413 qi(i)=qi(i)+d2t*(pihms(i)+pihmg(i)+pihmh(i)) !fix SEL
3415 y1(i)=d2t*(psaut(i)+psaci(i)+praci(i)+psfi(i) &
3416 +dgaci(i)+wgaci(i)+whaci(i))
3418 qi(i)=qi(i)-y1(i)
3419 !
3420 if (qi(i) .lt. 0.0) then
3421 y2(i)=1.
3422 if (y1(i) .ne. 0.0) y2(i)=qi(i)/y1(i)+1.
3423 psaut(i)=psaut(i)*y2(i)
3424 psaci(i)=psaci(i)*y2(i)
3425 praci(i)=praci(i)*y2(i)
3426 psfi(i)=psfi(i)*y2(i)
3427 dgaci(i)=dgaci(i)*y2(i)
3428 wgaci(i)=wgaci(i)*y2(i)
3429 whaci(i)=whaci(i)*y2(i) !4ice
3430 qi(i)=0.0
3431 endif
3433 wgacr(i)=qgacr(i)+qgacw(i)
3434 dlt3(i)=0.0
3435 if (qr(i) .lt. 1.e-4) dlt3(i)=1.
3436 dlt4(i)=1.
3438 if (qc(i) .gt. 5.e-4) dlt4(i)=0.0
3440 if (qs(i) .le. 1.e-4) dlt4(i)=1.
3442 if (tair(i) .ge. t0) then
3443 dlt3(i)=0.0
3444 dlt4(i)=0.0
3445 endif
3447 pr(i)=d2t*(qsacw(i)+praut(i)+pracw(i)+wgacr(i) &
3448 +qhacw(i)-qgacr(i))
3449 ps(i)=d2t*(psaut(i)+psaci(i)+dlt4(i)*psacw(i) &
3450 +psfw(i)+psfi(i))
3451 pg(i)=d2t*(dlt3(i)*praci(i)+dgaci(i) &
3452 +wgaci(i)+dgacw(i)+(1.-dlt4(i))*psacw(i) &
3453 +primh(i)) !4ice revised
3456 !* 7 * PRACS : ACCRETION OF QS BY QR ***7**
3457 !* 8 * PSACR : ACCRETION OF QR BY QS (QSACR FOR PSMLT) ***8**
3458 !* 2 * PGAUT : AUTOCONVERSION OF QS TO QG ***2**
3459 !* 18 * PGFR : FREEZING OF QR TO QG **18**
3461 qracs(i)=0.0
3462 y1(i)=abs(vr(i)-vs(i))
3463 y2(i)=zr(i)*zs(i)
3464 y3(i)=5./y2(i)
3465 y4(i)=.08*y3(i)*y3(i)
3466 y5(i)=.05*y3(i)*y4(i)
3468 pracs(i)=r2ig*r2is*(r7rf(i)*y1(i)*(y3(i)/zs(i)**5 & !improve4
3469 +y4(i)/zs(i)**3+y5(i)/zs(i))*ftns(i)*fros(i)) !improve4
3470 if(qs(i).le.cmin) pracs(i)=0.
3471 if(qr(i).le.cmin) pracs(i)=0.
3472 qracs(i)=r2ig*r2is*min(d2t*pracs(i), qs(i))
3473 psacr(i)=r2is*(r8rf(i)*y1(i)*(y3(i)/zr(i)**5 &
3474 +y4(i)/zr(i)**3+y5(i)/zr(i))*ftns(i))
3475 if(qs(i).le.cmin) psacr(i)=0.
3476 if(qr(i).le.cmin) psacr(i)=0.
3477 qsacr(i)=psacr(i)
3479 if (tair(i) .ge. t0) then
3480 pracs(i)=0.0
3481 psacr(i)=0.0
3482 else
3483 qsacr(i)=0.0
3484 qracs(i)=0.0
3485 endif
3487 !* 18 * pgfr : freezing of qr to qg **18**
3488 pgaut(i)=0.0
3489 pgfr(i)=0.0
3490 phfr(i)=0.0
3491 if (tair(i) .lt. t0) then
3492 y1(i)=exp(rn18a*(t0-tair(i)))
3493 if( qr(i).ge.0.)then
3494 temp = 1./zr(i)
3495 temp = temp*temp*temp*temp*temp*temp*temp
3496 phfr(i)=r2ih*max(r18r(i)*(y1(i)-1.)*temp, 0.0)
3497 ! phfr(i)=r2ih*max(r18r(i)*(y1(i)-1.)/zr(i)**7., 0.0)
3498 else
3499 temp = 1./zr(i)
3500 temp = temp*temp*temp*temp*temp*temp*temp
3501 pgfr(i)=r2ig*max(r18r(i)*(y1(i)-1.)*temp, 0.0)
3502 ! pgfr(i)=r2ig*max(r18r(i)*(y1(i)-1.)/zr(i)**7., 0.0)
3503 if(qr(i).le.cmin) pgfr(i)=0.
3504 endif
3505 endif
3507 !!
3508 ! endif ! for Processes 2, 7, 8, & 18
3510 !******** HANDLING THE NEGATIVE RAIN WATER (QR) *******************
3511 !******** HANDLING THE NEGATIVE SNOW (QS) *******************
3513 y1(i)=qr(i)/d2t
3514 y2(i)=-qh(i)/d2t
3515 piacr(i)=min(y1(i), piacr(i))
3516 dgacr(i)=min(y1(i), dgacr(i))
3517 dhacr(i)=min(y1(i), dhacr(i))
3518 ! whacr(i)=min(y1(i), whacr(i))
3519 ! whacr(i)=max(y2(i), whacr(i))
3520 psacr(i)=min(y1(i), psacr(i))
3521 pgfr(i)= min(y1(i), pgfr(i))
3522 phfr(i)= min(y1(i), phfr(i))
3523 del=0.
3524 IF(whacr(i) .LT. 0.) DEL=1.
3525 if(del.eq.1) whacr(i)=max(whacr(i),-dhacw(i)) !fix from JDC
3526 dhacr(i)=min((1.-del)*whacr(i),dhacr(i))
3527 if(del.eq.1) dhacw(i)=dhacw(i)+del*whacr(i) !fix from JDC
3528 y1(i)=(piacr(i)+dgacr(i)+dhacr(i)+psacr(i)+pgfr(i) &
3529 +phfr(i))*d2t
3530 qr(i)=qr(i)+pr(i)+qracs(i)-del*whacr(i)*d2t+qracg(i) !fix SEL
3531 qr(i)=qr(i)-y1(i) !fix SEL
3532 if (qr(i) .lt. 0.0) then
3533 y2(i)=1.
3534 if (y1(i) .ne. 0.0) y2(i)=qr(i)/y1(i)+1.
3535 piacr(i)=piacr(i)*y2(i)
3536 dgacr(i)=dgacr(i)*y2(i)
3537 dhacr(i)=dhacr(i)*y2(i)
3538 ! if(whacr(i).gt.0.) whacr(i)=whacr(i)*y2(i)
3539 pgfr(i)=pgfr(i)*y2(i)
3540 phfr(i)=phfr(i)*y2(i)
3541 psacr(i)=psacr(i)*y2(i)
3542 qr(i)=0.0
3543 endif !qr
3544 dlt2(i)=1.
3545 if (qr(i) .gt. 1.e-4) dlt2(i)=0.
3546 if (tair(i) .ge. t0) dlt2(i)=0.
3547 y1(i)=qs(i)/d2t
3548 pgacs(i)=min(y1(i), pgacs(i))
3549 dgacs(i)=min(y1(i), dgacs(i))
3550 wgacs(i)=min(y1(i), wgacs(i))
3551 whacs(i)=min(y1(i), whacs(i))
3552 pracs(i)=min(y1(i), pracs(i))
3553 pwacs(i)=min(y1(i), pwacs(i))
3555 prn(i)=d2t*(dlt3(i)*piacr(i)+dlt2(i)*dgacr(i)+pgfr(i) &
3556 +dlt2(i)*psacr(i))
3557 pracs(i)=(1.-dlt2(i))*pracs(i)
3558 pwacs(i)=(1.-dlt4(i))*pwacs(i)
3559 pracg(i)=(1.-dlt2(i))*pracg(i)
3560 psn(i)=d2t*(pgacs(i)+dgacs(i)+wgacs(i) &
3561 +pracs(i)+pwacs(i)+whacs(i))
3563 qs(i)=qs(i)+ps(i)-qracs(i)-psn(i)
3564 if (qs(i) .lt. 0.0) then
3565 y2(i)=1.
3566 if (psn(i) .ne. 0.) y2(i)=qs(i)/psn(i)+1.
3567 pgacs(i)=pgacs(i)*y2(i)
3568 dgacs(i)=dgacs(i)*y2(i)
3569 wgacs(i)=wgacs(i)*y2(i)
3570 whacs(i)=whacs(i)*y2(i)
3571 pracs(i)=pracs(i)*y2(i)
3572 pwacs(i)=pwacs(i)*y2(i)
3573 qs(i)=0.0
3574 endif
3575 psn(i)=d2t*(pgacs(i)+dgacs(i)+wgacs(i) &
3576 +pracs(i)+pwacs(i))
3577 qg(i)=qg(i)+pg(i)+prn(i)+psn(i)-qracg(i)
3578 qg(i)=qg(i)-d2t*(pracs(i)+whacg(i) &
3579 +pracg(i) & !fix from JDC
3580 +pg2h(i)) !fix from SEL
3582 if (qg(i) .lt. 0.0) then !fix from JDC
3583 y2(i)=1. !fix from JDC
3584 if (whacg(i)+pracg(i)+pg2h(i).ne. 0.) & !fix from JDC
3585 y2(i)=qg(i)/(d2t*(whacg(i)+pracg(i)+pg2h(i)))+1. !fix from JDC & SEL
3586 whacg(i)=whacg(i)*y2(i) !fix from JDC
3587 pracg(i)=pracg(i)*y2(i) !fix from JDC
3588 pg2h(i)=pg2h(i)*y2(i) !fix from SEL
3589 qg(i)=0.0 !fix from JDC
3590 endif
3592 qh(i)=qh(i)+d2t*(phfr(i)+(1.-dlt3(i))*piacr(i) &
3593 +(1.-dlt3(i))*praci(i)+(1.-dlt2(i))*psacr(i) &
3594 +pracs(i)+(1.-dlt2(i))*dgacr(i)+pracg(i) &
3595 +dhacw(i)+dhacr(i)-primh(i)+whaci(i) &
3596 +whacs(i)+whacg(i)+pg2h(i))
3598 y1(i)=d2t*(psacw(i)+psfw(i)+dgacw(i)+piacr(i) &
3599 +dgacr(i)+psacr(i)+pgfr(i)+phfr(i)+pihms(i) &
3600 +pihmg(i)+dhacw(i)+dhacr(i) +pihmh(i)) &
3601 -qracs(i)
3602 pt(i)=pt(i)+afcp(i)*y1(i)
3603 tair(i)=(pt(i)+tb0)*pi0(i)
3605 !* 11 * PSMLT : MELTING OF QS **11**
3606 !* 19 * PGMLT : MELTING OF QG TO QR **19**
3608 psmlt(i)=0.0
3609 pgmlt(i)=0.0
3610 phmlt(i)=0.0
3612 qhz2=qhwrf(i,k)
3613 qgz2=qgwrf(i,k)
3614 if (k .lt. kte-2 .and. tairc(i) .ge. -5) then
3615 qhz2=qhwrf(i,k+1)
3616 qgz2=qgwrf(i,k+1)
3617 endif
3618 call sgmap(1,qs(i),qg(i),qgz2,qh(i),qhz2,r00(i),tairc(i),ftns0(i))
3619 call sgmap(2,qs(i),qg(i),qgz2,qh(i),qhz2,r00(i),tairc(i),ftng0(i))
3620 call sgmap(3,qs(i),qg(i),qgz2,qh(i),qhz2,r00(i),tairc(i),fros0(i))
3621 call sgmap(4,qs(i),qg(i),qgz2,qh(i),qhz2,r00(i),tairc(i),ftnh0(i))
3623 if (tair(i).ge.t0) then
3624 tairc(i)=tair(i)-t0
3625 ftns(i)=1.
3626 ftng(i)=1.
3627 ftnh(i)=1.
3628 ftns(i)=ftns0(i)
3629 ftng(i)=ftng0(i)
3630 ftnh(i)=ftnh0(i)
3631 dd(i)=r11t*tairc(i)*(r101r(i)/zs(i)**2+r102rf(i) &
3632 /zs(i)**bsh5)*ftns(i)
3633 psmlt(i)=r2is*min(qs(i),max(dd(i),0.0))
3634 if (r00(i)*qg(i).gt.qrog2) then
3635 y2(i)=(r191r(i)/zg(i)**2+r192rf2(i)/zg(i)**bgh5_2)*ftng(i)
3636 else
3637 y2(i)=(r191r(i)/zg(i)**2+r192rf(i)/zg(i)**bgh5)*ftng(i)
3638 endif
3639 dd1(i)=tairc(i)*(r19t*y2(i)+r19at*(qgacw(i) &
3640 +qgacr(i)))
3641 pgmlt(i)=r2ig*min(qg(i),max(dd1(i),0.0))
3643 Y1(i)=TCA(i)*TAIRC(i)-ALVR(i)*DWV(i)*(RP0(i)-(QV(i)+QB0))
3644 Y3(i)=.78/ZH(i)**2+h19aq(i)*SCV(i)/ZH(i)**BHH5
3645 DDa0(i)=h19rt(i)*Y1(i)*Y3(i)*ftnh(i)+r19at*TAIRC(i)* &
3646 (QHACW(i)+QHACR(i))
3647 PHMLT(i)=r2ih*max(0.0, min(DDa0(i), QH(i)))
3649 pt(i)=pt(i)-afcp(i)*(psmlt(i)+pgmlt(i)+phmlt(i))
3650 tair(i)=(pt(i)+tb0)*pi0(i)
3651 qr(i)=qr(i)+psmlt(i)+pgmlt(i)+phmlt(i)
3652 qs(i)=qs(i)-psmlt(i)
3653 qg(i)=qg(i)-pgmlt(i)
3654 qh(i)=qh(i)-phmlt(i)
3656 endif !tair ! processes 11 & 19
3657 !
3658 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3659 !* 24 * PIHOM : HOMOGENEOUS FREEZING OF QC TO QI (T < T00) **24**
3660 !* 25 * PIDW : DEPOSITION GROWTH OF QC TO QI ( T0 < T <= T00) **25**
3661 !* 26 * PIMLT : MELTING OF QI TO QC (T >= T0) **26**
3662 !****** PIMM : IMMERSION FREEZING OF QC TO QI (T < T0) ******
3663 !****** PCFR : CONTACT NUCLEATION OF QC TO QI (T < T0) ******
3665 if (qc(i).le.cmin1) qc(i)=0.0
3666 if (qi(i).le.cmin1) qi(i)=0.0
3668 ftns(i)=1.
3669 ftng(i)=1.
3670 ftns0(i)=1.
3671 ftng0(i)=1.
3673 qhz2=qhwrf(i,k)
3674 qgz2=qgwrf(i,k)
3675 if (k .lt. kte-2 .and. tairc(i) .ge. -5) then
3676 qhz2=qhwrf(i,k+1)
3677 qgz2=qgwrf(i,k+1)
3678 endif
3679 call sgmap(1,qs(i),qg(i),qgz2,qh(i),qhz2,r00(i),tairc(i),ftns0(i))
3680 call sgmap(2,qs(i),qg(i),qgz2,qh(i),qhz2,r00(i),tairc(i),ftng0(i))
3682 if (tair(i).le.t00) then
3683 pihom(i)=qc(i)
3684 else
3685 pihom(i)=0.0
3686 endif
3687 if (tair(i).ge.t0) then
3688 pimlt(i)=qi(i)
3689 else
3690 pimlt(i)=0.0
3691 endif
3692 pidw(i)=0.0
3694 if (tair(i) .lt. t0 .and. tair(i) .gt. t00) then
3695 TAIRC(i)=TAIR(i)-T0
3696 y1(i) = max( min(tairc(i), -1.), -31.)
3697 it(i) = int(abs(y1(i)))
3698 y2(i)=aa1(it(i))
3699 y3(i) = aa2(it(i))
3700 if (tairc(i).le.-5.)then ! meyers
3701 rtair(i)=1./(tair(i)-c76)
3702 y5(i)=exp(c218-c580*rtair(i))
3703 qsi(i)=rp0(i)*y5(i)
3704 SSI(i)=(qv(i)+qb0)/qsi(i)-1.
3705 fssi=min(xssi,max(rssi,xssi*(tairc(i)+44.)/(44.0-38.0))) !max ssi f(tair)
3706 fssi=rssi
3707 wssi=ww1(i,k)/100.-2.0
3708 if(wssi.gt.0.) fssi=rssi+min(xssi,wssi*0.01)
3709 fssi=min(ssi(i),fssi)
3711 r_nci=max(1.e-3*exp(-.639+12.96*fssi), 0.528e-3) ! meyers et al.
3712 ! Roger found the bug on 2012/04/27
3713 if (r_nci.gt.15.) r_nci=15. !cap at 15000/liter
3715 ! Cooper curve
3716 if( tairc(i) .lt. -40.0 ) then
3717 c_nci = 5.0e-6*exp(0.304*40.0)
3718 else
3719 c_nci = 5.0e-6*exp(0.304*abs(tairc(i)))
3720 end if
3721 r_nci = c_nci
3722 ! JDC cap r_nci to the amount corresponding to crystal size of ami50
3723 r_nci = max(r_nci,r00(i)*qi(i)/ami50)
3725 #if (WRF_CHEM == 1)
3726 ! EMK...Only execute if GOCART and coupling is on
3727 if ( (chem_opt == 300 .or. chem_opt == 301 .or. &
3728 chem_opt == 302 .or. chem_opt == 303) .and. &
3729 (gsfcgce_gocart_coupling == 1) ) then
3730 !JJS 20110602 vvvvv
3731 ! Conversion rate of cloud water to ice in the Bergeron porcess based on Meyer + DeMott formulae
3732 !
3733 ! convert gocart aerosol mass conc to IN
3734 ! p0 need to be converted from g*cm/s2/cm2 to mb (hPa)
3735 ! call mass2icn(p0(i,k)*0.001,tair(i),aero(i,k,:), icn_out)
3736 ! icn_cgs(i,k) = max(min_icn, icn_out) * 1.e-3 !IN conc [#/cm3] <-- [#/Litter]
3737 ! icn_cgs(i,k) = icn_out * 1.e-3 !IN conc [#/cm3]
3738 ! r_nci = icn_cgs(i,k) !DeMotto's formuale
3739 r_nci = icn_out(i) * 1.e-3 !DeMotto's formuale
3740 !JJS 20110602 ^^^^^
3741 end if
3742 #endif
3744 dd(i)=(r00(i)*qi(i)/r_nci)**y3(i) !meyers
3745 PIDW(i)=min(rr0(i)*D2T*y2(i)*r_nci*dd(i), qc(i)) !meyers
3747 ! JDC water vapor diffusivity correction term
3748 esi(i) = qsi(i)/rp0(i)*c610
3749 y4(i) = 1./tair(i)
3750 ! term1 = y4(i)*(rn10a*y4(i)-rn10b)
3751 term1 = y4(i)*(rn20a*y4(i)-rn20b)
3752 term2 = rn10c*tair(i)/esi(i)
3753 fdwv = (term1+term2)/(term1+term2*dwv0/dwv(i))
3754 PIDW(i)=min(rr0(i)*D2T*y2(i)*r_nci*dd(i)*fdwv,qc(i)) !meyers
3755 endif !tairc
3757 pimm(i)=0.0
3758 pcfr(i)=0.0
3760 if (qc(i) .gt. 0.0) then
3761 y4(i) = 1./(tair(i)-c358)
3762 qsw(i)=rp0(i)*exp(c172-c409*y4(i))
3763 xncld=qc(i)/4.e-9 !cloud number
3764 esat=0.6112*exp(17.67*tairc(i)/(tairc(i)+243.5))*10.
3765 rv=0.622*esat/(p0(i,k)/1000.-esat)
3766 rlapse_m=980.616*(1.+2.5e6*rv/287./tair(i))/ &
3767 (1004.67+2.5e6*2.5e6*rv*0.622/(287.*tair(i)*tair(i)))
3768 delT=rlapse_m*ww1(i,k) !Roger
3769 if (delT.lt.0.) delT=0.
3770 Bhi=1.01e-2
3771 pimm(i)=xncld*Bhi*4.e-9*exp(-tairc(i))*delT*d2t*4.e-9
3773 xccld=xncld*r00(i) !cloud number concentration
3774 Xknud=7.37*tair(i)/(288.*Ra*p0(i,k)) ! Knudsen number
3775 alpha=1.257+0.400*exp(-1.10/Xknud) ! Cunningham correction (P&Klett)
3777 cunnF=1.+alpha*Xknud ! Cunningham correction (P&Klett)
3779 if (tairc(i).ge.0.) then
3780 dvair=(1.718+0.0049*tairc(i))*1.e-4
3781 !dynamic visc air (Prupp&Klett)
3782 else
3783 dvair=(1.718+0.0049*tairc(i)-1.2e-5*tairc(i)**2)*1.e-4
3784 !dynamic visc air (Prupp&Klett)
3785 endif !tairc
3786 DIFFar=1.3804e-16*tair(i)/6./cpi/dvair/Ra*cunnF !aerosol diffusion via P&Klett
3787 if (qv(i)+qb0-qsw(i).lt.0.) then !only when cloud evaporating
3788 pcfr(i)=4.e-9*4.*cpi*Rc*DIFFar*xccld*Cna*rr0(i)*d2t !Brownian part only via Cotton
3789 endif
3790 else !qc
3791 tairc(i)=tair(i)-t0
3792 y1(i)=max( min(tairc(i), -1.), -31.)
3793 it(i)=int(abs(y1(i)))
3794 y2(i)=aa1(it(i))
3795 y3(i)=aa2(it(i))
3796 y4(i)=exp(abs(.5*tairc(i)))
3797 dd(i)=(r00(i)*qi(i)/(r25a*y4(i)))**y3(i)
3798 pidw(i)=min(r25rt(i)*y2(i)*y4(i)*dd(i),qc(i))
3799 endif !qc
3800 endif !tair
3802 ENDDO
3803 !dir$ vector aligned
3804 DO i=1,irestrict
3806 ! STEVE: PLEASE CHECK
3807 y1(i)=pihom(i)+pidw(i)+pimm(i)+pcfr(i)-pimlt(i)
3809 if (y1(i).gt.qc(i)) then
3810 y1(i)=qc(i)
3811 y2(i)=1.
3812 y3(i)=pihom(i)+pidw(i)+pimm(i)+pcfr(i)
3813 if(y3(i).ne.0.) y2(i)=(qc(i)+pimlt(i))/y3(i)
3814 pihom(i)=pihom(i)*y2(i)
3815 pidw(i)=pidw(i)*y2(i)
3816 pimm(i)=pimm(i)*y2(i)
3817 pcfr(i)=pcfr(i)*y2(i)
3818 endif !y1
3820 pt(i)=pt(i)+afcp(i)*y1(i)
3821 tair(i)=(pt(i)+tb0)*pi0(i)
3822 qc(i)=qc(i)-y1(i)
3823 qi(i)=qi(i)+y1(i)
3825 !* 31 * pint : initiation of qi **31**
3826 !* 32 * pidep : deposition of qi **32**
3827 !
3828 ! CALCULATION OF PINT USES DIFFERENT VALUES OF THE INTERCEPT AND SLOPE FOR
3829 ! THE FLETCHER EQUATION. ALSO, ONLY INITIATE MORE ICE IF THE NEW NUMBER
3830 ! CONCENTRATION EXCEEDS THAT ALREADY PRESENT.
3832 tair(i)=(pt(i)+tb0)*pi0(i)
3833 if (tair(i) .lt. t0) then
3834 if (qi(i) .le. cmin) qi(i)=0.
3835 tairc(i)=tair(i)-t0
3836 rtair(i)=1./(tair(i)-c76)
3837 y2(i)=exp(c218-c580*rtair(i))
3838 qsi(i)=rp0(i)*y2(i)
3839 esi(i)=c610*y2(i)
3840 ssi(i)=(qv(i)+qb0)/qsi(i)-1.
3841 y1(i)=1./tair(i)
3842 y3(i)=SQRT(qi(i))
3843 ! JDC water vapor diffusivity correction term
3844 term1 = y1(i)*(rn20a*y1(i)-rn20b)
3845 term2 = rn10c*tair(i)/esi(i)
3846 dd(i) = term1+term2
3847 fdwv = dd(i)/(term1+term2*dwv0/dwv(i))
3849 dm(i)=max(qv(i)+qb0-qsi(i),0.0)
3850 rsub1(i)=cs580(i)*qsi(i)*rtair(i)*rtair(i)
3851 dep(i)=dm(i)/(1.+rsub1(i))
3852 if (tairc(i).le.-5.) then
3853 y4(i)=1./(tair(i)-c358)
3854 qsw(i)=rp0(i)*exp(c172-c409*y4(i))
3855 fssi=min(xssi,max(rssi,xssi*(tairc(i)+44.)/(44.-38.)))
3856 fssi=rssi
3857 wssi=ww1(i,k)/100.-2.0
3858 if(wssi.gt.0.) fssi=rssi+min(xssi,wssi*0.01)
3859 fssi=min(ssi(i),fssi)
3860 ! r_nci=min(1.e-3*exp(-.639+12.96*fssi),1.)
3861 r_nci=max(1.e-3*exp(-.639+12.96*fssi),0.528e-3)
3862 if (r_nci.gt.15.) r_nci=15.
3864 ! Cooper curve
3865 if( tairc(i) .lt. -40.0 ) then
3866 c_nci = 5.0e-6*exp(0.304*40.0)
3867 else
3868 c_nci = 5.0e-6*exp(0.304*abs(tairc(i)))
3869 end if
3870 r_nci = c_nci
3871 ! JDC cap r_nci to the amount corresponding to crystal size of ami50
3872 r_nci = max(r_nci,r00(i)*qi(i)/ami50)
3874 #if (WRF_CHEM == 1)
3875 ! EMK...Only execute if GOCART and coupling is on
3876 if ( (chem_opt == 300 .or. chem_opt == 301 .or. &
3877 chem_opt == 302 .or. chem_opt == 303) .and. &
3878 (gsfcgce_gocart_coupling == 1) ) then
3879 !JJS 20110602 vvvvv
3880 ! Conversion rate of cloud water to ice in the Bergeron porcess based on Meyer + DeMott formulae
3881 !
3882 ! convert gocart aerosol mass conc to IN
3883 ! p0 need to be converted from g*cm/s2/cm2 to mb (hPa)
3884 ! call mass2icn(p0(i,k)*0.001,tair(i),aero(i,k,:), icn_out)
3885 ! icn_cgs(i,k) = max(min_icn, icn_out) * 1.e-3 !IN conc [#/cm3] <-- [#/Litter]
3886 ! call mass2icn(p0(i,k)*0.001,tair(i),aero(i,k,:), icn_out,i,k)
3887 ! icn_out = min(1.e3, max(0.01e0 , icn_out) )
3888 ! icn_cgs(i,k) = icn_out * 1.e-3 !IN conc [#/cm3]
3889 ! r_nci = icn_cgs(i,k) !DeMotto's formuale
3890 r_nci = icn_out(i) * 1.e-3 !DeMotto's formuale
3891 !JJS 20110602 ^^^^^
3892 end if
3893 #endif
3895 pidep(i)=max(R32RT(i)*1.e4*fssi*sqrt(r_nci)*y3(i)/ & !meyers
3896 dd(i)*fdwv, -qi(i)) !fix SEL
3897 if(qi(i).le.cmin) pidep(i)=0.
3898 dd(i)=max(1.e-9*r_nci/r00(i)-qci(i,k)*1.e-9/ami50, 0.)
3899 pint(i)=max(min(dd(i),dm(i)),0.)
3900 pint(i)=min(pint(i)+pidep(i), dep(i))
3901 pt(i)=pt(i)+ascp(i)*pint(i)
3902 tair(i)=(pt(i)+tb0)*pi0(i)
3903 qv(i)=qv(i)-pint(i)
3904 qi(i)=qi(i)+pint(i)
3905 endif !taric
3906 endif !tair
3908 ! End of Process 31 & 32
3910 ! WRF satice has new_ice_sat option 0, 1 and 2
3911 ! Steve's satice has new_ice_sat option 0, 1, 2, 3 and 9
3912 ! option 0, 1 and 2 are identical in both satice
3913 ! I added option 3 below and wrapped them with "if (improve.eq.3)"
3915 !***** TAO ET AL (1989) SATURATION TECHNIQUE ***********************
3917 !!!!
3918 !!! new_ice_sat option 9 from Steve's satice
3919 !!! sequential, non-iterative, ssi
3920 !!!!
3922 dep(i)=0.0
3923 cnd(i)=0.0
3924 tair(i)=(pt(i)+tb0)*pi0(i)
3925 if (tair(i).ge.t00) THEN
3926 y1(i)=1./(tair(i)-c358)
3927 qsw(i)=rp0(i)*exp(c172-c409*y1(i))
3928 dd(i)=cp409(i)*y1(i)*y1(i)
3929 dm(i)=qv(i)+qb0-qsw(i)
3930 cnd(i)=dm(i)/(1.+avcp(i)*dd(i)*qsw(i))
3931 !c ****** condensation or evaporation of qc ******
3932 cnd(i)=max(-qc(i),cnd(i))
3934 ! JJS 20150831
3935 ! if(ww1(i,k).ge.-10.) & ! for 1-km grid
3936 ! thresh_evap is a function of dx and defined in the beginning of the subroutine
3938 if(ww1(i,k) .ge. thresh_evap) cnd(i)=max(0.,cnd(i)) !reduce spurious evap
3940 pt(i)=pt(i)+avcp(i)*cnd(i)
3941 tair(i)=(pt(i)+tb0)*pi0(i)
3942 qv(i)=qv(i)-cnd(i)
3943 qc(i)=qc(i)+cnd(i)
3944 endif
3945 if (tair(i).le.273.16) THEN
3946 !c ****** deposition or sublimation of qi ******
3947 y1(i)=1./(tair(i)-c358)
3948 qsw(i)=rp0(i)*exp(c172-c409*y1(i))
3949 y2(i)=1./(tair(i)-c76)
3950 qsi(i)=rp0(i)*exp(c218-c580*y2(i))
3951 fssi=min(xssi,max(rssi,xssi*(tair(i)-t0+44.0)/(44.0-38.0)))
3952 fssi=rssi
3953 wssi=ww1(i,k)/100.-2.0
3954 if(wssi.gt.0.) fssi=rssi+min(xssi,wssi*0.01)
3955 y3(i)=1.+min((qsw(i)-qsi(i))/qsi(i), fssi)
3956 y4(i)=qsi(i)*y3(i)
3957 if (tair(i).le.268.16.and.(qv(i)+qb0.gt.y4(i))) then
3958 dd1(i)=cp580(i)*y2(i)*y2(i)
3959 dep(i)=(qv(i)+qb0-y4(i))/(1.+ascp(i)*dd1(i)*y4(i))
3960 else if (qv(i)+qb0.lt.xsubi*qsi(i).and.qi(i).gt.cmin) then
3961 dd1(i)=cp580(i)*y2(i)*y2(i)
3962 dep(i)=(qv(i)+qb0-xsubi*qsi(i))/(1.+ascp(i)*dd1(i)*xsubi*qsi(i))
3963 dep(i)=max(-qi(i),dep(i))
3964 endif
3965 if(ww1(i,k).ge.0.) &
3966 dep(i)=max(0.,dep(i)) !reduce spurious sublimation
3967 pt(i)=pt(i)+ascp(i)*dep(i)
3968 tair(i)=(pt(i)+tb0)*pi0(i)
3969 qv(i)=qv(i)-dep(i)
3970 qi(i)=qi(i)+dep(i)
3971 endif
3973 !* 10 * PSDEP : DEPOSITION OR SUBLIMATION OF QS **10**
3974 !* 20 * PGSUB : SUBLIMATION OF QG **20**
3976 psdep(i)=0.0
3977 pgdep(i)=0.0
3978 phdep(i)=0.0
3979 pssub(i)=0.0
3980 pgsub(i)=0.0
3981 phsub(i)=0.0
3982 pvapg(i)=0.0
3983 pvaph(i)=0.0
3985 if (qc(i)+qi(i).gt.1.e-5) then
3986 dlt1(i)=1.
3987 else
3988 dlt1(i)=0.
3989 endif
3991 if (tair(i) .lt. t0) then
3992 rtair(i)=1./(tair(i)-c76)
3993 y2(i)=exp(c218-c580*rtair(i))
3994 qsi(i)=rp0(i)*y2(i)
3995 esi(i)=c610*y2(i)
3996 ftns(i)=1.
3997 ftng(i)=1.
3998 ftnh(i)=1.
4000 SSI(i)=(QV(i)+QB0)/QSI(i)-1.
4001 IF(SSI(i).GT.0.) DLT1(i)=1.
4002 IF(SSI(i).LE.0.) DLT1(i)=0.
4003 DM(i)=QV(i)+QB0-QSI(i)
4004 RSUB1(i)=cs580(i)*QSI(i)*RTAIR(i)*RTAIR(i)
4005 DD1(i)=DM(i)/(1.+RSUB1(i))
4006 Y3(i)=1./TAIR(i)
4007 ! JDC water vapor diffusivity correction
4008 term1 = y3(i)*(rn20a*y3(i)-rn20b)
4009 term2 = rn10c*tair(i)/esi(i)
4010 dd(i) = term1+term2
4011 fdwv = dd(i)/(term1+term2*dwv0/dwv(i))
4012 TAIRC(i)=TAIR(i)-T0
4013 !
4014 ftns(i)=1.
4015 ftng(i)=1.
4016 ftnh(i)=1.
4017 ftns0(i)=1.
4018 ftng0(i)=1.
4019 qhz2=qhwrf(i,k)
4020 qgz2=qgwrf(i,k)
4021 if (k .lt. kte-2 .and. tairc(i) .ge. -5) then
4022 qhz2=qhwrf(i,k+1)
4023 qgz2=qgwrf(i,k+1)
4024 endif
4025 call sgmap(1,qs(i),qg(i),qgz2,qh(i),qhz2,r00(i),tairc(i),ftns0(i))
4026 call sgmap(2,qs(i),qg(i),qgz2,qh(i),qhz2,r00(i),tairc(i),ftng0(i))
4027 call sgmap(3,qs(i),qg(i),qgz2,qh(i),qhz2,r00(i),tairc(i),fros0(i))
4028 call sgmap(4,qs(i),qg(i),qgz2,qh(i),qhz2,r00(i),tairc(i),ftnh0(i))
4030 ftns(i)=ftns0(i)
4031 ftng(i)=ftng0(i)
4032 ftnh(i)=ftnh0(i)
4034 Y4(i)=r10t*SSI(i)*(r101r(i)/ZS(i)**2+r102rf(i)/ZS(i)**BSH5) &
4035 /DD(i)*ftns(i)*fdwv
4036 PSDEP(i)=r2is*max(-QS(i), Y4(i))
4037 if(qs(i).le.cmin) psdep(i)=0.
4038 DD(i)=Y3(i)*(RN20A*Y3(i)-RN20B)+RN10C*TAIR(i)/ESI(i)
4039 Y2(i)=r191r(i)/ZG(i)**2+r192rf(i)/ZG(i)**BGH5
4040 if(r00(i)*qg(i).gt.qrog2) &
4041 y2(i)=(r191r(i)/zg(i)**2+r192rf2(i)/zg(i)**bgh5_2)
4042 PGDEP(i)=r2ig*MAX(-qg(i), R20T*SSI(i)*Y2(i)/DD(i) &
4043 *ftng(i)*fdwv)
4044 if(qg(i).le.cmin) pgdep(i)=0.
4045 Y1(i)=h10ar(i)/(TCA(i)*TAIR(i)**2)+1./(DWV(i)*QSI(i))
4046 Y2(i)=.78/ZH(i)**2+h20bq(i)*SCV(i)/ZH(i)**BHH5
4047 PHDEP(i)=r2ih*MAX(-qh(i), h20t(i)*ftnh(i)*SSI(i)*Y2(i)/Y1(i))
4048 if(qh(i).le.cmin) phdep(i)=0.
4049 PSSUB(i)=min(PSDEP(i), 0.)
4050 PSDEP(i)=max(PSDEP(i), 0.)
4051 PGSUB(i)=min(PGDEP(i), 0.)
4052 PGDEP(i)=max(PGDEP(i), 0.)
4053 PHSUB(i)=min(PHDEP(i), 0.)
4054 PHDEP(i)=max(PHDEP(i), 0.)
4056 ! ******************************************************************
4057 Y5(i)=min(0.,DD1(i))
4058 DD1(i)=max(0.,DD1(i))
4060 IF(DLT1(i).EQ.1.)THEN !Di
4061 Y1(i)=PSDEP(i)+PGDEP(i)+PHDEP(i)
4062 IF(Y1(i).ge.DD1(i))THEN
4063 PSDEP(i)=PSDEP(i)/Y1(i)*DD1(i) !...
4064 PGDEP(i)=PGDEP(i)/Y1(i)*DD1(i)
4065 PHDEP(i)=PHDEP(i)/Y1(i)*DD1(i)
4066 ENDIF
4067 ENDIF !Di
4069 if(qc(i).le.1.e-5) then
4070 vap_frac=2.0*min((tairc(i)/(t00-t0))**2,1.0)
4071 pvapg(i)=vap_frac*pgdep(i)
4072 pvaph(i)=vap_frac*phdep(i)
4073 endif
4074 if(pgdep(i).gt.0.) &
4075 pvapg(i)=min(pvapg(i),qg(i)+pgdep(i))
4076 if(phdep(i).gt.0.) &
4077 pvaph(i)=min(pvaph(i),qh(i)+phdep(i))
4080 IF(DLT1(i).EQ.0.)THEN
4081 Y1(i)=MAX(PSsub(i)+PGsub(i)+phsub(i), Y5(i))
4082 IF(Y5(i).gt.(PSsub(i)+PGsub(i)+phsub(i)))THEN
4083 Y3(i)=(PSsub(i)+PGsub(i)+phsub(i))
4084 IF(Y3(i).ne.0.0)THEN
4085 PSsub(i)=PSsub(i)/Y3(i)*Y5(i)
4086 PGsub(i)=PGsub(i)/Y3(i)*Y5(i)
4087 Phsub(i)=Phsub(i)/Y3(i)*Y5(i)
4088 ENDIF
4089 ENDIF
4090 ENDIF
4091 PSSUB(i)=-PSSUB(i)
4092 PGSUB(i)=-PGSUB(i)
4093 PHSUB(i)=-PHSUB(i)
4094 Y1(i)=PSDEP(i)+PGDEP(i)+PHDEP(i) &
4095 -PSsub(i)-PGsub(i)-PHsub(i)
4097 pt(i)=pt(i)+ascp(i)*y1(i)
4098 tair(i)=(pt(i)+tb0)*pi0(i)
4099 qv(i)=qv(i)-y1(i)
4100 qs(i)=qs(i)+psdep(i)-pssub(i)+pvapg(i)+pvaph(i)
4101 qg(i)=qg(i)+pgdep(i)-pgsub(i)-pvapg(i)
4102 qh(i)=qh(i)+phdep(i)-phsub(i)-pvaph(i)
4103 endif ! if (tair(i) .lt. t0)
4105 !* 23 * ERN : EVAPORATION OF QR (SUBSATURATION) **23**
4106 ! Steve did not make any improvement on this process
4107 ern(i)=0.0
4108 if (qr(i) .gt. 0.0) then
4109 tair(i)=(pt(i)+tb0)*pi0(i)
4110 rtair(i)=1./(tair(i)-c358)
4111 y2(i)=exp( c172-c409*rtair(i) )
4112 esw(i)=c610*y2(i)
4113 qsw(i)=rp0(i)*y2(i)
4114 ssw(i)=(qv(i)+qb0)/qsw(i)-1.
4115 dm(i)=qv(i)+qb0-qsw(i)
4116 rsub1(i)=cv409(i)*qsw(i)*rtair(i)*rtair(i)
4117 dd1(i)=max(-dm(i)/(1.+rsub1(i)),0.0)
4118 y3(i)=1./tair(i)
4119 ! JDC water vapor diffusivity correction term
4120 term1 = y3(i)*(rn30a*y3(i)-rn10b)
4121 term2 = rn10c*tair(i)/esw(i)
4122 dd(i) = term1+term2
4123 fdwv = dd(i)/(term1+term2*dwv0/dwv(i))
4124 ftnw=1.
4125 if (qr(i) .gt. cmin.and.tair(i).gt.t0) then !no need to check qc, no qc if ssw < 0
4126 bin_factor(i)=0.11*(1000.*qr(i))**(-1.27) + 0.98
4127 bin_factor(i)=min(bin_factor(i),1.30)
4128 ftnw=1./bin_factor(i)**3.35
4129 ftnwmin=r00(i)*qr(i)/draimax
4130 if(qr(i).le.0.001) ftnw=max(ftnw,ftnwmin/tnw)
4132 y4(i)=r00(i)*qr(i)
4133 y1(i)=sqrt(y4(i))
4134 y2(i)=sqrt(y1(i))
4135 zr(i)=zrc/y2(i)*ftnw**0.25
4136 endif !qr
4137 y1(i)=-r23t*ssw(i)*(r231r(i)/zr(i)**2+r232rf(i)/ &
4138 zr(i)**3)/dd(i)*ftnw*fdwv
4139 ern(i)=min(dd1(i),qr(i),max(y1(i),0.0))
4140 pt(i)=pt(i)-avcp(i)*ern(i)
4141 tair(i)=(pt(i)+tb0)*pi0(i)
4142 tairc(i)=tair(i)-t0
4143 qv(i)=qv(i)+ern(i)
4144 qr(i)=qr(i)-ern(i)
4145 endif !qr
4147 !!!!
4148 !! add processes 30 & 33 fpr pmltg and pmlts
4149 !!
4150 !!!!
4152 !* 30 * pmltg : evaporation of melting qg **30**
4153 !* 33 * pmlts : evaporation of melting qs **33**
4155 pmlts(i)=0.0
4156 pmltg(i)=0.0
4158 tair(i)=(pt(i)+tb0)*pi0(i)
4159 tairc(i)=tair(i)-t0
4161 ftns0(i)=1.
4162 ftng0(i)=1.
4164 qhz2=qhwrf(i,k)
4165 qgz2=qgwrf(i,k)
4166 if (k .lt. kte-2 .and. tairc(i) .ge. -5) then
4167 qhz2=qhwrf(i,k+1)
4168 qgz2=qgwrf(i,k+1)
4169 endif
4170 call sgmap(1,qs(i),qg(i),qgz2,qh(i),qhz2,r00(i),tairc(i),ftns0(i))
4171 call sgmap(2,qs(i),qg(i),qgz2,qh(i),qhz2,r00(i),tairc(i),ftng0(i))
4172 call sgmap(3,qs(i),qg(i),qgz2,qh(i),qhz2,r00(i),tairc(i),fros0(i))
4174 if (tair(i) .ge. t0) then
4175 ftns(i)=1.
4176 ftng(i)=1.
4177 ftns(i)=ftns0(i)
4178 ftng(i)=ftng0(i)
4179 rtair(i)=1./(t0-c358)
4180 y2(i)=exp( c172-c409*rtair(i) )
4181 esw(i)=c610*y2(i)
4182 qsw(i)=rp0(i)*y2(i)
4183 ssw(i)=1.-(qv(i)+qb0)/qsw(i)
4184 dm(i)=qsw(i)-qv(i)-qb0
4185 rsub1(i)=cv409(i)*qsw(i)*rtair(i)*rtair(i)
4186 dd1(i)=max(dm(i)/(1.+rsub1(i)),0.0)
4187 y3(i)=1./tair(i)
4188 ! JDC water vapor diffusivity correction term
4189 term1 = y3(i)*(rn30a*y3(i)-rn10b)
4190 term2 = rn10c*tair(i)/esw(i)
4191 dd(i) = term1+term2
4192 fdwv = dd(i)/(term1+term2*dwv0/dwv(i))
4193 y1(i)=ftng(i)*r30t*ssw(i)*(r191r(i)/zg(i)**2+r192rf(i) &
4194 /zg(i)**bgh5)/dd(i)*fdwv
4195 if(r00(i)*qg(i).gt.qrog2) &
4196 y1(i)=ftng(i)*r30t*ssw(i)*(r191r(i)/zg(i)**2+r192rf2(i) &
4197 /zg(i)**bgh5_2)/dd(i)*fdwv
4198 pmltg(i)=r2ig*min(qg(i),max(y1(i),0.0))
4199 y1(i)=ftns(i)*r33t*ssw(i)*(r331r(i)/zs(i)**2+r332rf(i) &
4200 /zs(i)**bsh5)/dd(i)*fdwv
4201 pmlts(i)=r2is*min(qs(i),max(y1(i),0.0))
4202 y1(i)=min(pmltg(i)+pmlts(i),dd1(i))
4203 pmltg(i)=y1(i)-pmlts(i)
4204 pt(i)=pt(i)-ascp(i)*y1(i)
4205 tair(i)=(pt(i)+tb0)*pi0(i)
4206 qv(i)=qv(i)+y1(i)
4207 qs(i)=qs(i)-pmlts(i)
4208 qg(i)=qg(i)-pmltg(i)
4209 endif !tair t0
4210 ! end Processes 30 and 33
4212 if (qc(i) .le. cmin) qc(i)=0.
4213 if (qr(i) .le. cmin) qr(i)=0.
4214 if (qi(i) .le. cmin) qi(i)=0.
4215 if (qs(i) .le. cmin) qs(i)=0.
4216 if (qg(i) .le. cmin) qg(i)=0.
4217 if (qh(i) .le. cmin) qh(i)=0.
4218 dpt(i,k)=pt(i)
4219 dqv(i,k)=qv(i)
4220 qcl(i,k)=qc(i)
4221 qrn(i,k)=qr(i)
4222 qci(i,k)=qi(i)
4223 qcs(i,k)=qs(i)
4224 qcg(i,k)=qg(i)
4225 qch(i,k)=qh(i)
4227 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
4228 !c henry: please take a look (start)
4229 !JJS modified by JJS on 5/1/2007 vvvvv
4231 dd(i)=max(-cnd(i), 0.)
4232 cnd(i)=max(cnd(i), 0.)
4233 dd1(i)=max(-dep(i), 0.) !bug fix by Di
4234 dep(i)=max(dep(i), 0.)
4235 del=0.
4236 IF(whacr(i) .LT. 0.) DEL=1.
4239 !!!!!!!!!!!DDDDDDDDDDDDDD double check
4240 sccc=cnd(i)
4241 seee=dd(i)+ern(i)
4242 sddd=dep(i)+amax1(pint(i),0.0)+psdep(i)+pgdep(i)+phdep(i)
4243 ssss=dd1(i)-amin1(pint(i),0.0)+pssub(i)+pgsub(i)+phsub(i)+pmlts(i)+pmltg(i)
4244 smmm=psmlt(i)+pgmlt(i)+pimlt(i)+qracs(i)+phmlt(i)+qracg(i) &
4245 -del*whacr(i)
4246 sfff=psacw(i)*d2t+piacr(i)*d2t+psfw(i)*d2t+pgfr(i)*d2t &
4247 +dgacw(i)*d2t+dgacr(i)*d2t+psacr(i)*d2t+pihom(i) &
4248 +pidw(i)+pimm(i)+pcfr(i)+pihms(i)*d2t &
4249 +pihmg(i)*d2t+phfr(i)*d2t+dhacw(i)*d2t+dhacr(i)*d2t+pihmh(i)*d2t
4251 ! for snapsot diabatic heating rate (deg K / s)
4252 physc(i,k) = avc * sccc / d2t !K/s
4253 physe(i,k) = avc * seee / d2t !K/s
4254 physd(i,k) = asc * sddd / d2t !K/s
4255 physs(i,k) = asc * ssss / d2t !K/s
4256 physf(i,k) = afc * sfff / d2t !K/s
4257 physm(i,k) = afc * smmm / d2t !K/s
4258 ! for accmulated diabatic heating (unit: deg K, in potential temperature)
4259 acphysc(i,k) = acphysc(i,k) + avcp(i) * sccc
4260 acphyse(i,k) = acphyse(i,k) + avcp(i) * seee
4261 acphysd(i,k) = acphysd(i,k) + ascp(i) * sddd
4262 acphyss(i,k) = acphyss(i,k) + ascp(i) * ssss
4263 acphysf(i,k) = acphysf(i,k) + afcp(i) * sfff
4264 acphysm(i,k) = acphysm(i,k) + afcp(i) * smmm
4266 !JJS modified by JJS on 5/1/2007 ^^^^^
4268 !JJS 2010/10/19 vvvvv
4269 ! radar reflectivity calculation
4271 call sgmap(1,qs(i),qg(i),qgz2,qh(i),qhz2,r00(i),tairc(i),ftns0(i)) !Di 20160114
4272 call sgmap(2,qs(i),qg(i),qgz2,qh(i),qhz2,r00(i),tairc(i),ftng0(i)) !Di 20160114
4273 call sgmap(3,qs(i),qg(i),qgz2,qh(i),qhz2,r00(i),tairc(i),fros0(i)) !Di 20160114
4274 call sgmap(4,qs(i),qg(i),qgz2,qh(i),qhz2,r00(i),tairc(i),ftnh0(i)) !Di 20160114
4276 a_1=1.e6*r00(i)*qr(i)
4277 a_2=1.e6*r00(i)*qs(i)
4278 a_3=1.e6*r00(i)*qg(i)
4279 a_4=1.e6*r00(i)*qh(i)
4281 ucor=3071.29/tnw**0.75
4282 ucos=687.97*roqs**0.25/tns**0.75
4283 ucog=687.97*roqg**0.25/tng**0.75
4284 ucog2=687.97*roqg2**.25/tng**.75
4285 ucoh=687.97*roqh**0.25/tnh**0.75
4286 uwet=4.464**0.95
4288 ftnw=1.
4289 if(qr(i).gt.cmin .and. qc(i).lt.cmin)then
4290 bin_factor(i)=0.11*(1000.*qr(i))**(-1.27) + 0.98
4291 ! bin_factor(i)=min(bin_factor(i),1.35)
4292 bin_factor(i)=min(bin_factor(i),1.30)
4293 ftnw=1./bin_factor(i)**3.50
4294 ftnwmin=r00(i)*qr(i)/draimax
4295 if(qr(i).le.0.001) ftnw=max(ftnw,ftnwmin/tnw)
4296 endif
4297 a_11=ucor*(max(0.,a_1))**1.75/ftnw**0.75
4299 ftns(i)=1.
4300 ftng(i)=1.
4301 ftns(i)=ftns0(i)**0.75
4302 ftng(i)=ftng0(i)**0.75
4303 ftnh(i)=ftnh0(i)**0.75 ! Di 20160114
4304 fros(i)=1.
4305 fros(i)=fros0(i)**.25
4306 a_22=ucos*(max(0.,a_2))**1.75/ftns(i)*fros(i)
4307 a_33=ucog*(max(0.,a_3))**1.75/ftng(i)
4308 if(a_3.ge.qrog2) a_33=ucog2*(max(0.,a_3))**1.75/ftng(i)
4309 a_44=ucoh*(max(0.,a_4))**1.75/ftnh(i) ! Di 20160114
4311 IF (TAIR(i).LT.273.16) THEN
4312 ZDRY = MAX(1.e-9,A_11+A_22+A_33+A_44) !rain,snow,graupel,hail,cloud ice,cloud water
4313 DBZ(i,k) = 10.*ALOG10(ZDRY)
4314 ELSE
4315 ZWET0 = A_11+UWET*(A_22+A_33+A_44)**.95
4316 ZWET = MAX(1.e-9,ZWET0)
4317 DBZ(i,k) = 10.*ALOG10(ZWET)
4318 ENDIF
4320 !JJS 2010/10/19 ^^^^^
4322 ENDDO
4323 !dir$ vector aligned
4324 DO i=1,irestrict
4326 #if ( WRF_CHEM == 1)
4327 ! EMK...Nuclei only calculated when coupling
4328 if ( (chem_opt == 300 .or. chem_opt == 301 .or. &
4329 chem_opt == 302 .or. chem_opt == 303) .and. &
4330 (gsfcgce_gocart_coupling == 1) ) then
4331 icn_diag(i,k) = icn_out(i) ! #/Litre
4332 nc_diag(i,k) = ccn_out(i) ! #/cm3
4333 else
4334 icn_diag(i,k) = 0
4335 nc_diag(i,k) = 0
4336 end if
4337 #endif
4339 !JJS 20140305 vvvvv Calculate effective radius for all cloud species
4340 ! eff_rad is a function of the slope parameter (Lambda)
4342 ! rain
4343 if (qrn(i,k) .lt. cmin) then
4344 re_rain_gsfc(i,k) = 0.e0
4345 else
4346 re_rain_gsfc(i,k) = eff_rad(zr(i))
4347 endif
4348 ! snow
4349 if (qcs(i,k) .lt. cmin) then
4350 re_snow_gsfc(i,k) = 0.e0
4351 else
4352 re_snow_gsfc(i,k) = eff_rad(zs(i))
4353 endif
4354 ! graupel
4355 if (qcg(i,k) .lt. cmin) then
4356 re_graupel_gsfc(i,k) = 0.e0
4357 else
4358 re_graupel_gsfc(i,k) = eff_rad(zg(i))
4359 endif
4360 ! hail
4361 if (qch(i,k) .lt. cmin) then
4362 re_hail_gsfc(i,k) = 0.e0
4363 else
4364 re_hail_gsfc(i,k) = eff_rad(zh(i))
4365 endif
4367 ! for cloud water
4369 if (qcl(i,k) .lt. cmin) then
4370 re_cloud_gsfc(i,k) = 0.e0
4371 else
4372 L_cloud = qcl(i,k) * rho(i,k) ! cloud water [g/cm3]
4373 #if (WRF_CHEM == 1)
4374 ! when running with WRF_Chem and using aerosol coupling in Goddard MP
4375 ! cpi: const_pi = 4.*atan(1.) ~ 3.1415
4376 ! roqr: 1.0 g/cm**3, liquid water density
4377 ! roqi: 0.9179, ice density
4378 if ( (chem_opt == 300 .or. chem_opt == 301 .or. &
4379 chem_opt == 302 .or. chem_opt == 303) .and. &
4380 (gsfcgce_gocart_coupling == 1) ) then
4381 ! for cloud water, estimate lambda (slope of gamma distribution)
4382 mu = min(15.e0, (1000.E0/ccn_out(i) + 2.e0))
4383 gamfac3 = ( gamma_toshi(mu+4.e0) / gamma_toshi(mu+3.e0) )
4384 gamfac1 = ( gamma_toshi(mu+4.e0) / gamma_toshi(mu+1.e0) )
4385 lambda = (4.e0/3.e0*cpi*roqr*ccn_out(i)/L_cloud* &
4386 gamfac1)**(1.e0/3.e0) ! [1/cm]
4387 re_cloud_gsfc(i,k) = 1.e0/lambda * gamfac3 * 1.e4 !effective radius [micron]
4388 else
4389 ! when running with WRF_Chem but no aerosol coupling in Goddard MP
4390 if (xland(i) .eq. 1.0) then
4391 ccn_ref = ccn_over_land
4392 else if (xland(i) .eq. 2.0) then
4393 ccn_ref = ccn_over_water
4394 else
4395 print *,' xland is not 1. or 2., run stopped'
4396 ! EMK NUWRF
4397 call wrf_error_fatal(' xland is not 1. or 2., run stopped')
4398 ! stop
4399 endif
4400 ! for cloud water, estimate lambda (slope of gamma distribution)
4401 mu = min(15.e0, (1000.E0/ccn_ref + 2.e0))
4402 gamfac3 = ( gamma_toshi(mu+4.e0) / gamma_toshi(mu+3.e0) )
4403 gamfac1 = ( gamma_toshi(mu+4.e0) / gamma_toshi(mu+1.e0) )
4404 lambda = (4.e0/3.e0*cpi*roqr*ccn_ref/L_cloud* &
4405 gamfac1)**(1.e0/3.e0) ! [1/cm]
4406 re_cloud_gsfc(i,k) = 1.e0/lambda * gamfac3 * 1.e4 !effective radius [micron]
4407 endif ! chem_opt and gsfcgce_gocart_coupling
4408 #else
4409 ! Not running with WRF_Chem
4410 ! ccn_over_land = 1500 ! [#/cm3] climatological value
4411 ! ccn_over_water = 150 ! [#/cm3] climatological value
4412 if (xland(i) .eq. 1.0) then
4413 ccn_ref = ccn_over_land
4414 else if (xland(i) .eq. 2.0) then
4415 ccn_ref = ccn_over_water
4416 else
4417 print *,' xland is not 1. or 2., run stopped'
4418 ! EMK NUWRF
4419 call wrf_error_fatal(' xland is not 1. or 2., run stopped')
4420 ! stop
4421 endif
4422 ! for cloud water, estimate lambda (slope of gamma distribution)
4423 mu = min(15.e0, (1000.E0/ccn_ref + 2.e0))
4424 gamfac3 = ( gamma_toshi(mu+4.e0) / gamma_toshi(mu+3.e0) )
4425 gamfac1 = ( gamma_toshi(mu+4.e0) / gamma_toshi(mu+1.e0) )
4426 lambda = (4.e0/3.e0*cpi*roqr*ccn_ref/L_cloud* &
4427 gamfac1)**(1.e0/3.e0) ! [1/cm]
4428 re_cloud_gsfc(i,k) = 1.e0/lambda * gamfac3 * 1.e4 !effective radius [micron]
4429 #endif
4430 endif ! qcl(i,k) < cmin test
4432 ! for cloud ice
4434 if (qci(i,k) .lt. cmin) then
4435 re_ice_gsfc(i,k) = 0.e0
4436 else
4437 #if (WRF_CHEM == 1)
4438 ! ! when running with WRF_Chem and using aerosol coupling in Goddard MP
4439 ! I_cloud = qci(i,k) * rho(i,k) ! cloud ice [g/cm3]
4440 ! if ( (chem_opt == 300 .or. chem_opt == 301 .or. &
4441 ! chem_opt == 302 .or. chem_opt == 303) .and. &
4442 ! (gsfcgce_gocart_coupling == 1) ) then
4443 ! ! for cloud ice, estimate lambda (slope of gamma distribution)
4444 ! mu = min(15.e0, (1000.E0/(icn_out*1.e-3) + 2.e0))
4445 ! gamfac3 = ( gamma_toshi(mu+4.e0) / gamma_toshi(mu+3.e0) )
4446 ! gamfac1 = ( gamma_toshi(mu+4.e0) / gamma_toshi(mu+1.e0) )
4447 ! lambda = (4.e0/3.e0*cpi*roqi*icn_out*1.e-3/I_cloud* &
4448 ! gamfac1)**(1.e0/3.e0) ! [1/cm]
4449 ! re_ice_gsfc(i,k) = 1.e0/lambda * gamfac3 * 1.e4 !effective radius [micron]
4450 ! else
4451 ! when running with WRF_Chem but no aerosol coupling in Goddard MP
4452 ! for cloud ice effective radius depends on temperature profile, formula from GCE
4453 re_ice_gsfc(i,k) = 125.e0 +(tair(i)-243.16)*5.e0 ! [micron]
4454 if (tair(i) .gt. 243.16) re_ice_gsfc(i,k) = 125.e0
4455 if (tair(i) .lt. 223.16) re_ice_gsfc(i,k) = 25.e0
4456 ! endif ! chem_opt and gsfcgce_gocart_coupling
4457 #else
4458 ! Not running with WRF_Chem
4459 ! for cloud ice effective radius depends on temperature profile, formula from GCE
4460 re_ice_gsfc(i,k) = 125.e0 +(tair(i)-243.16)*5.e0 ! [micron]
4461 if (tair(i) .gt. 243.16) re_ice_gsfc(i,k) = 125.e0
4462 if (tair(i) .lt. 223.16) re_ice_gsfc(i,k) = 25.e0
4463 #endif
4464 endif ! qci(i,k) < cmin test
4466 !JJS 20140305 ^^^^^ Calculate effective radius for all cloud species
4468 ENDDO
4470 1000 continue
4472 !JJS ****************************************************************
4473 !JJS convert from GCE grid back to WRF grid
4474 do k=kts,kte
4475 !dir$ vector aligned
4476 DO i=1,irestrict
4477 ptwrf(i,k) = dpt(i,k)
4478 qvwrf(i,k) = dqv(i,k)
4479 qlwrf(i,k) = qcl(i,k)
4480 qrwrf(i,k) = qrn(i,k)
4481 qiwrf(i,k) = qci(i,k)
4482 qswrf(i,k) = qcs(i,k)
4483 qgwrf(i,k) = qcg(i,k)
4484 qhwrf(i,k) = qch(i,k)
4485 ENDDO
4486 enddo !k
4488 ! ****************************************************************
4490 !+---+-----------------------------------------------------------------+
4491 ! EMK NUWRF...Replace Greg Thompson's dBZ values with those calculated
4492 ! above.
4493 IF ( PRESENT (diagflag) ) THEN
4494 if (diagflag .and. do_radar_ref == 1) then
4495 do k=kts,kte
4496 !dir$ vector aligned
4497 DO i=1,irestrict
4498 refl_10cm(i,k) = max(-35.,dbz(i,k))
4499 ENDDO
4500 end do
4502 endif
4503 ENDIF
4504 !+---+-----------------------------------------------------------------+
4506 END SUBROUTINE saticel_s
4508 SUBROUTINE auto_conversion( L, N, P , re)
4509 implicit none
4510 !-----------------------------------------------------------------------------------------------------
4511 ! Comments:
4512 ! This subroutine compute auto conversion rate folloing Li and Daum [2004], which account for
4513 ! total cloud liquid water, particle number concentrations, PSD lambda, broadening parameters.
4514 !
4515 ! History:
4516 ! 08/2010 Toshi Matsui@NASA GSFC : Initial.
4517 !
4518 !
4519 ! References:
4520 ! Liu, Y. and P. H. Daum, 2004: Parameterization of the autoconversion process. Part I: Analytical
4521 ! formulation of the Kessler-type parameterizations. J. Atmos. Sci, 61, 1539-1548.
4522 !-----------------------------------------------------------------------------------------------------
4523 real,intent(in) :: L ! cloud liquid water [g cm-3]
4524 real,intent(in) :: N ! total number concentration [# cm-3]
4525 real,intent(out) :: P ! auto conversion rate [g cm-3 s-1]
4526 real,intent(out) :: re ! cloud effective radius [micron]
4528 real :: mu ! mu of gamma PSD [-]
4529 real :: eta ! eta function [cm3 g-2 s-1]
4530 real :: beta, beta1, beta2 ! beta function [-]
4531 real :: gamfac , gfac1 , gfac2 ! gamma function [-]
4532 real :: R6_6power ! mean radius of the sixth moment [cm]
4533 real :: R6_thresh ! threshold of mean radius of the sixth moment [cm]
4534 real :: R6 ! mean radius of the sixth moment [cm]
4535 real :: Heaviside_func ! Heaviside step function (0 or 1)
4536 real :: lambda ! slope of gamma size ditribution [1/cm]
4537 ! real :: No ! intercept [cm-4]
4539 real,parameter :: Rc = 10.e0 * 1.e-4 ! threshold of particle radus (10 micron) [cm]
4540 real,parameter :: const_pi = 3.14159e0 ! pai
4541 real,parameter :: const_kappa = 1.9e11 ! coefficient for water droplet collection kernel [cm-3 s-1]
4542 ! from Long [1974, JAS].
4543 ! real,parameter :: const_kappa = 1.9e11*10000.e0 !10000 is to adjust the order to keseller
4546 real,parameter :: const_rho_liq = 1.e0 ! density of liquid water [g cm-3]
4547 real,parameter :: eta_func = ((3.e0/(4.e0*const_pi*const_rho_liq))**2) * const_kappa ! eta function [cm3 g-2 s-1]
4548 ! a part of (eq 27b)
4550 logical,parameter :: no_thresh = .true. ! logic to choose no threshold parameterization or not.
4552 !
4553 ! When no particel, no autoconversion.
4554 !
4555 ! EMK BUG FIX...Prevent overflow for small but non-zero values of L
4556 ! if( N <= 0.e0 .or. L <= 0.e0 ) then
4557 if( N <= 0.e0 .or. L <= 1.0e-32 ) then
4558 P = 0.e0
4559 return
4560 endif
4562 !
4563 ! check bad values of N and L
4564 !
4565 ! if( N < 0.e0 ) stop 'MSG auto_conversion: N is negative, it must be positive'
4566 ! if( L < 0.e0 ) stop 'MSG auto_conversion: L is negative, it must be positive'
4567 ! if( L > 1.e0 ) stop 'MSG auto_conversion: L is greater than 1g/cm3.'
4569 !
4570 ! Empirical fit of mu as a function of total particle number concentrations.
4571 ! From Martin et al. (1994), assign gamma shape parameter mu for cloud
4572 ! drops according to general dispersion characteristics.
4573 ! disp=~0.25 for Maritime and 0.45 for Continental.
4574 ! Since disp=SQRT((mu+2)/(mu+1) - 1), mu varies from 15 for Maritime (pristine air)
4575 ! to 2 for Continental (really dirty air). if mu = 0 --> expnential distribution (narrow dist)
4576 !
4578 ! orig
4579 mu = MIN(15.e0, (1000.E0/N + 2.e0))
4582 !
4583 ! gamma functions
4584 !
4585 gfac1 = gamma_toshi(mu+4.e0)
4586 gfac2 = gamma_toshi(mu+1.e0)
4587 gamfac = (gfac1/gfac2)
4589 !
4590 ! estimate lambda (slope of gamma distribution)
4591 !
4592 lambda = (4.e0/3.e0*const_pi*const_rho_liq*N/L*gamfac)**(1.e0/3.e0) ! [1/cm]
4595 THRESH: if( no_thresh ) then !-------------------------------------------
4597 !
4598 ! threshold of particle radius (mean radius of the sixth moment )
4599 !
4600 gfac1 = gamma_toshi(mu+7.e0)
4601 gfac2 = gamma_toshi(mu+1.e0)
4602 gamfac = (gfac1/gfac2)
4604 R6_6power = (1.e0 / lambda)**6.e0 * gamfac ![cm] (eq. A3)
4606 !
4607 ! auto conversion rate (eq. 26a)
4608 !
4609 P = const_kappa * N * R6_6power * L ! [g cm-3 s-1 ]
4610 ! [cm-3 s-1] * [#/cm3] * [cm6] * [g/cm3]
4613 else !with threshold ---------------------------------------------------
4615 !
4616 ! Estimate eta under gamma PSD
4617 !
4618 beta1 = (6.e0+mu)*(5.e0+mu)*(4.e0+mu)
4619 beta2 = (3.e0+mu)*(2.e0+mu)*(1.e0+mu)
4620 beta = beta1 / beta2
4622 eta = eta_func * beta ! eta function (eq 27b) [cm3 g-2 s-1]
4624 !
4625 ! threshold of particle radius (mean radius of the sixth moment )
4626 !
4627 R6_thresh = beta * Rc ![cm] (pg 1545)
4629 !
4630 ! mean radius of the sixth moment
4631 !
4632 gfac1 = gamma_toshi(6.e0+mu+1.e0)
4633 gfac2 = gamma_toshi(1.e0+mu)
4634 gamfac = (gfac1/gfac2)**(1.e0/6.e0)
4636 R6 = (1.e0 / lambda) * gamfac ![cm] (eq. A3)
4638 !
4639 ! Heaviside step function
4640 !
4641 if ( R6 - R6_thresh <= 0.e0 ) then
4642 Heaviside_func = 0.e0
4643 elseif( R6 - R6_thresh > 0.e0 ) then
4644 Heaviside_func = 1.e0
4645 else
4646 ! NUWRF EMK...User WRF's library to gracefully stop MPI.
4647 write(wrf_err_message,*)'MSG: auto_conversion: Strange value of R6= ',R6
4648 call wrf_error_fatal(trim(wrf_err_message))
4649 ! print*, 'MSG: auto_conversion: Strange value of R6= ', R6 ; stop
4650 endif
4652 !
4653 ! auto conversion rate [g cm-3 s-1 ] (eq. 27a)
4654 !
4655 P = eta * (1.e0/N) * (L**3) * Heaviside_func
4657 ! [cm3 g-2 s-1] * [cm3] * [g3/cm9]
4659 endif THRESH !------------------------------------------------------------
4662 ! optional
4664 !
4665 ! estimate effective radius
4666 !
4667 gfac1 = gamma_toshi(mu+4.e0)
4668 gfac2 = gamma_toshi(mu+3.e0)
4669 gamfac = (gfac1/gfac2)
4671 re = 1.e0/lambda * gamfac * 1.e4 !effective radius [micron]
4673 !
4674 ! estimate No
4675 !
4676 ! call gamma_function(mu+1.e0 ,gfac1)
4677 ! No = N * (lambda**(mu+1)) / gfac1
4679 END subroutine auto_conversion
4681 !DIR$ ATTRIBUTES FORCEINLINE :: gamma_toshi
4682 real function gamma_toshi(x)
4684 !---------------------------------------------------------------------------------------------------
4685 ! Comments:
4686 ! compute the gamma function T(x) for single precision floating point.
4687 ! input : x --- argument of a(x)
4688 ! ( x is not equal to 0,-1,-2,... )
4689 !
4690 ! History:
4691 ! 09/2009 Toshi Matsui@NASA GSFC ; Adapted to SDSU
4692 !
4693 ! References:
4694 !----------------------------------------------------------------------------------------------------
4695 implicit double precision (a-h,o-z)
4696 dimension g(26)
4697 data g/1.0d0,0.5772156649015329d0, &
4698 -0.6558780715202538d0, -0.420026350340952d-1, &
4699 0.1665386113822915d0,-.421977345555443d-1, &
4700 -.96219715278770d-2, .72189432466630d-2, &
4701 -.11651675918591d-2, -.2152416741149d-3, &
4702 .1280502823882d-3, -.201348547807d-4, &
4703 -.12504934821d-5, .11330272320d-5, &
4704 -.2056338417d-6, .61160950d-8, &
4705 .50020075d-8, -.11812746d-8, &
4706 .1043427d-9, .77823d-11, &
4707 -.36968d-11, .51d-12, &
4708 -.206d-13, -.54d-14, .14d-14, .1d-15/
4709 real :: x
4711 pi=3.141592653589793d0
4712 if (x.eq.int(x)) then
4713 if (x.gt.0.0d0) then
4714 ga=1.0d0
4715 m1=int(x)-1
4716 do k=2,m1
4717 ga=ga*k
4718 enddo
4719 else
4720 ga=1.0d+300
4721 endif
4722 else
4723 if (dabs(dble(x)).gt.1.0d0) then
4724 z=dabs(dble(x))
4725 m=int(z)
4726 r=1.0d0
4727 do k=1,m
4728 r=r*(z-k)
4729 enddo
4730 z=z-m
4731 else
4732 z=dble(x)
4733 endif
4734 gr=g(26)
4735 do k=25,1,-1
4736 gr=gr*z+g(k)
4737 enddo
4738 ga=1.0d0/(gr*z)
4739 if (dabs(dble(x)).gt.1.0d0) then
4740 ga=ga*r
4741 if (x.lt.0.0d0) ga=-pi/(x*ga*dsin(pi*x))
4742 endif
4743 endif
4745 gamma_toshi = real(ga)
4747 end function gamma_toshi
4749 SUBROUTINE Find_NaN_Inf_Double(Warning_MSG, real_input, i_in,j_in,k_in)
4750 implicit none
4752 real,intent(inout) :: real_input !anykind of Non-dimensional input Real parameters
4753 integer,intent(in) :: i_in, j_in, k_in
4754 character*(*),intent(in) :: Warning_MSG
4756 !
4757 ! Find Infinity
4758 !
4759 !if( exp(-abs(real_input)) == 0.) then ! this formulae is bit slow
4761 if( 1e+10/real_input == 0. ) then
4762 print*,'MSG Find_NaN_Inf: '//Warning_MSG//'Infinity at',i_in,j_in,k_in
4763 real_input = 0.
4764 return
4765 endif
4767 !
4768 ! Find NaN
4769 !
4770 if( real_input==0. .or. real_input>0. .or. real_input<0. .or. real_input>=0. .or. real_input<=0. ) then
4771 else
4772 print*,'MSG Find_NaN_Inf: '//Warning_MSG//'NaN at',i_in,j_in,k_in
4773 real_input = 0.
4774 return
4775 endif
4777 END SUBROUTINE Find_NaN_Inf_Double
4779 !+---+-----------------------------------------------------------------+
4781 subroutine refl10cm_gsfc (qv1d, qr1d, qs1d, qg1d, &
4782 t1d, p1d, dBZ, kts, kte, ii, jj)
4784 IMPLICIT NONE
4786 !..Sub arguments
4787 INTEGER, INTENT(IN):: kts, kte, ii, jj
4788 REAL, DIMENSION(kts:kte), INTENT(IN):: &
4789 qv1d, qr1d, qs1d, qg1d, t1d, p1d
4790 REAL, DIMENSION(kts:kte), INTENT(INOUT):: dBZ
4792 !..Local variables
4793 REAL, DIMENSION(kts:kte):: temp, pres, qv, rho
4794 REAL, DIMENSION(kts:kte):: rr, rs, rg
4796 DOUBLE PRECISION, DIMENSION(kts:kte):: ilamr, ilams, ilamg
4797 DOUBLE PRECISION, DIMENSION(kts:kte):: N0_r, N0_s, N0_g
4798 DOUBLE PRECISION:: lamr, lams, lamg
4799 LOGICAL, DIMENSION(kts:kte):: L_qr, L_qs, L_qg
4801 REAL, DIMENSION(kts:kte):: ze_rain, ze_snow, ze_graupel
4802 DOUBLE PRECISION:: fmelt_s, fmelt_g
4804 INTEGER:: i, k, k_0, kbot, n
4805 LOGICAL:: melti
4807 DOUBLE PRECISION:: cback, x, eta, f_d
4808 REAL, PARAMETER:: R=287.
4809 REAL, PARAMETER:: PIx=3.1415926536
4811 !+---+
4813 do k = kts, kte
4814 dBZ(k) = -35.0
4815 enddo
4817 !+---+-----------------------------------------------------------------+
4818 !..Put column of data into local arrays.
4819 !+---+-----------------------------------------------------------------+
4820 do k = kts, kte
4821 temp(k) = t1d(k)
4822 qv(k) = MAX(1.E-10, qv1d(k))
4823 pres(k) = p1d(k)
4824 rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622))
4826 if (qr1d(k) .gt. 1.E-9) then
4827 rr(k) = qr1d(k)*rho(k)
4828 N0_r(k) = xnor
4829 lamr = (xam_r*xcrg(3)*N0_r(k)/rr(k))**(1./xcre(1))
4830 ilamr(k) = 1./lamr
4831 L_qr(k) = .true.
4832 else
4833 rr(k) = 1.E-12
4834 L_qr(k) = .false.
4835 endif
4837 if (qs1d(k) .gt. 1.E-9) then
4838 rs(k) = qs1d(k)*rho(k)
4839 N0_s(k) = xnos
4840 lams = (xam_s*xcsg(3)*N0_s(k)/rs(k))**(1./xcse(1))
4841 ilams(k) = 1./lams
4842 L_qs(k) = .true.
4843 else
4844 rs(k) = 1.E-12
4845 L_qs(k) = .false.
4846 endif
4848 if (qg1d(k) .gt. 1.E-9) then
4849 rg(k) = qg1d(k)*rho(k)
4850 N0_g(k) = xnog
4851 lamg = (xam_g*xcgg(3)*N0_g(k)/rg(k))**(1./xcge(1))
4852 ilamg(k) = 1./lamg
4853 L_qg(k) = .true.
4854 else
4855 rg(k) = 1.E-12
4856 L_qg(k) = .false.
4857 endif
4858 enddo
4860 !+---+-----------------------------------------------------------------+
4861 !..Locate K-level of start of melting (k_0 is level above).
4862 !+---+-----------------------------------------------------------------+
4863 melti = .false.
4864 k_0 = kts
4865 do k = kte-1, kts, -1
4866 if ( (temp(k).gt.273.15) .and. L_qr(k) &
4867 .and. (L_qs(k+1).or.L_qg(k+1)) ) then
4868 k_0 = MAX(k+1, k_0)
4869 melti=.true.
4870 goto 195
4871 endif
4872 enddo
4873 195 continue
4875 !+---+-----------------------------------------------------------------+
4876 !..Assume Rayleigh approximation at 10 cm wavelength. Rain (all temps)
4877 !.. and non-water-coated snow and graupel when below freezing are
4878 !.. simple. Integrations of m(D)*m(D)*N(D)*dD.
4879 !+---+-----------------------------------------------------------------+
4881 do k = kts, kte
4882 ze_rain(k) = 1.e-22
4883 ze_snow(k) = 1.e-22
4884 ze_graupel(k) = 1.e-22
4885 if (L_qr(k)) ze_rain(k) = N0_r(k)*xcrg(4)*ilamr(k)**xcre(4)
4886 if (L_qs(k)) ze_snow(k) = (0.176/0.93) * (6.0/PIx)*(6.0/PIx) &
4887 * (xam_s/900.0)*(xam_s/900.0) &
4888 * N0_s(k)*xcsg(4)*ilams(k)**xcse(4)
4889 if (L_qg(k)) ze_graupel(k) = (0.176/0.93) * (6.0/PIx)*(6.0/PIx) &
4890 * (xam_g/900.0)*(xam_g/900.0) &
4891 * N0_g(k)*xcgg(4)*ilamg(k)**xcge(4)
4892 enddo
4895 !+---+-----------------------------------------------------------------+
4896 !..Special case of melting ice (snow/graupel) particles. Assume the
4897 !.. ice is surrounded by the liquid water. Fraction of meltwater is
4898 !.. extremely simple based on amount found above the melting level.
4899 !.. Uses code from Uli Blahak (rayleigh_soak_wetgraupel and supporting
4900 !.. routines).
4901 !+---+-----------------------------------------------------------------+
4903 if (melti .and. k_0.ge.kts+1) then
4904 do k = k_0-1, kts, -1
4906 !..Reflectivity contributed by melting snow
4907 if (L_qs(k) .and. L_qs(k_0) ) then
4908 fmelt_s = MAX(0.005d0, MIN(1.0d0-rs(k)/rs(k_0), 0.99d0))
4909 eta = 0.d0
4910 lams = 1./ilams(k)
4911 do n = 1, nrbins
4912 x = xam_s * xxDs(n)**xbm_s
4913 call rayleigh_soak_wetgraupel (x,DBLE(xocms),DBLE(xobms), &
4914 fmelt_s, melt_outside_s, m_w_0, m_i_0, lamda_radar, &
4915 CBACK, mixingrulestring_s, matrixstring_s, &
4916 inclusionstring_s, hoststring_s, &
4917 hostmatrixstring_s, hostinclusionstring_s)
4918 f_d = N0_s(k)*xxDs(n)**xmu_s * DEXP(-lams*xxDs(n))
4919 eta = eta + f_d * CBACK * simpson(n) * xdts(n)
4920 enddo
4921 ze_snow(k) = SNGL(lamda4 / (pi5 * K_w) * eta)
4922 endif
4924 !..Reflectivity contributed by melting graupel
4926 if (L_qg(k) .and. L_qg(k_0) ) then
4927 fmelt_g = MAX(0.005d0, MIN(1.0d0-rg(k)/rg(k_0), 0.99d0))
4928 eta = 0.d0
4929 lamg = 1./ilamg(k)
4930 do n = 1, nrbins
4931 x = xam_g * xxDg(n)**xbm_g
4932 call rayleigh_soak_wetgraupel (x,DBLE(xocmg),DBLE(xobmg), &
4933 fmelt_g, melt_outside_g, m_w_0, m_i_0, lamda_radar, &
4934 CBACK, mixingrulestring_g, matrixstring_g, &
4935 inclusionstring_g, hoststring_g, &
4936 hostmatrixstring_g, hostinclusionstring_g)
4937 f_d = N0_g(k)*xxDg(n)**xmu_g * DEXP(-lamg*xxDg(n))
4938 eta = eta + f_d * CBACK * simpson(n) * xdtg(n)
4939 enddo
4940 ze_graupel(k) = SNGL(lamda4 / (pi5 * K_w) * eta)
4941 endif
4943 enddo
4944 endif
4946 do k = kte, kts, -1
4947 dBZ(k) = 10.*log10((ze_rain(k)+ze_snow(k)+ze_graupel(k))*1.d18)
4948 enddo
4950 end subroutine refl10cm_gsfc
4952 !+---+-----------------------------------------------------------------+
4954 !JJS 20140225
4955 ! Calculate cloud droplet effective radius
4956 real function eff_rad(lambda)
4958 #ifndef NO_IEEE_MODULE
4959 use, intrinsic :: ieee_arithmetic
4960 #endif
4961 implicit none
4963 !---------------------------------------------------------------------------------------------------
4964 ! Comments:
4965 ! Compute drop effective radius from slope parameters (lambda) of expoential size distribution.
4966 !
4967 ! History:
4968 ! 02/2014 Toshi Matsui@NASA GSFC ; Initial
4969 !
4970 ! References:
4971 !----------------------------------------------------------------------------------------------------
4972 real,intent(in) :: lambda ! intercept parameter [1/cm]
4974 !
4975 ! for no particles.
4976 !
4977 #ifndef NO_IEEE_MODULE
4978 if ( lambda <= 0.e0 .or. ieee_is_nan(lambda) ) then
4979 #else
4980 if ( lambda <= 0.e0 ) then
4981 #endif
4982 eff_rad = 0.e0
4983 return
4984 endif
4986 !
4987 ! compute drop effective radius for exponential distribution N(D) = N0*exp(-lam*D)
4988 !
4989 eff_rad = 1.5e0 / (lambda*100.) * 1.0e+6 ! [micron]
4991 end function eff_rad
4993 END MODULE module_mp_gsfcgce_4ice_nuwrf