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