Line data Source code
1 : !===============================================================================
2 : ! CAMRA Aerosol Model
3 : !===============================================================================
4 : module aero_model
5 : use physics_buffer, only: physics_buffer_desc, pbuf_get_index, pbuf_add_field, dtype_r8
6 : use shr_kind_mod, only: r8 => shr_kind_r8
7 : use constituents, only: pcnst, cnst_name, cnst_get_ind
8 : use perf_mod, only: t_startf, t_stopf
9 : use ppgrid, only: pcols, pver, pverp
10 : use phys_control, only: phys_getopts, cam_physpkg_is
11 : use cam_abortutils, only: endrun
12 : use cam_logfile, only: iulog
13 : use physics_types, only: physics_state, physics_ptend, physics_ptend_init
14 : use camsrfexch, only: cam_in_t, cam_out_t
15 : use physics_buffer, only: pbuf_get_field, pbuf_set_field, dtype_r8
16 : use physconst, only: gravit, rair, rhoh2o
17 : use spmd_utils, only: masterproc
18 : use cam_history, only: outfld
19 : use chem_mods, only: gas_pcnst, adv_mass
20 : use mo_tracname, only: solsym
21 : use infnan, only: nan, assignment(=)
22 : use rad_constituents, only: rad_cnst_get_info, rad_cnst_get_info_by_bin, &
23 : rad_cnst_get_info_by_bin_spec, rad_cnst_get_bin_props_by_idx, &
24 : rad_cnst_get_bin_mmr_by_idx
25 : use mo_setsox, only: setsox, has_sox
26 : use carma_aerosol_properties_mod, only: carma_aerosol_properties
27 :
28 : use carma_intr, only: carma_get_group_by_name, carma_get_dry_radius, carma_get_wet_radius, carma_get_bin_rmass
29 : use carma_intr, only: carma_get_sad
30 :
31 : use aerosol_properties_mod, only: aero_name_len
32 :
33 : implicit none
34 : private
35 :
36 : public :: aero_model_readnl
37 : public :: aero_model_register
38 : public :: aero_model_init
39 : public :: aero_model_gasaerexch ! create, grow, change, and shrink aerosols.
40 : public :: aero_model_drydep ! aerosol dry deposition and sediment
41 : public :: aero_model_wetdep ! aerosol wet removal
42 : public :: aero_model_emissions ! aerosol emissions
43 : public :: aero_model_surfarea ! tropospheric aerosol wet surface area for chemistry
44 : public :: aero_model_strat_surfarea ! stub
45 :
46 : ! Misc private data
47 : character(len=32), allocatable :: fieldname(:) ! names for interstitial output fields
48 : character(len=32), allocatable :: fieldname_cw(:) ! names for cloud_borne output fields
49 :
50 : integer :: fracis_idx = 0
51 : integer :: prain_idx = 0
52 : integer :: rprddp_idx = 0
53 : integer :: rprdsh_idx = 0
54 : integer :: nevapr_shcu_idx = 0
55 : integer :: nevapr_dpcu_idx = 0
56 : integer :: nh3_ndx = 0
57 : integer :: nh4_ndx = 0
58 : integer :: h2so4_ndx = 0
59 :
60 : ! variables for table lookup of aerosol impaction/interception scavenging rates
61 : integer, parameter :: nimptblgrow_mind=-7, nimptblgrow_maxd=12
62 :
63 :
64 : ! description of bin aerosols
65 : integer, public, protected :: nspec_max = 0
66 : integer, public, protected :: nbins = 0
67 : integer, public, protected, allocatable :: nspec(:)
68 :
69 : ! local indexing for bins
70 : integer, allocatable :: bin_idx(:,:) ! table for local indexing of modal aero number and mmr
71 : integer :: ncnst_tot ! total number of mode number conc + mode species
72 : integer :: ncnst_extd ! twiece total number of mode number conc + mode species
73 :
74 : ! Indices for CARMA species in the ptend%q array. Needed for prognostic aerosol case.
75 : logical, allocatable :: bin_cnst_lq(:,:)
76 : integer, allocatable :: bin_cnst_idx(:,:)
77 :
78 :
79 : ! ptr2d_t is used to create arrays of pointers to 2D fields
80 : type ptr2d_t
81 : real(r8), pointer :: fld(:,:) => null()
82 : end type ptr2d_t
83 :
84 : logical :: lq(pcnst) = .false. ! set flags true for constituents with non-zero tendencies
85 : ! in the ptend object
86 :
87 : ! Namelist variables
88 : real(r8) :: sol_facti_cloud_borne = 1._r8
89 : real(r8) :: sol_factb_interstitial = 0.1_r8
90 : real(r8) :: sol_factic_interstitial = 0.4_r8
91 :
92 : logical :: convproc_do_aer
93 :
94 : type(carma_aerosol_properties), pointer :: aero_props =>null()
95 :
96 : contains
97 :
98 : !=============================================================================
99 : ! reads aerosol namelist options
100 : !=============================================================================
101 1536 : subroutine aero_model_readnl(nlfile)
102 :
103 : use namelist_utils, only: find_group_name
104 : use units, only: getunit, freeunit
105 : use aero_wetdep_cam, only: aero_wetdep_readnl
106 : use mpishorthand
107 :
108 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
109 :
110 : ! Local variables
111 : integer :: unitn, ierr
112 : character(len=*), parameter :: subname = 'aero_model_readnl'
113 :
114 : ! Namelist variables
115 : namelist /aerosol_nl/ sol_facti_cloud_borne, sol_factb_interstitial, sol_factic_interstitial
116 :
117 : !-----------------------------------------------------------------------------
118 :
119 : ! Read namelist
120 1536 : if (masterproc) then
121 2 : unitn = getunit()
122 2 : open( unitn, file=trim(nlfile), status='old' )
123 2 : call find_group_name(unitn, 'aerosol_nl', status=ierr)
124 2 : if (ierr == 0) then
125 0 : read(unitn, aerosol_nl, iostat=ierr)
126 0 : if (ierr /= 0) then
127 0 : call endrun(subname // ':: ERROR reading namelist')
128 : end if
129 : end if
130 2 : close(unitn)
131 2 : call freeunit(unitn)
132 : end if
133 :
134 : #ifdef SPMD
135 : ! Broadcast namelist variables
136 1536 : call mpibcast(sol_facti_cloud_borne, 1, mpir8, 0, mpicom)
137 1536 : call mpibcast(sol_factb_interstitial, 1, mpir8, 0, mpicom)
138 1536 : call mpibcast(sol_factic_interstitial, 1, mpir8, 0, mpicom)
139 : #endif
140 :
141 1536 : call aero_wetdep_readnl(nlfile)
142 :
143 1536 : end subroutine aero_model_readnl
144 :
145 : !=============================================================================
146 : !=============================================================================
147 1536 : subroutine aero_model_register()
148 :
149 1536 : use carma_flags_mod, only: carma_model
150 :
151 : integer :: m, l, i
152 : integer :: nsoa_vbs
153 : character(len=32) :: num_name
154 : character(len=32) :: num_name_cw
155 : character(len=32) :: spec_name_cw
156 :
157 : integer :: idx, ierr
158 :
159 1536 : call rad_cnst_get_info( 0, nbins=nbins)
160 4608 : allocate( nspec(nbins), stat=ierr )
161 1536 : if (ierr/=0) call endrun('aero_model_register: allocate error')
162 :
163 : ! add pbuf fields for interstitial (cloud borne) aerosols in CARMA
164 62976 : do m = 1, nbins
165 61440 : call rad_cnst_get_info_by_bin(0, m, num_name=num_name, num_name_cw=num_name_cw, nspec=nspec(m))
166 61440 : call pbuf_add_field(num_name,'global',dtype_r8,(/pcols,pver/), idx)
167 61440 : call pbuf_add_field(num_name_cw,'global',dtype_r8,(/pcols,pver/), idx)
168 278016 : do l = 1, nspec(m)
169 215040 : call rad_cnst_get_info_by_bin_spec(0, m, l, spec_name_cw=spec_name_cw)
170 276480 : call pbuf_add_field(spec_name_cw,'global',dtype_r8,(/pcols,pver/),idx)
171 : enddo
172 : enddo
173 :
174 : ! SOA information
175 : ! Define number of VBS bins (nsoa) based on number of SOAG chemistry species
176 1536 : nsoa_vbs = 0
177 336384 : do i = 1, pcnst
178 336384 : if (cnst_name(i)(:4) == 'SOAG') then
179 1536 : nsoa_vbs = nsoa_vbs + 1
180 : end if
181 : end do
182 1536 : if (masterproc) then
183 2 : write(iulog,*) 'nsoa_vbs = ', nsoa_vbs
184 : endif
185 :
186 : ! Define pbuf field for soa_fraction
187 7680 : call pbuf_add_field('FRACVBS','global',dtype_r8,(/pcols,pver,nbins,nsoa_vbs/), idx)
188 :
189 1536 : end subroutine aero_model_register
190 :
191 : !=============================================================================
192 : !=============================================================================
193 1536 : subroutine aero_model_init( pbuf2d )
194 :
195 : use mo_chem_utls, only: get_inv_ndx
196 : use cam_history, only: addfld, add_default, horiz_only
197 : use mo_chem_utls, only: get_rxt_ndx, get_spc_ndx
198 : use aero_wetdep_cam, only: aero_wetdep_init
199 : use mo_setsox, only: sox_inti
200 : use carma_aero_gasaerexch, only: carma_aero_gasaerexch_init
201 :
202 : use time_manager, only: is_first_step
203 : use constituents, only: cnst_set_convtran2
204 : use aero_deposition_cam, only: aero_deposition_cam_init
205 :
206 : ! args
207 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
208 :
209 :
210 : ! local vars
211 : character(len=*), parameter :: subrname = 'aero_model_init'
212 : integer :: m, n, ii, mm
213 : integer :: idxtmp = -1
214 :
215 : logical :: history_aerosol ! Output MAM or SECT aerosol tendencies
216 : logical :: history_chemistry, history_cesm_forcing
217 :
218 : integer :: l
219 :
220 : character(len=2) :: unit_basename ! Units 'kg' or '1'
221 : character(len=32) :: num_name
222 : character(len=32) :: num_name_cw
223 : character(len=32) :: spec_name_cw
224 :
225 : integer :: idx, ierr
226 : real(r8) :: nanval
227 :
228 1536 : aero_props => carma_aerosol_properties()
229 1536 : call aero_deposition_cam_init(aero_props)
230 :
231 1536 : if (is_first_step()) then
232 31488 : do m = 1, nbins
233 30720 : call rad_cnst_get_info_by_bin(0, m, num_name=num_name, num_name_cw=num_name_cw)
234 30720 : idx = pbuf_get_index(num_name)
235 30720 : call pbuf_set_field(pbuf2d, idx, 0.0_r8)
236 30720 : idx = pbuf_get_index(num_name_cw)
237 30720 : call pbuf_set_field(pbuf2d, idx, 0.0_r8)
238 139008 : do l = 1, nspec(m)
239 107520 : call rad_cnst_get_info_by_bin_spec(0, m, l, spec_name_cw=spec_name_cw)
240 107520 : idx = pbuf_get_index(spec_name_cw)
241 138240 : call pbuf_set_field(pbuf2d, idx, 0.0_r8)
242 : enddo
243 : enddo
244 : endif
245 :
246 : ! define pbuf field for soa_fraction
247 1536 : if (is_first_step()) then
248 768 : nanval = nan
249 768 : idx = pbuf_get_index('FRACVBS')
250 768 : call pbuf_set_field(pbuf2d, idx, nanval)
251 : end if
252 :
253 : ! aqueous chem initialization
254 1536 : call sox_inti()
255 :
256 1536 : h2so4_ndx = get_spc_ndx('H2SO4')
257 1536 : nh3_ndx = get_spc_ndx('NH3')
258 1536 : nh4_ndx = get_spc_ndx('NH4')
259 :
260 1536 : fracis_idx = pbuf_get_index('FRACIS')
261 1536 : prain_idx = pbuf_get_index('PRAIN')
262 1536 : rprddp_idx = pbuf_get_index('RPRDDP')
263 1536 : rprdsh_idx = pbuf_get_index('RPRDSH')
264 1536 : nevapr_shcu_idx = pbuf_get_index('NEVAPR_SHCU')
265 1536 : nevapr_dpcu_idx = pbuf_get_index('NEVAPR_DPCU')
266 :
267 : call phys_getopts(history_aerosol_out = history_aerosol, &
268 : history_chemistry_out=history_chemistry, &
269 : history_cesm_forcing_out=history_cesm_forcing, &
270 1536 : convproc_do_aer_out = convproc_do_aer)
271 :
272 1536 : call carma_aero_gasaerexch_init
273 :
274 62976 : nspec_max = maxval(nspec)
275 :
276 1536 : ncnst_tot = nspec(1)
277 61440 : do m = 2, nbins
278 61440 : ncnst_tot = ncnst_tot + nspec(m)
279 : end do
280 1536 : ncnst_extd = 2*ncnst_tot
281 :
282 : allocate( &
283 : bin_idx(nbins,nspec_max), &
284 : bin_cnst_lq(nbins,nspec_max), &
285 : bin_cnst_idx(nbins,nspec_max), &
286 : fieldname_cw(ncnst_tot), &
287 16896 : fieldname(ncnst_tot), stat=ierr )
288 1536 : if (ierr/=0) call endrun(subrname//' : allocate error')
289 :
290 1536 : ii = 0
291 62976 : do m = 1, nbins
292 278016 : do l = 1, nspec(m) ! loop through species
293 215040 : ii = ii + 1
294 215040 : bin_idx(m,l) = ii
295 :
296 215040 : if (l <= nspec(m) ) then ! species
297 215040 : call rad_cnst_get_info_by_bin_spec(0, m, l, spec_name=fieldname(ii), spec_name_cw=fieldname_cw(ii))
298 : else !number
299 0 : call rad_cnst_get_info_by_bin(0, m, num_name=fieldname(ii), num_name_cw=fieldname_cw(ii))
300 : end if
301 :
302 215040 : call cnst_get_ind(fieldname(ii), idxtmp, abort=.false.)
303 215040 : if (idxtmp.gt.0) then
304 215040 : bin_cnst_lq(m,l) = .true.
305 215040 : bin_cnst_idx(m,l) = idxtmp
306 215040 : lq(idxtmp) = .true.
307 215040 : call cnst_set_convtran2(idxtmp, .not.convproc_do_aer)
308 : else
309 0 : bin_cnst_lq(m,l) = .false.
310 0 : bin_cnst_idx(m,l) = 0
311 : end if
312 :
313 215040 : mm = ii
314 :
315 215040 : unit_basename = 'kg'
316 215040 : if (l == nspec(m) + 2) then ! number
317 0 : unit_basename = ' 1'
318 : end if
319 :
320 :
321 0 : call addfld( fieldname_cw(mm), (/ 'lev' /), 'A', unit_basename//'/kg ', &
322 430080 : trim(fieldname_cw(mm))//' in cloud water')
323 0 : call addfld (trim(fieldname_cw(mm))//'DDF', horiz_only, 'A', unit_basename//'/m2/s ', &
324 215040 : trim(fieldname_cw(mm))//' dry deposition flux at bottom (grav + turb)')
325 0 : call addfld (trim(fieldname_cw(mm))//'TBF', horiz_only, 'A', unit_basename//'/m2/s ', &
326 215040 : trim(fieldname_cw(mm))//' turbulent dry deposition flux')
327 0 : call addfld (trim(fieldname_cw(mm))//'GVF', horiz_only, 'A', unit_basename//'/m2/s ', &
328 215040 : trim(fieldname_cw(mm))//' gravitational dry deposition flux')
329 :
330 215040 : if ( history_aerosol.or. history_chemistry ) then
331 0 : call add_default( fieldname_cw(mm), 1, ' ' )
332 : endif
333 276480 : if ( history_aerosol ) then
334 0 : call add_default (trim(fieldname_cw(mm))//'GVF', 1, ' ')
335 0 : call add_default (trim(fieldname_cw(mm))//'TBF', 1, ' ')
336 0 : call add_default (trim(fieldname_cw(mm))//'DDF', 1, ' ')
337 : endif
338 : enddo
339 : enddo
340 :
341 125952 : do m = 1,gas_pcnst
342 :
343 124416 : unit_basename = 'kg' ! Units 'kg' or '1'
344 :
345 124416 : call addfld( 'GS_'//trim(solsym(m)),horiz_only, 'A', unit_basename//'/m2/s ', &
346 248832 : trim(solsym(m))//' gas chemistry/wet removal (for gas species)')
347 124416 : call addfld( 'AQ_'//trim(solsym(m)),horiz_only, 'A', unit_basename//'/m2/s ', &
348 248832 : trim(solsym(m))//' aqueous chemistry (for gas species)')
349 125952 : if ( history_aerosol ) then
350 0 : call add_default( 'AQ_'//trim(solsym(m)), 1, ' ')
351 : endif
352 :
353 : enddo
354 :
355 1536 : if (has_sox) then
356 62976 : do n = 1, nbins
357 278016 : do l = 1, nspec(n) ! not for total mass or number
358 215040 : mm = bin_idx(n, l)
359 : call addfld (&
360 0 : trim(fieldname_cw(mm))//'AQSO4',horiz_only, 'A','kg/m2/s', &
361 215040 : trim(fieldname_cw(mm))//' aqueous phase chemistry')
362 : call addfld (&
363 0 : trim(fieldname_cw(mm))//'AQH2SO4',horiz_only, 'A','kg/m2/s', &
364 215040 : trim(fieldname_cw(mm))//' aqueous phase chemistry')
365 276480 : if ( history_aerosol ) then
366 0 : call add_default (trim(fieldname_cw(mm))//'AQSO4', 1, ' ')
367 0 : call add_default (trim(fieldname_cw(mm))//'AQH2SO4', 1, ' ')
368 : endif
369 : end do
370 : end do
371 :
372 3072 : call addfld( 'XPH_LWC', (/ 'lev' /), 'A','kg/kg', 'pH value multiplied by lwc')
373 1536 : call addfld ('AQSO4_H2O2', horiz_only, 'A','kg/m2/s', 'SO4 aqueous phase chemistry due to H2O2')
374 1536 : call addfld ('AQSO4_O3', horiz_only, 'A','kg/m2/s', 'SO4 aqueous phase chemistry due to O3')
375 :
376 1536 : if ( history_aerosol ) then
377 0 : call add_default ('XPH_LWC', 1, ' ')
378 0 : call add_default ('AQSO4_H2O2', 1, ' ')
379 0 : call add_default ('AQSO4_O3', 1, ' ')
380 : endif
381 : endif
382 :
383 1536 : call aero_wetdep_init()
384 :
385 3072 : end subroutine aero_model_init
386 :
387 : !=============================================================================
388 : !=============================================================================
389 15978240 : subroutine aero_model_drydep ( state, pbuf, obklen, ustar, cam_in, dt, cam_out, ptend )
390 :
391 : ! args
392 : type(physics_state), intent(in) :: state ! Physics state variables
393 : real(r8), intent(in) :: obklen(:)
394 : real(r8), intent(in) :: ustar(:) ! sfc fric vel
395 : type(cam_in_t), target, intent(in) :: cam_in ! import state
396 : real(r8), intent(in) :: dt ! time step
397 : type(cam_out_t), intent(inout) :: cam_out ! export state
398 : type(physics_ptend), intent(out) :: ptend ! indivdual parameterization tendencies
399 : type(physics_buffer_desc), pointer :: pbuf(:)
400 :
401 1536 : endsubroutine aero_model_drydep
402 :
403 : !=============================================================================
404 : !=============================================================================
405 17660160 : subroutine aero_model_wetdep( state, dt, dlf, cam_out, ptend, pbuf)
406 :
407 : use aero_wetdep_cam, only: aero_wetdep_tend
408 :
409 : ! args
410 :
411 : type(physics_state), intent(in) :: state ! Physics state variables
412 : real(r8), intent(in) :: dt ! time step
413 : real(r8), intent(in) :: dlf(:,:) ! shallow+deep convective detrainment [kg/kg/s]
414 : type(cam_out_t), intent(inout) :: cam_out ! export state
415 : type(physics_ptend), intent(out) :: ptend ! indivdual parameterization tendencies
416 : type(physics_buffer_desc), pointer :: pbuf(:)
417 :
418 80640 : call aero_wetdep_tend(state, dt, dlf, cam_out, ptend, pbuf)
419 :
420 80640 : end subroutine aero_model_wetdep
421 :
422 : !-------------------------------------------------------------------------
423 : ! provides wet tropospheric aerosol surface area info for sectional aerosols
424 : ! called from mo_usrrxt
425 : !-------------------------------------------------------------------------
426 72960 : subroutine aero_model_surfarea( &
427 72960 : state, mmr, radmean, relhum, pmid, temp, strato_sad, sulfate, m, ltrop, &
428 72960 : dlat, het1_ndx, pbuf, ncol, sfc, dm_aer, sad_trop, reff_trop )
429 :
430 : ! dummy args
431 : type(physics_state), intent(in) :: state ! Physics state variables
432 : real(r8), intent(in) :: pmid(:,:)
433 : real(r8), intent(in) :: temp(:,:)
434 : real(r8), intent(in) :: mmr(:,:,:)
435 : real(r8), intent(in) :: radmean ! mean radii in cm
436 : real(r8), intent(in) :: strato_sad(:,:)
437 : integer, intent(in) :: ncol
438 : integer, intent(in) :: ltrop(:)
439 : real(r8), intent(in) :: dlat(:) ! degrees latitude
440 : integer, intent(in) :: het1_ndx
441 : real(r8), intent(in) :: relhum(:,:)
442 : real(r8), intent(in) :: m(:,:) ! total atm density (/cm^3)
443 : real(r8), intent(in) :: sulfate(:,:)
444 : type(physics_buffer_desc), pointer :: pbuf(:)
445 :
446 : real(r8), intent(inout) :: sfc(:,:,:)
447 : real(r8), intent(inout) :: dm_aer(:,:,:)
448 : real(r8), intent(inout) :: sad_trop(:,:) ! aerosol surface area density (cm2/cm3), zeroed above the tropopause
449 : real(r8), intent(out) :: reff_trop(:,:) ! aerosol effective radius (cm), zeroed above the tropopause
450 :
451 : ! local vars
452 145920 : integer :: beglev(ncol)
453 145920 : integer :: endlev(ncol)
454 :
455 1196544 : beglev(:ncol)=ltrop(:ncol)+1
456 1123584 : endlev(:ncol)=pver
457 72960 : call surf_area_dens( state, pbuf, ncol, mmr, beglev, endlev, sad_trop, reff_trop, sfc=sfc, dm_aer=dm_aer )
458 :
459 80640 : end subroutine aero_model_surfarea
460 :
461 : !-------------------------------------------------------------------------
462 : ! provides wet stratospheric aerosol surface area info for sectional aerosols
463 : ! called from mo_gas_phase_chemdr.F90
464 : !-------------------------------------------------------------------------
465 72960 : subroutine aero_model_strat_surfarea( state, ncol, mmr, pmid, temp, ltrop, pbuf, strato_sad, reff_strat )
466 :
467 : use ref_pres, only: clim_modal_aero_top_lev
468 :
469 : ! dummy args
470 : type(physics_state), intent(in) :: state ! Physics state variables
471 : integer, intent(in) :: ncol
472 : real(r8), intent(in) :: mmr(:,:,:)
473 : real(r8), intent(in) :: pmid(:,:)
474 : real(r8), intent(in) :: temp(:,:)
475 : integer, intent(in) :: ltrop(:) ! tropopause level indices
476 : type(physics_buffer_desc), pointer :: pbuf(:)
477 : real(r8), intent(out) :: strato_sad(:,:) ! aerosol surface area density (cm2/cm3), zeroed below the tropopause
478 : real(r8), intent(out) :: reff_strat(:,:) ! aerosol effective radius (cm), zeroed below the tropopause
479 :
480 : ! local vars
481 145920 : integer :: beglev(ncol)
482 145920 : integer :: endlev(ncol)
483 :
484 1123584 : beglev(:ncol) = clim_modal_aero_top_lev
485 1196544 : endlev(:ncol) = ltrop(:ncol)
486 :
487 72960 : call surf_area_dens( state, pbuf, ncol, mmr, beglev, endlev, strato_sad, reff_strat )
488 :
489 72960 : end subroutine aero_model_strat_surfarea
490 :
491 : !=============================================================================
492 : !=============================================================================
493 218880 : subroutine aero_model_gasaerexch( state, loffset, ncol, lchnk, troplev, delt, reaction_rates, &
494 72960 : tfld, pmid, pdel, mbar, relhum, &
495 145920 : zm, qh2o, cwat, cldfr, cldnum, &
496 72960 : airdens, invariants, del_h2so4_gasprod, &
497 72960 : vmr0, vmr, pbuf )
498 :
499 72960 : use carma_aero_gasaerexch, only : carma_aero_gasaerexch_sub
500 : use time_manager, only : get_nstep
501 : !-----------------------------------------------------------------------
502 : ! ... dummy arguments
503 : !-----------------------------------------------------------------------
504 : type(physics_state), intent(in) :: state ! Physics state variables
505 : integer, intent(in) :: loffset ! offset applied to modal aero "pointers"
506 : integer, intent(in) :: ncol ! number columns in chunk
507 : integer, intent(in) :: lchnk ! chunk index
508 : integer, intent(in) :: troplev(:)
509 : real(r8), intent(in) :: delt ! time step size (sec)
510 : real(r8), intent(in) :: reaction_rates(:,:,:) ! reaction rates
511 : real(r8), intent(in) :: tfld(:,:) ! temperature (K)
512 : real(r8), intent(in) :: pmid(:,:) ! pressure at model levels (Pa)
513 : real(r8), intent(in) :: pdel(:,:) ! pressure thickness of levels (Pa)
514 : real(r8), intent(in) :: mbar(:,:) ! mean wet atmospheric mass ( amu )
515 : real(r8), intent(in) :: relhum(:,:) ! relative humidity
516 : real(r8), intent(in) :: airdens(:,:) ! total atms density (molec/cm**3)
517 : real(r8), intent(in) :: invariants(:,:,:)
518 : real(r8), intent(in) :: del_h2so4_gasprod(:,:)
519 : real(r8), intent(in) :: zm(:,:)
520 : real(r8), intent(in) :: qh2o(:,:)
521 : real(r8), intent(in) :: cwat(:,:) ! cloud liquid water content (kg/kg)
522 : real(r8), intent(in) :: cldfr(:,:)
523 : real(r8), intent(in) :: cldnum(:,:) ! droplet number concentration (#/kg)
524 : real(r8), intent(in) :: vmr0(:,:,:) ! initial mixing ratios (before gas-phase chem changes)
525 : real(r8), intent(inout) :: vmr(:,:,:) ! mixing ratios ( vmr )
526 :
527 : type(physics_buffer_desc), pointer :: pbuf(:)
528 :
529 : ! local vars
530 :
531 : integer :: n, m, mm
532 : integer :: i,k,l
533 : integer :: nstep
534 :
535 72960 : type(ptr2d_t), allocatable :: raer(:) ! aerosol mass, number mixing ratios
536 72960 : type(ptr2d_t), allocatable :: qqcw(:)
537 :
538 145920 : real(r8) :: del_h2so4_aeruptk(ncol,pver)
539 :
540 : real(r8), pointer :: pblh(:) ! pbl height (m)
541 :
542 145920 : real(r8), dimension(ncol) :: wrk
543 : character(len=32) :: name
544 145920 : real(r8) :: dvmrcwdt(ncol,pver,ncnst_tot)
545 145920 : real(r8) :: dvmrdt(ncol,pver,gas_pcnst)
546 145920 : real(r8) :: delta_so4mass(ncol,pver,ncnst_tot)
547 145920 : real(r8) :: wetr_n(pcols,pver,nbins) ! wet radius from CARMA for different bin
548 145920 : real(r8) :: vmrcw(ncol,pver,ncnst_tot) ! cloud-borne aerosol (vmr)
549 145920 : real(r8) :: mmrcw(ncol,pver,ncnst_tot) ! cloud-borne aerosol (mmr)
550 145920 : real(r8) :: raervmr(ncol,pver,ncnst_tot) ! cloud-borne aerosol (vmr)
551 :
552 145920 : real(r8) :: aqso4(ncol,ncnst_tot) ! aqueous phase chemistry
553 145920 : real(r8) :: aqh2so4(ncol,ncnst_tot) ! aqueous phase chemistry
554 145920 : real(r8) :: aqso4_h2o2(ncol) ! SO4 aqueous phase chemistry due to H2O2
555 145920 : real(r8) :: aqso4_o3(ncol) ! SO4 aqueous phase chemistry due to O3
556 145920 : real(r8) :: xphlwc(ncol,pver) ! pH value multiplied by lwc
557 145920 : real(r8) :: nh3_beg(ncol,pver)
558 145920 : real(r8) :: mw_carma(ncnst_tot)
559 : real(r8), pointer :: fldcw(:,:)
560 : real(r8), pointer :: sulfeq(:,:,:)
561 : real(r8) :: wetr(pcols,pver) ! CARMA wet radius in cm
562 : real(r8) :: wetrho(pcols,pver) ! CARMA wet dens
563 72960 : real(r8), allocatable :: rmass(:) ! CARMA rmass
564 :
565 : real(r8) :: old_total_mass
566 : real(r8) :: new_total_mass
567 : real(r8) :: old_total_number
568 :
569 : character(len=32) :: spectype
570 :
571 : character(len=aero_name_len) :: bin_name, shortname
572 : integer :: igroup, ibin, rc, nchr, ierr
573 : character(len=*), parameter :: subname = 'aero_model_gasaerexch'
574 :
575 : !
576 : ! ... initialize nh3
577 : !
578 72960 : if ( nh3_ndx > 0 ) then
579 0 : nh3_beg = vmr(1:ncol,:,nh3_ndx)
580 : end if
581 : !
582 : ! do gas-aerosol exchange (h2so4, msa, nh3 condensation)
583 :
584 72960 : nstep = get_nstep()
585 :
586 : ! calculate tendency due to gas phase chemistry and processes
587 6376704000 : dvmrdt(:ncol,:,:) = (vmr(:ncol,:,:) - vmr0(:ncol,:,:)) / delt
588 5982720 : do m = 1, gas_pcnst
589 91010304 : wrk(:) = 0.0_r8
590 419592960 : do k = 1,pver
591 6376631040 : wrk(:ncol) = wrk(:ncol) + dvmrdt(:ncol,k,m)*adv_mass(m)/mbar(:ncol,k)*pdel(:ncol,k)/gravit
592 : end do
593 5909760 : name = 'GS_'//trim(solsym(m))
594 5982720 : call outfld( name, wrk(:ncol), ncol, lchnk )
595 : enddo
596 :
597 : !
598 : ! Aerosol processes ...
599 : !
600 : allocate( &
601 : rmass(nbins), &
602 : raer(ncnst_tot), &
603 20866560 : qqcw(ncnst_tot), stat=ierr )
604 72960 : if (ierr /= 0) call endrun(subname//': allocate error')
605 :
606 10287360 : mw_carma(:) = 0.0_r8
607 2991360 : do m = 1, nbins ! main loop over aerosol bins
608 : ! dryr is the dry bin radius
609 : ! wetr is the dry bin radius
610 : ! Note: taken here from CARMA pbuf field which may be not any more consistent with changed fields after carma was applied
611 : ! Need to add new code that recalcuates dryr and wetr
612 : ! get bin info
613 2918400 : call rad_cnst_get_info_by_bin(0, m, nspec=nspec(m), bin_name=bin_name)
614 :
615 2918400 : nchr = len_trim(bin_name)-2
616 2918400 : shortname = bin_name(:nchr)
617 :
618 2918400 : call carma_get_group_by_name(shortname, igroup, rc)
619 2918400 : if (rc/=0) then
620 0 : call endrun(subname//': ERROR in carma_get_group_by_name')
621 : end if
622 :
623 2918400 : read(bin_name(nchr+1:),*) ibin
624 :
625 2918400 : call carma_get_wet_radius(state, igroup, ibin, wetr, wetrho, rc) ! m
626 2918400 : if (rc/=0) then
627 0 : call endrun(subname//': ERROR in carma_get_wet_radius')
628 : end if
629 3148953600 : wetr(:ncol,:) = wetr(:ncol,:) * 1.e2_r8 ! cm
630 :
631 2918400 : call carma_get_bin_rmass(igroup, ibin, rmass(m), rc) ! grams
632 2918400 : if (rc/=0) then
633 0 : call endrun(subname//': ERROR in carma_get_bin_rmass')
634 : end if
635 :
636 3475814400 : wetr_n(:,:,m) = wetr(:,:)
637 :
638 : ! Init pointers to mode number and specie mass mixing ratios in
639 : ! intersitial and cloud borne phases.
640 19042560 : do l = 1, nspec(m)
641 10214400 : mm = bin_idx(m, l)
642 10214400 : if (l <= nspec(m)) then
643 10214400 : call rad_cnst_get_bin_props_by_idx(0, m, l,spectype=spectype)
644 10214400 : call rad_cnst_get_bin_mmr_by_idx(0, m, l, 'a', state, pbuf, raer(mm)%fld)
645 10214400 : call rad_cnst_get_bin_mmr_by_idx(0, m, l, 'c', state, pbuf, qqcw(mm)%fld) ! cloud-borne aerosol
646 10214400 : if (trim(spectype) == 'sulfate') then
647 2918400 : mw_carma(mm) = 96._r8
648 : end if
649 10214400 : if (trim(spectype) == 'black-c') then
650 1459200 : mw_carma(mm) = 12._r8
651 : end if
652 10214400 : if (trim(spectype) == 'p-organic') then
653 1459200 : mw_carma(mm) = 12._r8
654 : end if
655 10214400 : if (trim(spectype) == 's-organic') then
656 1459200 : mw_carma(mm) = 250._r8
657 : end if
658 10214400 : if (trim(spectype) == 'dust') then
659 1459200 : mw_carma(mm) = 12._r8
660 : end if
661 10214400 : if (trim(spectype) == 'seasalt') then
662 1459200 : mw_carma(mm) = 57._r8
663 : end if
664 : end if
665 11021337600 : mmrcw(:ncol,:,mm) = qqcw(mm)%fld(:ncol,:)
666 11021337600 : vmrcw(:ncol,:,mm) = qqcw(mm)%fld(:ncol,:)
667 11024256000 : raervmr(:ncol,:,mm) = raer(mm)%fld(:ncol,:)
668 : end do
669 : end do
670 :
671 : ! qqcw2vrm is different from what is done in MAM, here we pass in the fields set by the qqcw and raer pointer
672 : ! for all the CARMA aerosols, species, mmr, and number, vmrcw (kg/kg) -> vmr
673 72960 : call mmr2vmr_carma ( lchnk, vmrcw, mbar, mw_carma, ncol, loffset, rmass )
674 :
675 6376704000 : dvmrdt(:ncol,:,:) = vmr(:ncol,:,:) ! all adveced species no aerosols
676 11021410560 : dvmrcwdt(:ncol,:,:) = vmrcw(:ncol,:,:) ! cloud borne carma aerosol species
677 : ! aqueous chemistry ...
678 :
679 72960 : if( has_sox ) then
680 : call setsox( state, &
681 : ncol, &
682 : lchnk, &
683 : loffset, &
684 : delt, &
685 : pmid, &
686 : pdel, &
687 : tfld, &
688 : mbar, &
689 : cwat, &
690 : cldfr, &
691 : cldnum, &
692 : airdens, &
693 : invariants, &
694 : vmrcw, &
695 : vmr, &
696 : xphlwc, &
697 : aqso4, &
698 : aqh2so4, &
699 : aqso4_h2o2, &
700 : aqso4_o3 &
701 72960 : )
702 :
703 2991360 : do n = 1, nbins
704 13205760 : do l = 1, nspec(n) ! not for total mass or number
705 10214400 : mm = bin_idx(n, l)
706 10214400 : call outfld( trim(fieldname_cw(mm))//'AQSO4', aqso4(:ncol,mm), ncol, lchnk)
707 13132800 : call outfld( trim(fieldname_cw(mm))//'AQH2SO4', aqh2so4(:ncol,mm), ncol, lchnk)
708 : end do
709 : end do
710 :
711 72960 : call outfld( 'AQSO4_H2O2', aqso4_h2o2(:ncol), ncol, lchnk)
712 72960 : call outfld( 'AQSO4_O3', aqso4_o3(:ncol), ncol, lchnk)
713 72960 : call outfld( 'XPH_LWC', xphlwc(:ncol,:), ncol, lchnk )
714 :
715 : endif
716 :
717 : ! Tendency due to aqueous chemistry
718 6376704000 : dvmrdt = (vmr - dvmrdt) / delt
719 11021410560 : dvmrcwdt = (vmrcw - dvmrcwdt) / delt
720 :
721 5982720 : do m = 1, gas_pcnst
722 91010304 : wrk(:) = 0.0_r8
723 419592960 : do k = 1,pver
724 6376631040 : wrk(:ncol) = wrk(:ncol) + dvmrdt(:ncol,k,m) * adv_mass(m)/mbar(:ncol,k)*pdel(:ncol,k)/gravit
725 : end do
726 5909760 : name = 'AQ_'//trim(solsym(m))
727 5982720 : call outfld( name, wrk(:ncol), ncol, lchnk )
728 : enddo
729 :
730 : ! do gas-aerosol exchange (h2so4, msa, nh3 condensation)
731 :
732 72960 : if (h2so4_ndx > 0) then
733 78723840 : del_h2so4_aeruptk(1:ncol,:) = vmr(1:ncol,:,h2so4_ndx)
734 : else
735 0 : del_h2so4_aeruptk(:,:) = 0.0_r8
736 : endif
737 :
738 : ! need to transform raer to raervmr from CARMA, routine requires vmr, note number wil not be changed here
739 72960 : call mmr2vmr_carma ( lchnk, raervmr, mbar, mw_carma, ncol, loffset, rmass)
740 :
741 : call carma_aero_gasaerexch_sub( state, &
742 : pbuf, lchnk, ncol, nstep, &
743 : loffset, delt, mbar , &
744 : tfld, pmid, pdel, &
745 : qh2o, troplev, &
746 : vmr, raervmr, &
747 72960 : wetr_n )
748 :
749 : ! note vmr2qqcw does not change qqcw pointer (different than in MAM)
750 72960 : call vmr2mmr_carma ( lchnk, vmrcw, mbar, mw_carma, ncol, loffset, rmass )
751 :
752 : !vmrcw in kg/kg
753 : ! change pointer value for total mmr and number. In order to do this correctly
754 : ! only mass has to be added to each bin (not number). This will require redistributing
755 : ! mass to different bins. Here, we change both mass and number until we have a better
756 : ! solution.
757 11021410560 : delta_so4mass(:,:,:) = 0.0_r8
758 2991360 : do m = 1, nbins
759 13205760 : do l = 1, nspec(m) ! for sulfate only
760 10214400 : mm = bin_idx(m, l)
761 : ! sulfate mass that needs to be added to the total mass
762 10214400 : call rad_cnst_get_bin_props_by_idx(0, m, l,spectype=spectype)
763 13132800 : if (trim(spectype) == 'sulfate') then
764 : ! only do loop if vmrcw has changed
765 207206400 : do k=1,pver
766 3148953600 : do i=1,ncol
767 3146035200 : if (vmrcw(i,k,mm) .gt. mmrcw(i,k,mm) .and. mmrcw(i,k,mm) /= 0.0_r8) then
768 44420676 : delta_so4mass(i,k,mm) = ( vmrcw(i,k,mm) - mmrcw(i,k,mm) )
769 : else
770 2897326524 : delta_so4mass(i,k,mm) = 0.0_r8
771 : end if
772 : end do
773 : end do
774 : end if
775 : end do
776 : end do
777 :
778 2991360 : do m = 1, nbins
779 13205760 : do l = 1, nspec(m) ! for sulfate only
780 10214400 : mm = bin_idx(m, l)
781 11021337600 : qqcw(mm)%fld(:ncol,:) = vmrcw(:ncol,:,mm)
782 11024256000 : call outfld( trim(fieldname_cw(mm)), qqcw(mm)%fld(:ncol,:), ncol, lchnk)
783 : end do
784 : end do
785 :
786 :
787 72960 : end subroutine aero_model_gasaerexch
788 :
789 : !=============================================================================
790 : !=============================================================================
791 72960 : subroutine aero_model_emissions( state, cam_in )
792 :
793 : ! Arguments:
794 :
795 : type(physics_state), intent(in) :: state ! Physics state variables
796 : type(cam_in_t), intent(inout) :: cam_in ! import state
797 :
798 72960 : end subroutine aero_model_emissions
799 :
800 :
801 : !===============================================================================
802 : !===============================================================================
803 : ! private methods
804 :
805 :
806 : !=============================================================================
807 : !=============================================================================
808 145920 : subroutine surf_area_dens( state, pbuf, ncol, mmr, beglev, endlev, sad, reff, sfc, dm_aer )
809 : use mo_constants, only: pi
810 : use carma_intr, only: carma_effecitive_radius
811 :
812 : ! dummy args
813 : type(physics_state), intent(in) :: state ! Physics state variables
814 : type(physics_buffer_desc), pointer :: pbuf(:)
815 : integer, intent(in) :: ncol
816 : real(r8), intent(in) :: mmr(:,:,:)
817 : integer, intent(in) :: beglev(:)
818 : integer, intent(in) :: endlev(:)
819 : real(r8), intent(out) :: sad(:,:) ! bulk surface area density in cm2/cm3 from beglev to endlev, zero elsewhere
820 : real(r8), intent(out) :: reff(:,:) ! bulk effective radius in cm from beglev to endlev, zero elsewhere
821 : real(r8), optional, intent(out) :: sfc(:,:,:) ! surface area density per bin
822 : real(r8), optional, intent(out) :: dm_aer(:,:,:) ! diameter per bin
823 :
824 : ! local vars
825 : real(r8) :: reffaer(pcols,pver) ! bulk effective radius in cm
826 :
827 291840 : real(r8) :: sad_bin(pcols,pver,nbins)
828 : integer :: icol, ilev, ibin, ispec !!, reff_pbf_ndx
829 : real(r8) :: chm_mass, tot_mass
830 : character(len=32) :: spectype
831 : real(r8) :: wetr(pcols,pver) ! CARMA bin wet radius in cm
832 : real(r8) :: wetrho(pcols,pver) ! CARMA bin wet density
833 : real(r8) :: sad_carma(pcols,pver) ! CARMA bin wet surface area density in cm2/cm3
834 145920 : real(r8), pointer :: aer_bin_mmr(:,:)
835 :
836 : character(len=aero_name_len) :: bin_name, shortname
837 : integer :: igroup, indxbin, rc, nchr
838 :
839 173790720 : sad = 0._r8
840 173790720 : reff = 0._r8
841 :
842 : !
843 : ! Compute surface aero for each bin.
844 : ! Total over all bins as the surface area for chemical reactions.
845 : !
846 :
847 145920 : reffaer = carma_effecitive_radius(state)
848 :
849 173790720 : sad = 0._r8
850 6951774720 : sad_bin = 0._r8
851 173790720 : reff = 0._r8
852 :
853 5982720 : do ibin=1,nbins ! loop over aerosol bins
854 5836800 : call rad_cnst_get_info_by_bin(0, ibin, bin_name=bin_name)
855 :
856 5836800 : nchr = len_trim(bin_name)-2
857 5836800 : shortname = bin_name(:nchr)
858 :
859 5836800 : call carma_get_group_by_name(shortname, igroup, rc)
860 :
861 5836800 : read(bin_name(nchr+1:),*) indxbin
862 :
863 5836800 : call carma_get_wet_radius(state, igroup, indxbin, wetr, wetrho, rc) ! m
864 6297907200 : wetr(:ncol,:) = wetr(:ncol,:) * 1.e2_r8 ! cm
865 5836800 : call carma_get_sad(state, igroup, indxbin, sad_carma, rc)
866 :
867 5836800 : if (present(dm_aer)) then
868 3148953600 : dm_aer(:ncol,:,ibin) = 2._r8 * wetr(:ncol,:) ! convert wet radius (cm) to wet diameter (cm)
869 : endif
870 6309726720 : sad_bin(:ncol,:,ibin) = sad_carma(:ncol,:) ! cm^2/cm^3
871 : end do
872 :
873 2247168 : do icol = 1,ncol
874 75790848 : do ilev = beglev(icol),endlev(icol)
875 3015290880 : do ibin=1,nbins ! loop over aerosol bins
876 : !
877 : ! compute a mass weighting of the number
878 : !
879 2941747200 : tot_mass = 0._r8
880 2941747200 : chm_mass = 0._r8
881 13237862400 : do ispec=1,nspec(ibin)
882 :
883 10296115200 : call rad_cnst_get_bin_mmr_by_idx(0, ibin, ispec, 'a', state, pbuf, aer_bin_mmr)
884 :
885 10296115200 : tot_mass = tot_mass + aer_bin_mmr(icol,ilev)
886 :
887 10296115200 : call rad_cnst_get_bin_props_by_idx(0, ibin, ispec, spectype=spectype)
888 :
889 : if ( trim(spectype) == 'sulfate' .or. &
890 : trim(spectype) == 's-organic' .or. &
891 : trim(spectype) == 'p-organic' .or. &
892 10296115200 : trim(spectype) == 'black-c' .or. &
893 2941747200 : trim(spectype) == 'ammonium') then
894 7354368000 : chm_mass = chm_mass + aer_bin_mmr(icol,ilev)
895 : end if
896 :
897 : end do
898 3015290880 : if ( tot_mass > 0._r8 ) then
899 : ! surface area density
900 2861412147 : sad_bin(icol,ilev,ibin) = chm_mass / tot_mass * sad_bin(icol,ilev,ibin) ! cm^2/cm^3
901 : else
902 80335053 : sad_bin(icol,ilev,ibin) = 0._r8
903 : end if
904 : end do
905 3015290880 : sad(icol,ilev) = sum(sad_bin(icol,ilev,:))
906 75644928 : reff(icol,ilev) = reffaer(icol,ilev)
907 :
908 : end do
909 : end do
910 :
911 145920 : if (present(sfc)) then
912 3475887360 : sfc(:,:,:) = sad_bin(:,:,:)
913 : endif
914 :
915 291840 : end subroutine surf_area_dens
916 :
917 : !=============================================================================
918 145920 : subroutine mmr2vmr_carma(lchnk, vmr, mbar, mw_carma, ncol, im, rmass)
919 : !-----------------------------------------------------------------
920 : ! ... Xfrom from mass to volume mixing ratio
921 : !-----------------------------------------------------------------
922 :
923 : implicit none
924 :
925 : !-----------------------------------------------------------------
926 : ! ... Dummy args
927 : !-----------------------------------------------------------------
928 : integer, intent(in) :: lchnk, ncol, im
929 : real(r8), intent(in) :: mbar(ncol,pver)
930 : real(r8), intent(in) :: rmass(nbins)
931 : real(r8), intent(in) :: mw_carma(ncnst_tot)
932 : real(r8), intent(inout) :: vmr(ncol,pver,ncnst_tot)
933 :
934 : !-----------------------------------------------------------------
935 : ! ... Local variables
936 : !-----------------------------------------------------------------
937 : integer :: k, m, mm, l
938 :
939 5982720 : do m = 1, nbins
940 26411520 : do l = 1, nspec(m) ! for each species, not total mmr or number, information of mw are missing
941 20428800 : mm = bin_idx(m, l)
942 1456281600 : do k=1,pver
943 22042675200 : vmr(:ncol,k,mm) = mbar(:ncol,k) * vmr(:ncol,k,mm) / mw_carma(mm)
944 : end do
945 : end do
946 : end do
947 :
948 145920 : end subroutine mmr2vmr_carma
949 : !=============================================================================
950 :
951 : !=============================================================================
952 72960 : subroutine vmr2mmr_carma ( lchnk, vmr, mbar, mw_carma, ncol, im, rmass )
953 : !-----------------------------------------------------------------
954 : ! ... Xfrom from volume to mass mixing ratio
955 : !-----------------------------------------------------------------
956 :
957 : implicit none
958 :
959 : !-----------------------------------------------------------------
960 : ! ... Dummy args
961 : !-----------------------------------------------------------------
962 : integer, intent(in) :: lchnk, ncol, im
963 : real(r8), intent(in) :: mbar(ncol,pver)
964 : real(r8), intent(in) :: rmass(nbins)
965 : real(r8), intent(inout) :: vmr(ncol,pver,ncnst_tot)
966 : real(r8), intent(in) :: mw_carma(ncnst_tot)
967 :
968 : !-----------------------------------------------------------------
969 : ! ... Local variables
970 : !-----------------------------------------------------------------
971 : integer :: k, m, mm, l
972 : !-----------------------------------------------------------------
973 : ! ... The non-group species
974 : !-----------------------------------------------------------------
975 2991360 : do m = 1, nbins
976 13205760 : do l = 1, nspec(m) ! for each species, not total mmr or number, information of mw are missing
977 10214400 : mm = bin_idx(m, l)
978 728140800 : do k=1,pver
979 11021337600 : vmr(:ncol,k,mm) = mw_carma(mm) * vmr(:ncol,k,mm) / mbar(:ncol,k)
980 : end do
981 : end do
982 : end do
983 :
984 72960 : end subroutine vmr2mmr_carma
985 :
986 0 : end module aero_model
|
