| blob | 0afa3434f6e0351be97650c4ca0ee6723622e9ce |
1 MODULE module_mosaic_wetscav
3 CONTAINS
5 !===========================================================================
6 !===========================================================================
7 subroutine wetscav_cbmz_mosaic (id,ktau,dtstep,ktauc,config_flags, &
8 dtstepc,alt,t_phy,p8w,t8w,p_phy,chem,rho_phy,cldfra, &
9 qlsink,precr,preci,precs,precg,qsrflx, &
10 gas_aqfrac, numgas_aqfrac, &
11 ids,ide, jds,jde, kds,kde, &
12 ims,ime, jms,jme, kms,kme, &
13 its,ite, jts,jte, kts,kte )
15 ! wet removal by grid-resolved precipitation
16 ! scavenging of cloud-phase aerosols and gases by collection, freezing, ...
17 ! scavenging of interstitial-phase aerosols by impaction
18 ! scavenging of gas-phase gases by mass transfer and reaction
20 !----------------------------------------------------------------------
21 USE module_configure, only: grid_config_rec_type
22 USE module_state_description
23 USE module_data_mosaic_asect
25 !----------------------------------------------------------------------
26 IMPLICIT NONE
28 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
30 INTEGER, INTENT(IN ) :: &
31 ids,ide, jds,jde, kds,kde, &
32 ims,ime, jms,jme, kms,kme, &
33 its,ite, jts,jte, kts,kte, &
34 id, ktau, ktauc, numgas_aqfrac
35 REAL, INTENT(IN ) :: dtstep,dtstepc
36 !
37 ! all advected chemical species
38 !
39 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), &
40 INTENT(INOUT ) :: chem
42 ! fraction of gas species in cloud water
43 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, numgas_aqfrac ), &
44 INTENT(IN ) :: gas_aqfrac
46 !
47 !
48 ! input from meteorology
49 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
50 INTENT(IN ) :: &
51 alt, &
52 t_phy, &
53 p_phy, &
54 t8w,p8w, &
55 qlsink,precr,preci,precs,precg, &
56 rho_phy,cldfra
57 REAL, DIMENSION( ims:ime, jms:jme, num_chem ), &
58 INTENT(OUT ) :: qsrflx ! column change due to scavening
60 call wetscav (id,ktau,dtstep,ktauc,config_flags, &
61 dtstepc,alt,t_phy,p8w,t8w,p_phy,chem,rho_phy,cldfra, &
62 qlsink,precr,preci,precs,precg,qsrflx, &
63 gas_aqfrac, numgas_aqfrac, &
64 ntype_aer, nsize_aer, ncomp_aer, &
65 massptr_aer, dens_aer, numptr_aer, &
66 maxd_acomp, maxd_asize,maxd_atype, maxd_aphase, ai_phase, cw_phase, &
67 volumcen_sect, volumlo_sect, volumhi_sect, &
68 waterptr_aer, dens_water_aer, &
69 scavimptblvol, scavimptblnum, nimptblgrow_mind, nimptblgrow_maxd, dlndg_nimptblgrow, &
70 ids,ide, jds,jde, kds,kde, &
71 ims,ime, jms,jme, kms,kme, &
72 its,ite, jts,jte, kts,kte )
74 end subroutine wetscav_cbmz_mosaic
77 !===========================================================================
78 !===========================================================================
79 subroutine wetscav_mozart_mosaic (id,ktau,dtstep,ktauc,config_flags, &
80 dtstepc,alt,t_phy,p8w,t8w,p_phy,chem,rho_phy,cldfra, &
81 qlsink,precr,preci,precs,precg,qsrflx, &
82 gas_aqfrac, numgas_aqfrac, &
83 ids,ide, jds,jde, kds,kde, &
84 ims,ime, jms,jme, kms,kme, &
85 its,ite, jts,jte, kts,kte )
87 ! wet removal by grid-resolved precipitation
88 ! scavenging of cloud-phase aerosols and gases by collection, freezing, ...
89 ! scavenging of interstitial-phase aerosols by impaction
90 ! scavenging of gas-phase gases by mass transfer and reaction
92 !----------------------------------------------------------------------
93 USE module_configure, only: grid_config_rec_type
94 USE module_state_description
95 USE module_data_mosaic_asect
97 !----------------------------------------------------------------------
98 IMPLICIT NONE
100 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
102 INTEGER, INTENT(IN ) :: &
103 ids,ide, jds,jde, kds,kde, &
104 ims,ime, jms,jme, kms,kme, &
105 its,ite, jts,jte, kts,kte, &
106 id, ktau, ktauc, numgas_aqfrac
107 REAL, INTENT(IN ) :: dtstep,dtstepc
108 !
109 ! all advected chemical species
110 !
111 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), &
112 INTENT(INOUT ) :: chem
114 ! fraction of gas species in cloud water
115 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, numgas_aqfrac ), &
116 INTENT(IN ) :: gas_aqfrac
118 !
119 !
120 ! input from meteorology
121 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
122 INTENT(IN ) :: &
123 alt, &
124 t_phy, &
125 p_phy, &
126 t8w,p8w, &
127 qlsink,precr,preci,precs,precg, &
128 rho_phy,cldfra
129 REAL, DIMENSION( ims:ime, jms:jme, num_chem ), &
130 INTENT(OUT ) :: qsrflx ! column change due to scavening
132 call wetscav (id,ktau,dtstep,ktauc,config_flags, &
133 dtstepc,alt,t_phy,p8w,t8w,p_phy,chem,rho_phy,cldfra, &
134 qlsink,precr,preci,precs,precg,qsrflx, &
135 gas_aqfrac, numgas_aqfrac, &
136 ntype_aer, nsize_aer, ncomp_aer, &
137 massptr_aer, dens_aer, numptr_aer, &
138 maxd_acomp, maxd_asize,maxd_atype, maxd_aphase, ai_phase, cw_phase, &
139 volumcen_sect, volumlo_sect, volumhi_sect, &
140 waterptr_aer, dens_water_aer, &
141 scavimptblvol, scavimptblnum, nimptblgrow_mind, nimptblgrow_maxd, dlndg_nimptblgrow, &
142 ids,ide, jds,jde, kds,kde, &
143 ims,ime, jms,jme, kms,kme, &
144 its,ite, jts,jte, kts,kte )
146 end subroutine wetscav_mozart_mosaic
148 !===========================================================================
149 !===========================================================================
150 subroutine wetscav (id,ktau,dtstep,ktauc,config_flags, &
151 dtstepc,alt,t_phy,p8w,t8w,p_phy,chem,rho_phy,cldfra, &
152 qlsink,precr,preci,precs,precg,qsrflx, &
153 gas_aqfrac, numgas_aqfrac, &
154 ntype_aer, nsize_aer, ncomp_aer, &
155 massptr_aer, dens_aer, numptr_aer, &
156 maxd_acomp, maxd_asize,maxd_atype, maxd_aphase, ai_phase, cw_phase, &
157 volumcen_sect, volumlo_sect, volumhi_sect, &
158 waterptr_aer, dens_water_aer, &
159 scavimptblvol, scavimptblnum, nimptblgrow_mind, nimptblgrow_maxd, dlndg_nimptblgrow, &
160 ids,ide, jds,jde, kds,kde, &
161 ims,ime, jms,jme, kms,kme, &
162 its,ite, jts,jte, kts,kte )
164 ! wet removal by grid-resolved precipitation
165 ! scavenging of cloud-phase aerosols and gases by collection, freezing, ...
166 ! scavenging of interstitial-phase aerosols by impaction
167 ! scavenging of gas-phase gases by mass transfer and reaction
169 !----------------------------------------------------------------------
170 USE module_configure, only: grid_config_rec_type
171 USE module_state_description
172 USE module_model_constants, only: g, rhowater, mwdry
173 USE module_dep_simple, only: is_aerosol
175 IMPLICIT NONE
177 !======================================================================
179 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
181 INTEGER, INTENT(IN ) :: &
182 ids,ide, jds,jde, kds,kde, &
183 ims,ime, jms,jme, kms,kme, &
184 its,ite, jts,jte, kts,kte, &
185 id, ktau, ktauc, numgas_aqfrac
186 REAL, INTENT(IN ) :: dtstep,dtstepc
187 !
188 ! all advected chemical species
189 !
190 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), &
191 INTENT(INOUT ) :: chem
193 ! fraction of gas species in cloud water
194 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, numgas_aqfrac ), &
195 INTENT(IN ) :: gas_aqfrac
197 !
198 ! input from meteorology
199 !
200 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
201 INTENT(IN ) :: &
202 alt, &
203 t_phy, &
204 p_phy, &
205 t8w,p8w, &
206 qlsink,precr,preci,precs,precg, &
207 rho_phy,cldfra
208 integer, intent(in) :: maxd_atype, maxd_asize, maxd_acomp, maxd_aphase
209 integer, intent(in) :: ai_phase, cw_phase
210 integer, intent(in) :: ntype_aer
211 integer, intent(in) :: nsize_aer( maxd_atype ), & ! number of size bins
212 ncomp_aer( maxd_atype ), & ! number of chemical components
213 massptr_aer( maxd_acomp, maxd_asize, maxd_atype, maxd_aphase), & ! index for mixing ratio
214 waterptr_aer( maxd_asize, maxd_atype ), & ! index for aerosol water
215 numptr_aer( maxd_asize, maxd_atype, maxd_aphase ) ! index for the number mixing ratio
216 real, intent(in) :: dens_aer( maxd_acomp, maxd_atype ), & ! material density (g/cm3)
217 dens_water_aer ! water density (g/m3)
218 real, intent(in) :: volumcen_sect( maxd_asize, maxd_atype ), & ! single-particle volumes (cm3)
219 volumlo_sect( maxd_asize, maxd_atype ), & ! at section center, lower boundary,
220 volumhi_sect( maxd_asize, maxd_atype ) ! upper boundary
222 real, intent(in) :: dlndg_nimptblgrow
223 integer, intent(in) :: nimptblgrow_mind, nimptblgrow_maxd
224 real, intent(in) :: scavimptblnum(4, nimptblgrow_mind:nimptblgrow_maxd, maxd_asize, maxd_atype), &
225 scavimptblvol(4, nimptblgrow_mind:nimptblgrow_maxd, maxd_asize, maxd_atype)
226 REAL, DIMENSION( ims:ime, jms:jme, num_chem ), &
227 INTENT(OUT ) :: qsrflx ! column change due to scavening
229 !
230 ! LOCAL VAR
231 !
232 integer :: i,j,k,l,m,n
233 integer :: lmasscw,lnumcw
234 real scale
235 logical :: isprx(ims:ime,kms:kme,jms:jme)
237 real :: pdel(ims:ime,kms:kme,jms:jme)
238 real :: delpf(kms:kme)
239 real :: delpfhalf
240 real :: dump
241 real :: fac_pwght_ls(kms:kme) !
242 real :: fapincld, fapoutcld, faptot
243 real :: fapmin_ls !
244 real :: fapx(ims:ime,kms:kme,jms:jme)
245 real :: fracscav
246 real :: pf_above, pf_below, pdel_fac
247 real :: pf_ls(kms:kme)
248 real :: pfoutcld
249 real :: pfsmall_ls ! l-s precip fluxes (kg/m2/s) smaller than this
250 ! are ignored (treated as zero)
251 real :: pfsmall_min
252 real :: pfx(ims:ime,kms:kme,jms:jme)
253 real :: pfx_inrain(ims:ime,kms:kme,jms:jme)
254 real :: sumfa, sumpf, sumpffa
255 REAL :: dqdt( ims:ime, kms:kme, jms:jme, num_chem )
256 logical dotend(num_chem)
258 dotend(:) = .false.
259 dqdt(:,:,:,:) = 0.0
260 qsrflx(:,:,:) = 0.0
263 ! scavenging cloud-phase aerosol assuming precip falls to surface
264 do 100 j=jts,jte
265 do k=kts,kte
266 do i=its,ite
267 pdel(i,k,j)=p8w(i,k,j)-p8w(i,k+1,j)
268 end do
269 end do
271 do 100 k=kts,kte
272 do 100 i=its,ite
273 scale=1.0-dtstepc*qlsink(i,k,j)
274 scale=max(0.0,min(1.0,scale)) ! make sure 0 <= scale <= 1
275 if (qlsink(i,k,j) > 0.0) then
276 pdel_fac = pdel(i,k,j)/g
277 do n=1,ntype_aer
278 do m=1,nsize_aer(n)
279 do l=1,ncomp_aer(n)
280 lmasscw=massptr_aer(l,m,n,cw_phase)
281 if (lmasscw < param_first_scalar) cycle
282 qsrflx(i,j,lmasscw)=qsrflx(i,j,lmasscw)+chem(i,k,j,lmasscw)*(scale-1.)*pdel_fac ! aerosol mass (ug/m2)
283 chem(i,k,j,lmasscw)=chem(i,k,j,lmasscw)*scale
284 end do ! comp
285 lnumcw=numptr_aer(m,n,cw_phase)
286 if (lnumcw < param_first_scalar) cycle
287 qsrflx(i,j,lnumcw)=qsrflx(i,j,lnumcw)+chem(i,k,j,lnumcw)*(scale-1.)*pdel_fac ! aerosol number (1/m2)
288 chem(i,k,j,lnumcw)=chem(i,k,j,lnumcw)*scale
289 end do ! size
290 end do ! type
291 end if ! qlsink > 0
292 100 continue
295 ! scavenging of gases in cloud-water
296 ! only if not MOZART (--> Neu and Prather used)
297 if ( .NOT. (config_flags%chem_opt == mozart_mosaic_4bin_aq_kpp) ) then
298 do 290 l = param_first_scalar, min( num_chem, numgas_aqfrac )
299 if ( is_aerosol(l) ) goto 290
300 do 270 j = jts,jte
301 do 270 k = kts,kte
302 do 270 i = its,ite
303 fracscav = dtstepc*qlsink(i,k,j)*gas_aqfrac(i,k,j,l)
304 if (fracscav > 0.0) then ! make sure fracscav > 0
305 fracscav = min(1.0,fracscav) ! make sure fracscav <= 1
306 scale = 1.0 - fracscav
307 pdel_fac = pdel(i,k,j)/(g*mwdry)
308 qsrflx(i,j,l) = qsrflx(i,j,l)+chem(i,k,j,l)*(scale-1.)*pdel_fac ! mmol/m2
309 chem(i,k,j,l) = chem(i,k,j,l)*scale
310 end if
311 270 continue
312 290 continue
313 end if
316 ! below-cloud scavenging
318 ! precip rates -- 1.0 kgwtr/m2/s = 1.0e-3 m3wtr/m2/s = 1.0e-3 m/s
319 ! 7.06e-5 kg/m2/s = 7.06e-8 m/s = 0.01 inch/h
320 ! 3.17e-5 kg/m2/s = 3.17e-8 m/s = 1 m/y = global annual average
321 ! 1.00e-7 kg/m2/s = 1.00e-10 m/s is a very small precip rate!
323 fapmin_ls = 1.0e-3
324 pfsmall_ls = 1.0e-7
326 isprx(:,:,:) = .false.
327 pfx(:,:,:) = 0.0
328 pfx_inrain(:,:,:) = 0.0
329 fapx(:,:,:) = 0.0
332 do 5900 j=jts,jte
333 do 5900 i=its,ite
335 !----------------------------------------------------------------------
336 ! compute l-s precip rates
337 !
338 ! pf_ls(k) = precip flux at center of level
339 pf_below = 0.0
340 do k = kte,kts,-1
341 pf_above = pf_below
342 pf_below = precr(i,k,j) + preci(i,k,j) + precs(i,k,j) + precg(i,k,j) ! total precip rate
343 if (pf_below < pfsmall_ls) pf_below = 0.0
344 delpf(k) = pf_below - pf_above
345 pf_ls(k) = 0.5*(pf_below + pf_above)
346 if (pf_ls(k) < pfsmall_ls) pf_ls(k) = 0.0
347 end do
349 ! compute fac_pwght_ls which is an average of cloud fractional area in and
350 ! above the current level, weighted by precip production in each level
351 ! basically this reflect the cloud area at levels where precip is produced
352 do k = kte, kts,-1
353 if (k == kte) then
354 ! compute change from (k-1/k) interface to level k center
355 delpfhalf = 0.5*delpf(k)
356 if (delpfhalf > 0.0) then
357 fac_pwght_ls(k) = max( cldfra(i,k,j), fapmin_ls )
358 sumpffa = delpfhalf * fac_pwght_ls(k)
359 sumpf = delpfhalf
360 else
361 fac_pwght_ls(k) = fapmin_ls
362 sumpffa = 0.0
363 sumpf = 0.0
364 end if
365 else
366 ! compute change from level (k+1) center to (k+1/k) interface
367 delpfhalf = 0.5*delpf(k+1)
368 if (delpfhalf > 0.0) then
369 sumpffa = sumpffa + delpfhalf*max( cldfra(i,k+1,j), fapmin_ls )
370 sumpf = sumpf + delpfhalf
371 fac_pwght_ls(k) = max( (sumpffa/sumpf), fapmin_ls )
372 else
373 fac_pwght_ls(k) = fac_pwght_ls(k+1)
374 sumpffa = max( (sumpffa + delpfhalf*fac_pwght_ls(k)), 0.0 )
375 sumpf = max( (sumpf + delpfhalf), 0.0 )
376 end if
377 ! compute change from (k-1/k) interface to level k center
378 delpfhalf = 0.5*delpf(k)
379 if (delpfhalf > 0.0) then
380 sumpffa = sumpffa + delpfhalf*max( cldfra(i,k,j), fapmin_ls )
381 sumpf = sumpf + delpfhalf
382 fac_pwght_ls(k) = max( (sumpffa/sumpf), fapmin_ls )
383 else
384 ! here, fac_pwght_ls(k) is unchanged from its value computed just above
385 sumpffa = max( (sumpffa + delpfhalf*fac_pwght_ls(k)), 0.0 )
386 sumpf = max( (sumpf + delpfhalf), 0.0 )
387 end if
388 end if
389 end do
392 !----------------------------------------------------------------------
393 !----------------------------------------------------------------------
394 ! do the scavenging
395 !
396 !----------------------------------------------------------------------
397 ! loop through levels
398 !
399 do 4900 k = kte, kts,-1
402 !----------------------------------------------------------------------
403 ! for each cloud/precip type (ls, dp, sh), compute
404 ! fapx = fractional area with precip
405 ! pfx = precip flux based on entire grid-cell area (kg/m2/s)
406 ! pfx_inrain = precip flux within the precip fractional area (kg/m2/s)
407 !
408 sumpf = 0.0
409 sumfa = 0.0
411 ! l-s cloud
412 ! assume faptot = total (in + out-of-cloud) precip area
413 ! = 0.5*fac_pwght_ls(k)
414 ! then fapincld = in-cloud precip area
415 ! = min( faptot, cloud area )
416 ! fapoutcld = out-of-cloud precip area
417 ! = max( 0.0, faptot-fapincld )
418 ! also pfoutcld = out-of-cloud precip flux = fraction of total precip flux
419 ! = pf_ls(k)*(fapoutcld/faptot)
420 if (pf_ls(k) > 0.0) then
421 faptot= 0.5*fac_pwght_ls(k)
422 fapincld = min( faptot, cldfra(i,k,j) )
423 fapoutcld = max( 0.0, faptot-fapincld )
424 pfoutcld = pf_ls(k)*(fapoutcld/faptot)
425 if (pfoutcld >= pfsmall_ls) then
426 isprx(i,k,j) = .true.
427 pfx(i,k,j) = pfoutcld
428 fapx(i,k,j) = fapoutcld
429 pfx_inrain(i,k,j) = pf_ls(k)/faptot
430 sumpf = sumpf + pfx(i,k,j)
431 sumfa = sumfa + fapx(i,k,j)
432 end if
433 end if
435 4900 continue ! "k = 1, pver"
437 5900 continue ! "i = 1, ncol"
440 call aerimpactscav (ims,ime,kms,kme,jms,jme,its,ite,kts,kte,jts,jte, num_chem, &
441 ntype_aer, nsize_aer, ncomp_aer, massptr_aer, dens_aer, numptr_aer, &
442 maxd_acomp, maxd_asize,maxd_atype, maxd_aphase, ai_phase, &
443 volumcen_sect, volumlo_sect, volumhi_sect, &
444 waterptr_aer, dens_water_aer, &
445 scavimptblvol, scavimptblnum, &
446 nimptblgrow_mind, nimptblgrow_maxd, dlndg_nimptblgrow, &
447 dtstepc, t_phy, p_phy, pdel, chem, &
448 isprx, fapx, pfx, pfx_inrain, &
449 dqdt, dotend, qsrflx )
451 ! only if not MOZART (--> Neu and Prather used)
452 if ( .NOT. (config_flags%chem_opt == mozart_mosaic_4bin_aq_kpp) ) then
453 call gasrainscav (ims,ime,kms,kme,jms,jme,its,ite,kts,kte,jts,jte, num_chem, &
454 config_flags, &
455 dtstepc, t_phy, p_phy, pdel, chem, &
456 isprx, fapx, pfx, pfx_inrain, &
457 dqdt, dotend, qsrflx )
458 end if
460 ! update chem
462 do n=1,num_chem
463 if(dotend(n))then
464 do 6000 j=jts,jte
465 do 6000 k=kts,kte
466 do 6000 i=its,ite
467 chem(i,k,j,n) = chem(i,k,j,n) + dqdt(i,k,j,n)*dtstepc
468 6000 continue
469 end if
470 end do
472 end subroutine wetscav
476 !===========================================================================
477 !===========================================================================
478 subroutine aerimpactscav (ims,ime,kms,kme,jms,jme,its,ite,kts,kte,jts,jte,num_chem, &
479 ntype_aer, nsize_aer, ncomp_aer, massptr_aer, dens_aer, numptr_aer, &
480 maxd_acomp, maxd_asize,maxd_atype, maxd_aphase, ai_phase, &
481 volumcen_sect, volumlo_sect, volumhi_sect, &
482 waterptr_aer, dens_water_aer, &
483 scavimptblvol, scavimptblnum, nimptblgrow_mind, nimptblgrow_maxd, dlndg_nimptblgrow, &
484 deltat, t, pmid, pdel, chem, &
485 isprx, fapx, pfx, pfx_inrain, &
486 dqdt, dotend, qsrflx )
489 !-----------------------------------------------------------------------
490 !
491 ! Purpose:
492 ! Does below cloud scavenging of aerosols through impaction-interception
493 !
494 ! Removal rates for aersol number, (surface - still to do), and volume
495 ! are computed for each mode using precalculated lookup tables.
496 ! The tables account for variables in wet-dgnum, wet density,
497 ! air temperature, and air pressure.
498 !
499 ! Authors: R. Easter
500 !
501 !-----------------------------------------------------------------------
502 USE module_model_constants, only: g,rhowater, mwdry
503 USE module_state_description, only: param_first_scalar
505 implicit none
507 !-----------------------------------------------------------------------
508 !
509 ! Input arguments
510 !
511 ! abbreviations & acronyms
512 ! TMR = tracer mixing ratio
513 ! l-s = large scale
514 !
515 integer, intent(in) :: num_chem ! number of chemical species
516 integer, intent(in) :: ims,ime ! column dimension
517 integer, intent(in) :: kms,kme ! level dimension
518 integer, intent(in) :: jms,jme ! column dimension
519 integer, intent(in) :: its,ite ! column identifier
520 integer, intent(in) :: kts,kte ! level identifier
521 integer, intent(in) :: jts,jte ! column identifier
522 real, intent(in) :: deltat ! model timestep
524 real, intent(in) :: t(ims:ime,kms:kme,jms:jme) ! temperature
525 real, intent(in) :: pmid(ims:ime,kms:kme,jms:jme) ! pressure at model levels
526 real, intent(in) :: pdel(ims:ime,kms:kme,jms:jme) ! pressure thickness of levels
527 real, intent(in) :: chem(ims:ime,kms:kme,jms:jme,num_chem) ! chem array
529 logical, intent(in) :: isprx(ims:ime,kms:kme,jms:jme) ! true if precip at i,k
530 real, intent(in) :: fapx(ims:ime,kms:kme,jms:jme) ! frac. area for precip
531 real, intent(in) :: pfx(ims:ime,kms:kme,jms:jme) ! grid-avg precip
532 ! flux (kg/m2/s)
533 real, intent(in) :: pfx_inrain(ims:ime,kms:kme,jms:jme) ! in-rain-area precip flux (kg/m2/s)
535 real, intent(inout) :: dqdt(ims:ime,kms:kme,jms:jme,num_chem) ! TMR tendency array
536 logical, intent(inout) :: dotend(num_chem) ! flag for doing scav
537 real, intent(inout) :: qsrflx(ims:ime,jms:jme,num_chem)
538 ! changes to column tracer burdens by wet scavenging over current timestep
539 ! this routine adds on the contribution from aerosol impaction and brownian
540 ! diffusion scavenging by precipitation
541 integer, intent(in) :: maxd_atype, maxd_asize, maxd_acomp, maxd_aphase
542 integer, intent(in) :: ai_phase
543 integer, intent(in) :: ntype_aer
544 integer, intent(in) :: nsize_aer( maxd_atype ), & ! number of size bins
545 ncomp_aer( maxd_atype ), & ! number of chemical components
546 massptr_aer( maxd_acomp, maxd_asize, maxd_atype, maxd_aphase), & ! index for mixing ratio
547 waterptr_aer( maxd_asize, maxd_atype ), & ! index for aerosol water
548 numptr_aer( maxd_asize, maxd_atype, maxd_aphase ) ! index for the number mixing ratio
549 real, intent(in) :: volumcen_sect( maxd_asize, maxd_atype ), & ! single-particle volumes (cm3)
550 volumlo_sect( maxd_asize, maxd_atype ), & ! at section center, lower boundary,
551 volumhi_sect( maxd_asize, maxd_atype ) ! upper boundary
552 real, intent(in) :: dens_aer( maxd_acomp, maxd_atype ), & ! material density (g/cm3)
553 dens_water_aer ! water density (g/m3)
555 real, intent(in) :: dlndg_nimptblgrow
556 integer, intent(in) :: nimptblgrow_mind, nimptblgrow_maxd
557 real, intent(in) :: scavimptblnum(4, nimptblgrow_mind:nimptblgrow_maxd, maxd_asize, maxd_atype), &
558 scavimptblvol(4, nimptblgrow_mind:nimptblgrow_maxd, maxd_asize, maxd_atype)
560 !--------------------------Local Variables------------------------------
562 integer :: i,j ! Work index
563 integer :: jgrow, jp ! Work index
564 integer :: k ! Work index
565 integer :: l, ll, m, n ! Work index
567 logical :: ispr_anywhere
569 real :: dry_mass, dry_volu, dry_volu_1p
570 real :: duma
571 real :: dumfhi, dumflo
572 real :: dumimpactamt0, dumimpactamt3, dumimpactamtsum
573 real :: dumimpactratea0, dumimpactratea3
574 real :: dumimpactrateb0, dumimpactrateb3
575 real :: dumdgratio, dumrate, dumwetdens
576 real :: dumlogdens, dumlogptot, dumlogtemp
577 real :: dumnumb
578 real :: dumscavratenum, dumscavratevol
579 real :: pdel_dt_fac
580 real :: pf_to_prmmh
581 real :: scavimptbl1, scavimptbl2, scavimptbl3, scavimptbl4
582 real :: xgrow
583 real :: wet_mass, wet_volu
584 real, save :: third = 0.33333333
586 !-----------------------------------------------------------------------
587 !
589 ! if (ncol .ne. -987654321) return
591 ! precip rates -- 1.0 kgwtr/m2/s = 1.0e-3 m3wtr/m2/s = 1.0e-3 m/s
592 ! = 1.0 mm/s = 3600 mm/h
594 ispr_anywhere = .false.
596 do 5900 i = its,ite
597 do 5900 j = jts,jte
599 do 4900 k = kte,kts,-1
601 ! skip this level if no precip
602 if ( isprx(i,k,j) ) then
603 ispr_anywhere = .true.
604 else
605 goto 4900
606 end if
609 dumlogtemp = log( t(i,k,j) )
610 ! dumlogptot = log( pressure in dynes/cm2 )
611 dumlogptot = log( 10.0*pmid(i,k,j) )
613 ! compute removal amounts for each aerosol mode
614 do 3900 n=1,ntype_aer
615 do 3900 m=1,nsize_aer(n)
616 dry_volu = 0.
617 dry_mass = 0
618 do ll = 1, ncomp_aer(n)
619 l = massptr_aer(ll,m,n,ai_phase)
620 if (l < param_first_scalar) cycle
621 duma = max( chem(i,k,j,l)*1.0e-6, 0.0 )
622 dry_mass = dry_mass + duma ! g-AP/kg-air
623 dry_volu = dry_volu + duma/dens_aer(ll,n) ! cm3-AP/kg-air
624 end do
625 ! If no aerosol is present at this size and type, there is nothing
626 ! to scavenge so go onto the next size bin. wig, 25-Oct-2005
627 if( dry_volu < 1.0e-26 ) goto 3900
629 ! wet volume
630 l = waterptr_aer(m,n)
631 if (l >= param_first_scalar) then
632 duma = max( chem(i,k,j,l)*1.0e-6, 0.0 )
633 else
634 duma = 0.0
635 end if
636 wet_mass = dry_mass + duma
637 wet_volu = dry_volu + duma/dens_water_aer
639 ! calc dry volume of 1 particle
640 dumnumb = chem(i,k,j,numptr_aer(m,n,ai_phase))
641 if (dry_volu > volumhi_sect(m,n)*dumnumb) then
642 dry_volu_1p = volumhi_sect(m,n)
643 else if (dry_volu < volumlo_sect(m,n)*dumnumb) then
644 dry_volu_1p = volumlo_sect(m,n)
645 else
646 dry_volu_1p = dry_volu/dumnumb
647 end if
648 ! interpolate table values using
649 ! dumdgratio = [(actual-wet-diam)/(base-dry-diam)]
650 ! = [(actual-wet-volu)/(base-dry-volu)]**1/3
651 if (wet_volu <= 1.0e9*dry_volu) then
652 duma = wet_volu/dry_volu
653 else
654 duma = 1.0e9
655 end if
656 dumdgratio = ((dry_volu_1p/volumcen_sect(m,n))*duma)**third
658 ! dumwetdens = wet aerosol density in g/cm3
659 dumwetdens = wet_mass/wet_volu
660 dumlogdens = log( dumwetdens )
661 dumimpactamt3 = 0
662 dumimpactamt0 = 0
663 !
664 ! compute impaction scavenging removal amount for volume
665 !
667 if ((dumdgratio .ge. 0.99) .and. (dumdgratio .le. 1.01)) then
668 scavimptbl1 = scavimptblvol(1,0,m,n)
669 scavimptbl2 = scavimptblvol(2,0,m,n)
670 scavimptbl3 = scavimptblvol(3,0,m,n)
671 scavimptbl4 = scavimptblvol(4,0,m,n)
673 else
674 xgrow = log( dumdgratio ) / dlndg_nimptblgrow
675 jgrow = int( xgrow )
676 if (xgrow .lt. 0.) jgrow = jgrow - 1
677 if (jgrow .lt. nimptblgrow_mind) then
678 jgrow = nimptblgrow_mind
679 xgrow = jgrow
680 else
681 jgrow = min( jgrow, nimptblgrow_maxd-1 )
682 end if
684 dumfhi = xgrow - jgrow
685 dumflo = 1. - dumfhi
687 scavimptbl1 = dumflo*scavimptblvol(1,jgrow,m,n) + &
688 dumfhi*scavimptblvol(1,jgrow+1,m,n)
690 scavimptbl2 = dumflo*scavimptblvol(2,jgrow,m,n) + &
691 dumfhi*scavimptblvol(2,jgrow+1,m,n)
693 scavimptbl3 = dumflo*scavimptblvol(3,jgrow,m,n) + &
694 dumfhi*scavimptblvol(3,jgrow+1,m,n)
696 scavimptbl4 = dumflo*scavimptblvol(4,jgrow,m,n) + &
697 dumfhi*scavimptblvol(4,jgrow+1,m,n)
698 end if
700 ! apply temperature and pressure corrections
701 dumimpactratea3 = exp( scavimptbl1 + scavimptbl2*dumlogtemp + &
702 scavimptbl3*dumlogptot + scavimptbl4*dumlogdens )
704 ! dumimpactratea3 = impaction scav rate (1/h) for precip = 1 mm/h
705 ! dumimpactrateb3 = impaction scav rate (1/s) for precip = pfx_inrain
706 ! (dumimpactratea3/3600) = impaction scav rate (1/s) for precip = 1 mm/h
707 ! (pfx_inrain*3600) = in-rain-area precip rate (mm/h)
708 ! dumimpactrateb3 = (dumimpactratea3/3600) * (pfx_inrain*3600)
709 ! dumimpactamt3 = fraction of aerosol removed from entire grid cell
710 ! in time deltat
711 dumimpactamtsum = 0.0
712 dumimpactrateb3 = dumimpactratea3 * pfx_inrain(i,k,j)
713 dumimpactamt3 = (1. - exp(-deltat*dumimpactrateb3)) * fapx(i,k,j)
714 dumimpactamtsum = dumimpactamtsum + dumimpactamt3
715 dumimpactamt3 = min( dumimpactamtsum, 1.0 )
717 ! diagnostic output
718 ! dump = 10.0*pmid(i,k)
719 ! call calc_1_impact_rate( dgncur_awet(i,k,j,m,n), &
720 ! sigmag_amode(n), dumwetdens, &
721 ! t(i,k), dump, &
722 ! dumscavratenum, dumscavratevol, lun )
723 !
724 ! dumr = -9.99e35
725 ! if (dumscavratevol > 1.0e-35) &
726 ! dumr = (dumimpactratea3/dumscavratevol) - 1.0
727 ! write(lun,9440) nstep, lchnk, i, k, jp, &
728 ! (dumtemp-273.16), dumpress*.001
729 ! write(lun,9442) 'rhowet, dgnwet, dgratio, xgrow', &
730 ! dumwetdens, dgncur_awet(i,k,n), dumdgratio, xgrow
731 ! write(lun,9442) 'exa&approx vol rt, relerr, amt', &
732 ! dumscavratevol, dumimpactratea3, dumr, dumimpactamt3
733 ! write(lun,9442) 'pfx_inrain(1:3) ', &
734 ! (pfx_inrain(jp,i,k), jp=1,3)
735 ! write(lun,9442) 'fapx(1:3) ', &
736 ! (fapx(jp,i,k), jp=1,3)
737 !9440 format( / 'ns,lc,i,k,jp, T(C), p(mb)', i6, 2i4, 2i3, 2f7.1 )
738 !9442 format( a, 4(1pe11.3) )
739 ! end diagnostic output
742 !
743 ! compute impaction scavenging removal amount to number
744 !
745 if (numptr_aer(m,n,ai_phase) < param_first_scalar) then
746 dumimpactamt0 = dumimpactamt3
747 goto 3700
748 end if
750 ! interpolate table values using log of (actual-wet-size)/(base-dry-size)
751 if ((dumdgratio .ge. 0.99) .and. (dumdgratio .le. 1.01)) then
752 scavimptbl1 = scavimptblnum(1,0,m,n)
753 scavimptbl2 = scavimptblnum(2,0,m,n)
754 scavimptbl3 = scavimptblnum(3,0,m,n)
755 scavimptbl4 = scavimptblnum(4,0,m,n)
757 else
758 scavimptbl1 = dumflo*scavimptblnum(1,jgrow,m,n) + &
759 dumfhi*scavimptblnum(1,jgrow+1,m,n)
761 scavimptbl2 = dumflo*scavimptblnum(2,jgrow,m,n) + &
762 dumfhi*scavimptblnum(2,jgrow+1,m,n)
764 scavimptbl3 = dumflo*scavimptblnum(3,jgrow,m,n) + &
765 dumfhi*scavimptblnum(3,jgrow+1,m,n)
767 scavimptbl4 = dumflo*scavimptblnum(4,jgrow,m,n) + &
768 dumfhi*scavimptblnum(4,jgrow+1,m,n)
769 end if
771 ! apply temperature and pressure corrections
772 dumimpactratea0 = exp( scavimptbl1 + scavimptbl2*dumlogtemp + &
773 scavimptbl3*dumlogptot + scavimptbl4*dumlogdens )
775 dumimpactamt0 = 0.0
776 dumimpactrateb0 = dumimpactratea0 * pfx_inrain(i,k,j)
777 dumimpactamt0 = dumimpactamt0 + &
778 (1. - exp( -deltat*dumimpactrateb0 )) * fapx(i,k,j)
779 dumimpactamt0 = min( dumimpactamt0, 1.0 )
781 ! diagnostic output
782 ! dumr = -9.99e35
783 ! if (dumscavratenum > 1.0e-35) &
784 ! dumr = (dumimpactratea0/dumscavratenum) - 1.0
785 ! write(lun,9442) 'exa&approx num rt, relerr, amt', &
786 ! dumscavratenum, dumimpactratea0, dumr, dumimpactamt0
787 ! end diagnostic output
790 3700 continue
792 !
793 ! compute tendencies
794 !
795 pdel_dt_fac = deltat*pdel(i,k,j)/g
796 dumrate = -dumimpactamt3/(deltat*(1.0 + 1.0e-8))
797 dumrate = min(0.0,max(-1.0/deltat,dumrate)) ! make sure -1 <= dumrate*deltat <= 0
798 do ll = 1, ncomp_aer(n)
799 l = massptr_aer(ll,m,n,ai_phase)
800 if (l < param_first_scalar) cycle
801 dqdt(i,k,j,l) = chem(i,k,j,l)*dumrate
802 qsrflx(i,j,l) = qsrflx(i,j,l) + dqdt(i,k,j,l)*pdel_dt_fac ! aerosol mass (ug/m2)
803 end do
804 l = waterptr_aer(m,n)
805 if (l >= param_first_scalar) then
806 dqdt(i,k,j,l) = chem(i,k,j,l)*dumrate
807 qsrflx(i,j,l) = qsrflx(i,j,l) + dqdt(i,k,j,l)*pdel_dt_fac ! aerosol water mass (ug/m2)
808 end if
809 l = numptr_aer(m,n,ai_phase)
810 if (l >= param_first_scalar) then
811 dumrate = -dumimpactamt0/(deltat*(1.0 + 1.0e-8))
812 dqdt(i,k,j,l) = chem(i,k,j,l)*dumrate
813 qsrflx(i,j,l) = qsrflx(i,j,l) + dqdt(i,k,j,l)*pdel_dt_fac ! aerosol number (1/m2)
814 end if
818 3900 continue ! "n = 1, ntot_amode"
820 4900 continue ! "k = 1, pver"
822 5900 continue ! "i = 1, ncol"
825 ! set dotend's
826 if ( ispr_anywhere ) then
827 do n=1,ntype_aer
828 do m=1,nsize_aer(n)
829 do ll = 1, ncomp_aer(n)
830 if (massptr_aer(ll,m,n,ai_phase) >= param_first_scalar) &
831 dotend(massptr_aer(ll,m,n,ai_phase)) = .true.
832 end do
833 if (waterptr_aer(m,n) >= param_first_scalar) &
834 dotend(waterptr_aer(m,n)) = .true.
835 if (numptr_aer(m,n,ai_phase) >= param_first_scalar) &
836 dotend(numptr_aer(m,n,ai_phase)) = .true.
837 end do
838 end do
839 end if
842 return
843 end subroutine aerimpactscav
847 !===========================================================================
848 !===========================================================================
849 subroutine initwet( ntype_aer, nsize_aer, ncomp_aer, massptr_aer, &
850 dens_aer, numptr_aer, maxd_acomp, maxd_asize,maxd_atype, &
851 maxd_aphase, dcen_sect, sigmag_aer, waterptr_aer, dens_water_aer, &
852 scavimptblvol, scavimptblnum, nimptblgrow_mind, nimptblgrow_maxd, &
853 dlndg_nimptblgrow)
854 !-----------------------------------------------------------------------
855 !
856 ! Purpose:
857 ! Computes lookup table for aerosol impaction/interception scavenging rates
858 !
859 ! Authors: R. Easter
860 !
861 !-----------------------------------------------------------------------
862 implicit none
864 integer, intent(in) :: maxd_atype, maxd_asize, maxd_acomp, maxd_aphase
865 integer, intent(in) :: ntype_aer
866 integer, intent(in) :: nsize_aer( maxd_atype ) ! number of size bins
867 integer, intent(in) :: ncomp_aer( maxd_atype ) ! number of chemical components
868 integer, intent(in) :: massptr_aer( maxd_acomp, maxd_asize, maxd_atype, maxd_aphase) ! index for mixing ratio
869 integer, intent(in) :: waterptr_aer( maxd_asize, maxd_atype ) ! index for aerosol water
870 integer, intent(in) :: numptr_aer( maxd_asize, maxd_atype, maxd_aphase ) ! index for the number mixing ratio
871 real, intent(in) :: dens_aer( maxd_acomp, maxd_atype ) ! material density (g/cm3)
872 real, intent(in) :: dens_water_aer ! water density (g/m3)
873 real, intent(in) :: sigmag_aer(maxd_asize, maxd_atype)
874 real, intent(in) :: dcen_sect( maxd_asize, maxd_atype )
875 ! "base" volume-mean diameter (cm) == section center diameter
876 ! scav rates for aerosol number and volume are pre-calculated for a range of volume-mean diameters
877 ! as well as aerosol densities and air pressure and temperature
878 ! for specific time/location/aerosol properties, these rates are interpolated, and
879 ! the size interpolation uses [(actual volume-mean diameter)/(base volume-mean diameter)]
881 real, intent(out) :: dlndg_nimptblgrow
882 integer, intent(in) :: nimptblgrow_mind, nimptblgrow_maxd
883 real, intent(out) :: scavimptblnum(4, nimptblgrow_mind:nimptblgrow_maxd, maxd_asize, maxd_atype)
884 real, intent(out) :: scavimptblvol(4, nimptblgrow_mind:nimptblgrow_maxd, maxd_asize, maxd_atype)
886 ! local variables
887 integer nnfit_maxd
888 parameter (nnfit_maxd=27)
890 integer i, j, m, n, jgrow, jdens, jpress, jtemp, ll, mode, nnfit
891 integer lunerr
893 real dg0, dg0_base, press, rhodryaero, rhowetaero, rmserr, &
894 scavratenum, scavratevol, sigmag, &
895 temp, wetdiaratio, wetvolratio
896 real aafitnum(4), xxfitnum(4,nnfit_maxd), yyfitnum(nnfit_maxd)
897 real aafitvol(4), xxfitvol(4,nnfit_maxd), yyfitvol(nnfit_maxd)
900 lunerr = 6
901 dlndg_nimptblgrow = log( 1.25d00 )
903 do 7900 n=1,ntype_aer
904 do 7900 m=1,nsize_aer(n)
906 sigmag = sigmag_aer(m,n)
907 dg0_base = dcen_sect(m,n)*exp( -1.5*((log(sigmag))**2) )
909 ! for setting up the lookup table, use the dry density of the first
910 ! chemical component of the aerosol type (which currently will be so4)
911 rhodryaero = dens_aer(1,n)
913 do 7800 jgrow = nimptblgrow_mind, nimptblgrow_maxd
915 wetdiaratio = exp( jgrow*dlndg_nimptblgrow )
916 dg0 = dg0_base*wetdiaratio
918 wetvolratio = exp( jgrow*dlndg_nimptblgrow*3. )
919 rhowetaero = 1.0 + (rhodryaero-1.0)/wetvolratio
920 rhowetaero = min( rhowetaero, rhodryaero )
922 !
923 ! compute impaction scavenging rates at 9 temp-press pairs and save
924 !
925 nnfit = 0
927 do 5900 jtemp = -1, 1
928 temp = 273.16 + 25.*jtemp
930 do 5900 jpress = -1, 1
931 press = 0.75e6 + 0.25e6*jpress
933 do 5900 jdens = 0, 2
934 rhowetaero = rhodryaero**(0.5*jdens)
936 call calc_1_impact_rate( &
937 dg0, sigmag, rhowetaero, temp, press, &
938 scavratenum, scavratevol, lunerr )
940 nnfit = nnfit + 1
941 if (nnfit .gt. nnfit_maxd) then
942 call wrf_error_fatal('*** subr. aerosol_impact_setup -- nnfit too big')
943 end if
945 xxfitnum(1,nnfit) = 1.
946 xxfitnum(2,nnfit) = log( temp )
947 xxfitnum(3,nnfit) = log( press )
948 xxfitnum(4,nnfit) = log( rhowetaero )
949 yyfitnum(nnfit) = log( scavratenum )
951 xxfitvol(1,nnfit) = 1.
952 xxfitvol(2,nnfit) = log( temp )
953 xxfitvol(3,nnfit) = log( press )
954 xxfitvol(4,nnfit) = log( rhowetaero )
955 yyfitvol(nnfit) = log( scavratevol )
957 5900 continue
959 !
960 ! do linear regression
961 ! log(scavrate) = a1 + a2*log(temp) + a3*log(press) + a4*log(wetdens)
962 !
963 call mlinft( xxfitnum, yyfitnum, aafitnum, nnfit, 4, 4, rmserr )
964 call mlinft( xxfitvol, yyfitvol, aafitvol, nnfit, 4, 4, rmserr )
966 do i = 1, 4
967 scavimptblnum(i,jgrow,m,n) = aafitnum(i)
968 scavimptblvol(i,jgrow,m,n) = aafitvol(i)
969 end do
972 7800 continue
973 7900 continue
975 return
976 end subroutine initwet
980 !===========================================================================
981 !===========================================================================
982 subroutine calc_1_impact_rate( &
983 dg0, sigmag, rhoaero, temp, press, &
984 scavratenum, scavratevol, lunerr )
985 !
986 ! routine computes a single impaction scavenging rate
987 ! for precipitation rate of 1 mm/h
988 !
989 ! dg0 = geometric mean diameter of aerosol number size distrib. (cm)
990 ! sigmag = geometric standard deviation of size distrib.
991 ! rhoaero = density of aerosol particles (g/cm^3)
992 ! temp = temperature (K)
993 ! press = pressure (dyne/cm^2)
994 ! scavratenum = number scavenging rate (1/h)
995 ! scavratevol = volume or mass scavenging rate (1/h)
996 ! lunerr = logical unit for error message
997 !
998 implicit none
1000 ! subr. parameters
1001 integer lunerr
1002 real dg0, sigmag, rhoaero, temp, press, scavratenum, scavratevol
1004 ! local variables
1005 integer nrainsvmax
1006 parameter (nrainsvmax=50)
1007 real rrainsv(nrainsvmax), xnumrainsv(nrainsvmax),&
1008 vfallrainsv(nrainsvmax)
1010 integer naerosvmax
1011 parameter (naerosvmax=51)
1012 real aaerosv(naerosvmax), &
1013 ynumaerosv(naerosvmax), yvolaerosv(naerosvmax)
1015 integer i, ja, jr, na, nr
1016 real a, aerodiffus, aeromass, ag0, airdynvisc, airkinvisc
1017 real anumsum, avolsum, boltzmann, cair, chi
1018 real d, dr, dum, dumfuchs, dx
1019 real ebrown, eimpact, eintercept, etotal, freepath, gravity
1020 real pi, precip, precipmmhr, precipsum
1021 real r, rainsweepout, reynolds, rhi, rhoair, rhowater, rlo, rnumsum
1022 real scavsumnum, scavsumnumbb
1023 real scavsumvol, scavsumvolbb
1024 real schmidt, sqrtreynolds, sstar, stokes, sx
1025 real taurelax, vfall, vfallstp
1026 real x, xg0, xg3, xhi, xlo, xmuwaterair
1029 rlo = .005
1030 rhi = .250
1031 dr = 0.005
1032 nr = 1 + nint( (rhi-rlo)/dr )
1033 if (nr .gt. nrainsvmax) then
1034 call wrf_error_fatal('*** subr. calc_1_impact_rate -- nr > nrainsvmax')
1035 end if
1037 precipmmhr = 1.0
1038 precip = precipmmhr/36000. ! cm/s
1040 ag0 = dg0/2.
1041 sx = log( sigmag )
1042 xg0 = log( ag0 )
1043 xg3 = xg0 + 3.*sx*sx
1045 xlo = xg3 - 4.*sx
1046 xhi = xg3 + 4.*sx
1047 dx = 0.2*sx
1049 dx = max( 0.2*sx, 0.01 )
1050 xlo = xg3 - max( 4.*sx, 2.*dx )
1051 xhi = xg3 + max( 4.*sx, 2.*dx )
1053 na = 1 + nint( (xhi-xlo)/dx )
1054 if (na .gt. naerosvmax) then
1055 call wrf_error_fatal('*** subr. calc_1_impact_rate -- na > naerosvmax')
1056 end if
1058 ! air molar density
1059 cair = press/(8.31436e7*temp)
1060 ! air mass density
1061 rhoair = 28.966*cair
1062 ! molecular freepath
1063 freepath = 2.8052e-10/cair
1064 ! boltzmann constant
1065 boltzmann = 1.3807e-16
1066 ! water density
1067 rhowater = 1.0
1068 ! gravity
1069 gravity = 980.616
1070 ! air dynamic viscosity
1071 airdynvisc = 1.8325e-4 * (416.16/(temp+120.)) * &
1072 ((temp/296.16)**1.5)
1073 ! air kinemaic viscosity
1074 airkinvisc = airdynvisc/rhoair
1075 ! ratio of water viscosity to air viscosity (from Slinn)
1076 xmuwaterair = 60.0
1078 pi = 3.1415926536
1080 !
1081 ! compute rain drop number concentrations
1082 ! rrainsv = raindrop radius (cm)
1083 ! xnumrainsv = raindrop number concentration (#/cm^3)
1084 ! (number in the bin, not number density)
1085 ! vfallrainsv = fall velocity (cm/s)
1086 !
1087 precipsum = 0.
1088 do i = 1, nr
1089 r = rlo + (i-1)*dr
1090 rrainsv(i) = r
1091 xnumrainsv(i) = exp( -r/2.7e-2 )
1093 d = 2.*r
1094 if (d .le. 0.007) then
1095 vfallstp = 2.88e5 * d**2.
1096 else if (d .le. 0.025) then
1097 vfallstp = 2.8008e4 * d**1.528
1098 else if (d .le. 0.1) then
1099 vfallstp = 4104.9 * d**1.008
1100 else if (d .le. 0.25) then
1101 vfallstp = 1812.1 * d**0.638
1102 else
1103 vfallstp = 1069.8 * d**0.235
1104 end if
1106 vfall = vfallstp * sqrt(1.204e-3/rhoair)
1107 vfallrainsv(i) = vfall
1108 precipsum = precipsum + vfall*(r**3)*xnumrainsv(i)
1109 end do
1110 precipsum = precipsum*pi*1.333333
1112 rnumsum = 0.
1113 do i = 1, nr
1114 xnumrainsv(i) = xnumrainsv(i)*(precip/precipsum)
1115 rnumsum = rnumsum + xnumrainsv(i)
1116 end do
1118 !
1119 ! compute aerosol concentrations
1120 ! aaerosv = particle radius (cm)
1121 ! fnumaerosv = fraction of total number in the bin (--)
1122 ! fvolaerosv = fraction of total volume in the bin (--)
1123 !
1124 anumsum = 0.
1125 avolsum = 0.
1126 do i = 1, na
1127 x = xlo + (i-1)*dx
1128 a = exp( x )
1129 aaerosv(i) = a
1130 dum = (x - xg0)/sx
1131 ynumaerosv(i) = exp( -0.5*dum*dum )
1132 yvolaerosv(i) = ynumaerosv(i)*1.3333*pi*a*a*a
1133 anumsum = anumsum + ynumaerosv(i)
1134 avolsum = avolsum + yvolaerosv(i)
1135 end do
1137 do i = 1, na
1138 ynumaerosv(i) = ynumaerosv(i)/anumsum
1139 yvolaerosv(i) = yvolaerosv(i)/avolsum
1140 end do
1143 !
1144 ! compute scavenging
1145 !
1146 scavsumnum = 0.
1147 scavsumvol = 0.
1148 !
1149 ! outer loop for rain drop radius
1150 !
1151 do 5900 jr = 1, nr
1153 r = rrainsv(jr)
1154 vfall = vfallrainsv(jr)
1156 reynolds = r * vfall / airkinvisc
1157 sqrtreynolds = sqrt( reynolds )
1159 !
1160 ! inner loop for aerosol particle radius
1161 !
1162 scavsumnumbb = 0.
1163 scavsumvolbb = 0.
1165 do 5500 ja = 1, na
1167 a = aaerosv(ja)
1169 chi = a/r
1171 dum = freepath/a
1172 dumfuchs = 1. + 1.246*dum + 0.42*dum*exp(-0.87/dum)
1173 taurelax = 2.*rhoaero*a*a*dumfuchs/(9.*rhoair*airkinvisc)
1175 aeromass = 4.*pi*a*a*a*rhoaero/3.
1176 aerodiffus = boltzmann*temp*taurelax/aeromass
1178 schmidt = airkinvisc/aerodiffus
1179 stokes = vfall*taurelax/r
1181 ebrown = 4.*(1. + 0.4*sqrtreynolds*(schmidt**0.3333333)) / &
1182 (reynolds*schmidt)
1184 dum = (1. + 2.*xmuwaterair*chi) / &
1185 (1. + xmuwaterair/sqrtreynolds)
1186 eintercept = 4.*chi*(chi + dum)
1188 dum = log( 1. + reynolds )
1189 sstar = (1.2 + dum/12.) / (1. + dum)
1190 eimpact = 0.
1191 if (stokes .gt. sstar) then
1192 dum = stokes - sstar
1193 eimpact = (dum/(dum+0.6666667)) ** 1.5
1194 end if
1196 etotal = ebrown + eintercept + eimpact
1197 etotal = min( etotal, 1.0 )
1199 rainsweepout = xnumrainsv(jr)*4.*pi*r*r*vfall
1201 scavsumnumbb = scavsumnumbb + rainsweepout*etotal*ynumaerosv(ja)
1202 scavsumvolbb = scavsumvolbb + rainsweepout*etotal*yvolaerosv(ja)
1204 5500 continue
1206 scavsumnum = scavsumnum + scavsumnumbb
1207 scavsumvol = scavsumvol + scavsumvolbb
1208 5900 continue
1210 scavratenum = scavsumnum*3600.
1211 scavratevol = scavsumvol*3600.
1215 return
1216 end subroutine calc_1_impact_rate
1220 !===========================================================================
1221 !===========================================================================
1222 subroutine gasrainscav (ims,ime,kms,kme,jms,jme,its,ite,kts,kte,jts,jte,num_chem, &
1223 config_flags, &
1224 deltat, t, pmid, pdel, &
1225 chem, &
1226 isprx, fapx, pfx, pfx_inrain, &
1227 dqdt, dotend, qsrflx )
1230 !-----------------------------------------------------------------------
1231 !
1232 ! Purpose:
1233 ! Does below cloud scavenging of gases by rain
1234 !
1235 ! Currently does
1236 ! Irreversible uptake of h2so4 and msa
1237 ! Reactive uptake of so2 and h2o2. This is assumed to be rate limited
1238 ! by the mass transfer to rain (and not by aqueous reaction)
1239 !
1240 ! Authors: R. Easter
1241 !
1242 !-----------------------------------------------------------------------
1243 USE module_model_constants, only: g,rhowater, mwdry
1244 use module_configure, only: grid_config_rec_type, &
1245 param_first_scalar, &
1246 p_so2, p_h2o2, p_sulf,p_h2so4, p_msa, &
1247 p_hno3, p_hcl, p_nh3
1249 implicit none
1251 !-----------------------------------------------------------------------
1252 !
1253 ! Input arguments
1254 !
1255 ! abbreviations & acronyms
1256 ! TMR = tracer mixing ratio
1257 ! l-s = large scale
1258 ! dp-cnv = deep convective
1259 ! sh-cnv = shallow convective
1260 !
1261 integer, intent(in) :: num_chem ! number of chemical species
1262 integer, intent(in) :: ims,ime ! column dimension
1263 integer, intent(in) :: kms,kme ! level dimension
1264 integer, intent(in) :: jms,jme ! column dimension
1265 integer, intent(in) :: its,ite ! column identifier
1266 integer, intent(in) :: kts,kte ! level identifier
1267 integer, intent(in) :: jts,jte ! column identifier
1268 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
1270 real, intent(in) :: deltat ! model timestep
1272 real, intent(in) :: t(ims:ime,kms:kme,jms:jme) ! temperature
1273 real, intent(in) :: pmid(ims:ime,kms:kme,jms:jme) ! pressure at model levels
1274 real, intent(in) :: pdel(ims:ime,kms:kme,jms:jme) ! pressure thickness of levels
1275 real, intent(in) :: chem(ims:ime,kms:kme,jms:jme,num_chem) ! TMR array including chemistry
1277 logical, intent(in) :: isprx(ims:ime,kms:kme,jms:jme) ! true if precip at i,k
1278 real, intent(in) :: fapx(ims:ime,kms:kme,jms:jme) ! frac. area for precip
1279 real, intent(in) :: pfx(ims:ime,kms:kme,jms:jme) ! grid-avg precip
1280 ! flux (kg/m2/s)
1281 real, intent(in) :: pfx_inrain(ims:ime,kms:kme,jms:jme) ! precip flux (kg/m2/s)
1283 real, intent(out) :: dqdt(ims:ime,kms:kme,jms:jme,num_chem) ! TMR tendency array
1284 logical, intent(inout) :: dotend(num_chem) ! flag for doing scav
1285 real, intent(inout) :: qsrflx(ims:ime,jms:jme,num_chem)
1286 ! changes to column tracer burdens by wet scavenging over current timestep
1287 ! this routine adds on the contribution from gas scavenging
1288 ! by mass transfer to rain
1290 !--------------------------Local Variables------------------------------
1292 integer :: i, j, k ! x, y, z work index
1293 integer :: jp ! precip type index
1294 integer :: p1st
1295 logical :: ispr_anywhere
1297 integer, parameter :: ng = 7
1298 integer, parameter :: ig_so2 = 1
1299 integer, parameter :: ig_h2o2 = 2
1300 integer, parameter :: ig_h2so4 = 3
1301 integer, parameter :: ig_msa = 4
1302 integer, parameter :: ig_hno3 = 5
1303 integer, parameter :: ig_hcl = 6
1304 integer, parameter :: ig_nh3 = 7
1305 integer :: ig, lg, lg_ptr(ng)
1307 real :: amtscav(ng), amtscav_sub(ng)
1308 real :: fracgas(ng)
1309 real :: fracscav(ng), fracscav_sub(ng)
1310 real :: deltatinv
1311 real :: dum, dumamt, dumprecipmmh, dumpress, dumtemp
1312 real :: pdel_fac
1313 real :: r_gc(ng)
1314 real :: scavrate_hno3
1315 real :: scavrate(ng), scavrate_factor(ng)
1317 !-----------------------------------------------------------------------
1318 !
1320 ! if (ncol .ne. -987654321) return
1322 ! precip rates -- 1.0 kgwtr/m2/s = 1.0e-3 m3wtr/m2/s = 1.0e-3 m/s
1323 ! = 1.0 mm/s = 3600 mm/h
1325 ispr_anywhere = .false.
1326 deltatinv = 1.0/(deltat*(1.0d0 + 1.0d-15))
1328 p1st = param_first_scalar
1330 lg_ptr(ig_so2 ) = p_so2
1331 lg_ptr(ig_h2o2 ) = p_h2o2
1332 lg_ptr(ig_h2so4) = p_sulf
1333 lg_ptr(ig_h2so4) = p_h2so4
1334 lg_ptr(ig_msa ) = p_msa
1335 lg_ptr(ig_hno3 ) = p_hno3
1336 lg_ptr(ig_hcl ) = p_hcl
1337 lg_ptr(ig_nh3 ) = p_nh3
1339 scavrate_factor(ig_so2 ) = 1.08
1340 scavrate_factor(ig_h2o2 ) = 1.38
1341 scavrate_factor(ig_h2so4) = 0.80
1342 scavrate_factor(ig_msa ) = 0.80
1343 scavrate_factor(ig_hno3 ) = 1.00
1344 scavrate_factor(ig_hcl ) = 1.15
1345 scavrate_factor(ig_nh3 ) = 1.59
1348 do 5900 j = jts,jte
1349 do 5900 i = its,ite
1351 do 4900 k = kts,kte
1353 ! skip this level if no precip
1354 if ( isprx(i,k,j) ) then
1355 ispr_anywhere = .true.
1356 else
1357 goto 4900
1358 end if
1361 ! skip this level if below freezing
1362 dumtemp = t(i,k,j)
1363 if (dumtemp .le. 273.16) goto 4900
1364 dumpress = 10.0*pmid(i,k,j)
1366 do ig = 1, ng
1367 fracscav(ig) = 0.0
1368 fracgas(ig) = 1.0
1369 lg = lg_ptr(ig)
1370 if (lg .ge. p1st) then
1371 r_gc(ig) = max( chem(i,k,j,lg), 0.0 )
1372 ! activate this after gas_aqfrac is added to arguments
1373 ! if (lg .le. numgas_aqfrac) &
1374 ! fracgas(ig) = gas_aqfrac(lg)
1375 else
1376 r_gc(ig) = 0.0
1377 end if
1378 end do
1380 if ( .not. isprx(i,k,j) ) goto 3600
1382 ! precip rate in mm/h over rainy portion of the subarea
1383 dumprecipmmh = pfx_inrain(i,k,j)*3600.0
1385 ! rain scavenging rate for hno3 (power law fit to schwarz and levine,
1386 ! with temperature and pressure adjustments) -- units are (1/s)
1387 scavrate_hno3 = 6.262e-5*(dumprecipmmh**0.7366) &
1388 * ((dumtemp/298.0)**1.12) &
1389 * ((1.013e6/dumpress)**.75)
1391 do ig = 1, ng
1392 scavrate(ig) = scavrate_hno3*scavrate_factor(ig)
1393 fracscav_sub(ig) = (1. - exp(-scavrate(ig)*deltat)) &
1394 *fracgas(ig)*fapx(i,k,j)
1395 amtscav_sub(ig) = r_gc(ig)*min( fracscav_sub(ig), 1.0 )
1396 end do
1398 ! for so2 & h2o2, assume aqueous oxidation is fast, so reactive
1399 ! uptake is limited by the smaller of the two mass transfer rates
1400 dumamt = min( amtscav_sub(ig_so2), amtscav_sub(ig_h2o2) )
1401 fracscav_sub(ig_so2 ) = dumamt/max( r_gc(ig_so2 ), 1.0e-30 )
1402 fracscav_sub(ig_h2o2) = dumamt/max( r_gc(ig_h2o2), 1.0e-30 )
1403 amtscav_sub(ig_so2 ) = r_gc(ig_so2 )*min( fracscav_sub(ig_so2 ), 1.0 )
1404 amtscav_sub(ig_h2o2) = r_gc(ig_h2o2)*min( fracscav_sub(ig_h2o2), 1.0 )
1406 ! for nh3, limit uptake by uptake of all acid gases combined
1407 dumamt = 2.0*amtscav_sub(ig_so2) &
1408 + 2.0*amtscav_sub(ig_h2so4) + amtscav_sub(ig_msa) &
1409 + amtscav_sub(ig_hno3) + amtscav_sub(ig_hcl)
1410 dumamt = min( dumamt, amtscav_sub(ig_nh3) )
1411 fracscav_sub(ig_nh3) = dumamt/max( r_gc(ig_nh3), 1.0e-30 )
1412 amtscav_sub(ig_nh3 ) = r_gc(ig_nh3 )*min( fracscav_sub(ig_nh3 ), 1.0 )
1414 do ig = 1, ng
1415 fracscav(ig) = fracscav(ig) + fracscav_sub(ig)
1416 end do
1418 ! diagnostic output
1419 ! write(lun,9440) nstep, lchnk, i, k, jp, &
1420 ! (dumtemp-273.16), dumpress*.001
1421 ! write(lun,9442) 'pfx, pfx_inrain, fapx ', &
1422 ! pfx(jp,i,k), pfx_inrain(jp,i,k), fapx(jp,i,k)
1423 ! write(lun,9442) 'scavrate_so2, h2o2, msa, h2so4 ', &
1424 ! scavrate(ig_so2), scavrate(ig_h2o2), &
1425 ! scavrate(ig_msa), scavrate(ig_h2so4)
1426 ! write(lun,9442) 'rso2gc, rso2g, rh2o2gc, rh2o2g ', &
1427 ! r_gc(ig_so2), r_gc(ig_so2)*fracgas(ig)so2), &
1428 ! r_gc(ig_h2o2), r_gc(ig_h2o2)*fracgas(ig)h2o2),
1429 ! write(lun,9442) 'amtscav_sub so2, h2o2 ', &
1430 ! amtscav_sub(ig_so2), amtscav_sub(ig_h2o2)
1431 ! write(lun,9442) 'fracscav_sub so2, h2o2, msa, h2so4 ', &
1432 ! fracscav_sub(ig_so2), fracscav_sub(ig_h2o2), &
1433 ! fracscav_sub(ig_msa), fracscav_sub(ig_h2so4)
1434 !9440 format( / 'ns,lc,i,k,jp, T(C), p(mb)', i6, 2i4, 2i3, 2f7.1 )
1435 !9442 format( a, 4(1pe11.3) )
1436 ! end diagnostic output
1439 3600 continue
1441 !
1442 ! compute tendencies
1443 !
1444 pdel_fac = (pdel(i,k,j)/(g*mwdry))
1446 do ig = 1, ng
1447 fracscav(ig) = max(0.0,min(1.0,fracscav(ig))) ! make sure 0 <= fracscav <= 1
1448 amtscav(ig) = fracscav(ig)*r_gc(ig)
1449 lg = lg_ptr(ig)
1450 if (lg .ge. p1st) then
1451 dqdt(i,k,j,lg) = -deltatinv*amtscav(ig)
1452 qsrflx(i,j,lg) = qsrflx(i,j,lg) - pdel_fac*amtscav(ig) ! mmol/m2
1453 end if
1454 end do
1456 4900 continue ! "k = 1, pver"
1458 5900 continue ! "i = 1, ncol"
1461 ! set dotend's
1462 if ( ispr_anywhere ) then
1463 do ig = 1, ng
1464 if (lg_ptr(ig) .ge. p1st) dotend(lg_ptr(ig)) = .true.
1465 end do
1466 end if
1469 return
1470 end subroutine gasrainscav
1474 !===========================================================================
1475 !===========================================================================
1476 subroutine mlinft( x, y, a, n, m, mmaxd, rmserr )
1477 !
1478 ! fits y = a(1)*x(1) + a(2)*x(2) + ... + a(m)*x(m)
1479 !
1480 ! x - array containing x values
1481 ! x(i,k) is parameter i, observation k
1482 ! y - array containing y values
1483 ! y(k) is observation
1484 ! a - array !ontaining the regression coefficients
1485 ! n - number of observations
1486 ! m - number of parameters
1487 ! mmaxd - first dimension of the x array
1488 ! rmserr - root mean square residual
1489 ! rmserr = sqrt( avg-sq-err )
1490 ! avg-sq-err = (sum of residuals squared)/(number of values)
1491 ! residual = y - (a1*x1 + a2*x2 + ... + am*xm)
1492 !
1493 implicit none
1495 ! subr. parameters
1496 integer n, m, mmaxd
1497 real x(mmaxd,n), y(n), a(mmaxd), rmserr
1499 ! local variables
1500 integer i, j, jflag, k
1501 real aa(10,10), bb(10), errsq, resid, ydum
1503 if (n .le. 1) then
1504 a(1) = 1.e30
1505 rmserr = 0.
1506 return
1507 end if
1509 do 2900 i = 1, m
1510 do 2100 j = 1, m
1511 aa(i,j) = 0.0
1512 2100 continue
1513 bb(i) = 0.0
1515 do 2500 k = 1, n
1516 do 2300 j = 1, m
1517 aa(i,j) = aa(i,j) + x(i,k)*x(j,k)
1518 2300 continue
1519 bb(i) = bb(i) + x(i,k)*y(k)
1520 2500 continue
1522 2900 continue
1524 ! do 4100 i = 1, m
1525 ! write(13,9300) i, (aa(i,j), j=1,m), bb(i)
1526 !4100 continue
1527 !9300 format( i5, 5f15.2 )
1529 ! subr linsolv( a, x, b, n, m1, m2, jflag )
1530 call linsolv( aa, a, bb, m, 10, 10, jflag )
1533 errsq = 0.
1534 do 3300 k = 1, n
1535 ydum = 0.0
1536 do 3100 i = 1, m
1537 ydum = ydum + a(i)*x(i,k)
1538 3100 continue
1539 resid = ydum - y(k)
1540 errsq = errsq + resid*resid
1541 3300 continue
1542 rmserr = sqrt( errsq/n )
1544 return
1545 end subroutine mlinft
1549 !===========================================================================
1550 !===========================================================================
1551 subroutine linsolv( a, x, b, n, m1, m2, jflag )
1552 !
1553 ! solves linear eqn system a*x = b using gaussian-elimination
1554 ! with partial pivoting
1555 !
1556 ! n = order of the system
1557 ! m1, m2 = fortran dimensions of a array
1558 ! jflag = completion flag
1559 ! 1 - system solved successfully
1560 ! 0 - system is singular or close to it, and could not be solved.
1561 ! computation was halted to avoid overflow or divide by zero.
1562 !
1563 ! *** note *** rsmall should be defined as close to but somewhat larger
1564 ! than the smallest single precision real on the computer.
1565 !
1566 ! initial coding on 29-aug-86 by r.c. easter
1567 ! change on 4-feb-89 by r.c.easter -- added jflag to parameter list
1568 ! and checking for singularity
1569 !
1570 implicit none
1572 ! subr. parameters
1573 integer n, m1, m2, jflag
1574 real a(m1,m2), x(n), b(n)
1576 ! local variables
1577 integer i, imax, iup, j, k
1578 real amax, asmall, dmy, rsmall
1579 parameter (rsmall = 1.0e-16)
1581 jflag = 0
1583 !
1584 ! reduce coef. matrix to upper triangular form
1585 !
1586 do 1900 k = 1, n
1587 !
1588 ! find pivot element, and
1589 ! move pivot row into row k if necessary
1590 !
1591 imax = k
1592 amax = abs( a(imax,k) )
1593 do 1200 i = k+1, n
1594 if (abs(a(i,k)) .gt. amax) then
1595 imax = i
1596 amax = abs(a(i,k))
1597 end if
1598 1200 continue
1599 if (amax .eq. 0.) return
1601 if (imax .ne. k) then
1602 do 1400 j = k, n
1603 dmy = a(imax,j)
1604 a(imax,j) = a(k,j)
1605 a(k,j) = dmy
1606 1400 continue
1607 dmy = b(imax)
1608 b(imax) = b(k)
1609 b(k) = dmy
1610 end if
1612 !
1613 ! reduce
1614 !
1615 asmall = abs(a(k,k))
1616 do 1700 i = k+1, n
1617 if (a(i,k) .ne. 0.0) then
1618 if (asmall .le. abs(rsmall*a(i,k))) return
1619 dmy = a(i,k)/a(k,k)
1620 a(i,k) = 0.0
1621 do 1600 j = k+1, n
1622 a(i,j) = a(i,j) - dmy*a(k,j)
1623 1600 continue
1624 b(i) = b(i) - dmy*b(k)
1625 end if
1626 1700 continue
1628 1900 continue
1630 !
1631 ! backsolve
1632 !
1633 do 2900 iup = 1, n
1634 i = n + 1 - iup
1635 dmy = b(i)
1636 do 2500 j = i+1, n
1637 dmy = dmy - a(i,j)*x(j)
1638 2500 continue
1639 asmall = abs(a(i,i))
1640 if (abs(a(i,i)) .le. abs(rsmall*dmy)) return
1641 x(i) = dmy/a(i,i)
1642 2900 continue
1644 jflag = 1
1646 return
1647 end subroutine linsolv
1649 END MODULE module_mosaic_wetscav