Line data Source code
1 : module aerosol_optics_cam
2 : use shr_kind_mod, only: r8 => shr_kind_r8
3 : use shr_kind_mod, only: cl => shr_kind_cl
4 : use cam_logfile, only: iulog
5 : use radconstants, only: nswbands, nlwbands, idx_sw_diag, idx_uv_diag, idx_nir_diag
6 : use radconstants, only: get_lw_spectral_boundaries
7 : use phys_prop, only: ot_length
8 : use physics_types,only: physics_state
9 : use physics_buffer,only: physics_buffer_desc
10 : use ppgrid, only: pcols, pver
11 : use physconst, only: rga, rair
12 : use cam_abortutils, only: endrun
13 : use spmd_utils, only: masterproc
14 : use rad_constituents, only: n_diag, rad_cnst_get_call_list
15 : use cam_history, only: addfld, add_default, outfld, horiz_only, fieldname_len
16 : use cam_history_support, only: fillvalue
17 :
18 : use tropopause, only : tropopause_findChemTrop
19 :
20 : use aerosol_properties_mod, only: aerosol_properties
21 : use modal_aerosol_properties_mod, only: modal_aerosol_properties
22 :
23 : use aerosol_state_mod, only: aerosol_state
24 : use modal_aerosol_state_mod,only: modal_aerosol_state
25 :
26 : use aerosol_optics_mod, only: aerosol_optics
27 : use refractive_aerosol_optics_mod, only: refractive_aerosol_optics
28 :
29 : implicit none
30 :
31 : private
32 :
33 : public :: aerosol_optics_cam_readnl
34 : public :: aerosol_optics_cam_init
35 : public :: aerosol_optics_cam_final
36 : public :: aerosol_optics_cam_sw
37 : public :: aerosol_optics_cam_lw
38 :
39 : type aero_props_t
40 : class(aerosol_properties), pointer :: obj => null()
41 : end type aero_props_t
42 : type aero_state_t
43 : class(aerosol_state), pointer :: obj => null()
44 : end type aero_state_t
45 :
46 : type(aero_props_t), allocatable :: aero_props(:) ! array of aerosol properties objects to allow for
47 : ! multiple aerosol representations in the same sim
48 : ! such as MAM and CARMA
49 :
50 : ! refractive index for water read in read_water_refindex
51 : complex(r8) :: crefwsw(nswbands) = -huge(1._r8) ! complex refractive index for water visible
52 : complex(r8) :: crefwlw(nlwbands) = -huge(1._r8) ! complex refractive index for water infrared
53 : character(len=cl) :: water_refindex_file = 'NONE' ! full pathname for water refractive index dataset
54 :
55 : logical :: modal_active = .false.
56 : integer :: num_aero_models = 0
57 : integer :: lw10um_indx = -1 ! wavelength index corresponding to 10 microns
58 : real(r8), parameter :: lw10um = 10._r8 ! microns
59 :
60 : character(len=4) :: diag(0:n_diag) = (/' ','_d1 ','_d2 ','_d3 ','_d4 ','_d5 ', '_d6 ','_d7 ','_d8 ','_d9 ','_d10'/)
61 :
62 : type out_name
63 : character(len=fieldname_len), allocatable :: name(:) ! nbins
64 : end type out_name
65 :
66 : type(out_name), allocatable :: burden_fields(:) ! num_aero_models
67 : type(out_name), allocatable :: aodbin_fields(:)
68 : type(out_name), allocatable :: aoddust_fields(:)
69 : type(out_name), allocatable :: burdendn_fields(:) ! num_aero_models
70 : type(out_name), allocatable :: aodbindn_fields(:)
71 : type(out_name), allocatable :: aoddustdn_fields(:)
72 :
73 : contains
74 :
75 : !===============================================================================
76 1536 : subroutine aerosol_optics_cam_readnl(nlfile)
77 : use namelist_utils, only : find_group_name
78 : use spmd_utils, only : mpicom, masterprocid, mpi_character, mpi_success
79 :
80 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
81 :
82 : integer :: unitn, ierr
83 : character(len=cl) :: errmsg
84 : character(len=*), parameter :: subname = 'aerosol_optics_cam_readnl'
85 :
86 : ! ===================
87 : ! Namelist definition
88 : ! ===================
89 : namelist /aerosol_optics_nl/ water_refindex_file
90 :
91 : ! =============
92 : ! Read namelist
93 : ! =============
94 1536 : if (masterproc) then
95 2 : open( newunit=unitn, file=trim(nlfile), status='old' )
96 2 : call find_group_name(unitn, 'aerosol_optics_nl', status=ierr)
97 2 : if (ierr == 0) then
98 2 : read(unitn, aerosol_optics_nl, iostat=ierr)
99 2 : if (ierr /= 0) then
100 0 : write(errmsg,'(2a,i10)') subname,':: ERROR reading namelist, error code: ',ierr
101 0 : call endrun(errmsg)
102 : end if
103 : end if
104 2 : close(unitn)
105 : end if
106 :
107 : ! ============================
108 : ! Broadcast namelist variables
109 : ! ============================
110 1536 : call mpi_bcast(water_refindex_file, len(water_refindex_file), mpi_character, masterprocid, mpicom, ierr)
111 1536 : if (ierr/=mpi_success) then
112 0 : call endrun(subname // ':: ERROR mpi_bcast '//trim(water_refindex_file))
113 : end if
114 :
115 1536 : if (masterproc) then
116 2 : write(iulog,*) subname,': water_refindex_file = ',trim(water_refindex_file)
117 : end if
118 :
119 1536 : end subroutine aerosol_optics_cam_readnl
120 :
121 : !===============================================================================
122 1536 : subroutine aerosol_optics_cam_init
123 : use rad_constituents, only: rad_cnst_get_info
124 : use phys_control, only: phys_getopts
125 : use ioFileMod, only: getfil
126 :
127 : character(len=*), parameter :: prefix = 'aerosol_optics_cam_init: '
128 : integer :: nmodes=0, iaermod, istat, ilist, i
129 :
130 : logical :: call_list(0:n_diag)
131 : real(r8) :: lwavlen_lo(nlwbands), lwavlen_hi(nlwbands)
132 : integer :: m, n
133 :
134 : character(len=fieldname_len) :: fldname
135 : character(len=128) :: lngname
136 : logical :: history_aero_optics ! output aerosol optics diagnostics
137 : logical :: history_amwg ! output the variables used by the AMWG diag package
138 : logical :: history_dust ! output dust diagnostics
139 :
140 : character(len=cl) :: locfile
141 :
142 : call phys_getopts(history_amwg_out = history_amwg, &
143 : history_aero_optics_out = history_aero_optics, &
144 1536 : history_dust_out = history_dust )
145 :
146 1536 : num_aero_models = 0
147 :
148 1536 : call rad_cnst_get_info(0, nmodes=nmodes)
149 1536 : modal_active = nmodes>0
150 :
151 1536 : if (modal_active) then
152 1536 : num_aero_models = num_aero_models+1 ! count aerosol models
153 : end if
154 :
155 1536 : if (num_aero_models>0) then
156 6144 : allocate(aero_props(num_aero_models), stat=istat)
157 1536 : if (istat/=0) then
158 0 : call endrun(prefix//'array allocation error: aero_props')
159 : end if
160 : end if
161 :
162 1536 : iaermod = 0
163 :
164 1536 : if (modal_active) then
165 1536 : iaermod = iaermod+1
166 1536 : aero_props(iaermod)%obj => modal_aerosol_properties()
167 : end if
168 :
169 1536 : if (water_refindex_file=='NONE') then
170 0 : call endrun(prefix//'water_refindex_file must be specified')
171 : else
172 1536 : call getfil(water_refindex_file, locfile)
173 1536 : call read_water_refindex(locfile)
174 : end if
175 :
176 1536 : call get_lw_spectral_boundaries(lwavlen_lo, lwavlen_hi, units='um')
177 26112 : do i = 1,nlwbands
178 26112 : if ((lwavlen_lo(i)<=lw10um) .and. (lwavlen_hi(i)>=lw10um)) then
179 1536 : lw10um_indx = i ! index corresponding to 10 microns
180 : end if
181 : end do
182 1536 : call rad_cnst_get_call_list(call_list)
183 :
184 18432 : do ilist = 0, n_diag
185 18432 : if (call_list(ilist)) then
186 : call addfld ('EXTINCT'//diag(ilist), (/ 'lev' /), 'A','/m',&
187 3072 : 'Aerosol extinction 550 nm, day only', flag_xyfill=.true.)
188 : call addfld ('EXTINCTUV'//diag(ilist), (/ 'lev' /), 'A','/m',&
189 3072 : 'Aerosol extinction 350 nm, day only', flag_xyfill=.true.)
190 : call addfld ('EXTINCTNIR'//diag(ilist), (/ 'lev' /), 'A','/m',&
191 3072 : 'Aerosol extinction 1020 nm, day only', flag_xyfill=.true.)
192 : call addfld ('ABSORB'//diag(ilist), (/ 'lev' /), 'A','/m',&
193 3072 : 'Aerosol absorption, day only', flag_xyfill=.true.)
194 : call addfld ('AODVIS'//diag(ilist), horiz_only, 'A',' ', &
195 1536 : 'Aerosol optical depth 550 nm', flag_xyfill=.true.)
196 : call addfld ('AODVISst'//diag(ilist), horiz_only, 'A',' ', &
197 1536 : 'Stratospheric aerosol optical depth 550 nm, day only', flag_xyfill=.true.)
198 : call addfld ('AODNIRst'//diag(ilist), horiz_only, 'A',' ', &
199 1536 : 'Stratospheric aerosol optical depth 1020 nm, day only',flag_xyfill=.true.)
200 : call addfld ('AODUVst'//diag(ilist), horiz_only, 'A',' ', &
201 1536 : 'Stratospheric aerosol optical depth 350 nm, day only', flag_xyfill=.true.)
202 : call addfld ('AODUV'//diag(ilist), horiz_only, 'A',' ', &
203 1536 : 'Aerosol optical depth 350 nm, day only', flag_xyfill=.true.)
204 : call addfld ('AODNIR'//diag(ilist), horiz_only, 'A',' ', &
205 1536 : 'Aerosol optical depth 1020 nm, day only',flag_xyfill=.true.)
206 : call addfld ('AODABS'//diag(ilist), horiz_only, 'A',' ', &
207 1536 : 'Aerosol absorption optical depth 550 nm, day only', flag_xyfill=.true.)
208 : call addfld ('AODxASYM'//diag(ilist), horiz_only, 'A',' ', &
209 1536 : 'Aerosol optical depth 550 * asymmetry factor, day only', flag_xyfill=.true.)
210 : call addfld ('EXTxASYM'//diag(ilist), (/ 'lev' /), 'A',' ', &
211 3072 : 'extinction 550 nm * asymmetry factor, day only', flag_xyfill=.true.)
212 : call addfld ('AODTOT'//diag(ilist), horiz_only, 'A','1',&
213 1536 : 'Aerosol optical depth summed over all sw wavelengths', flag_xyfill=.true.)
214 :
215 : call addfld ('EXTINCTdn'//diag(ilist), (/ 'lev' /), 'A','/m',&
216 3072 : 'Aerosol extinction 550 nm, day only', flag_xyfill=.true.)
217 : call addfld ('EXTINCTUVdn'//diag(ilist), (/ 'lev' /), 'A','/m',&
218 3072 : 'Aerosol extinction 350 nm, day only', flag_xyfill=.true.)
219 : call addfld ('EXTINCTNIRdn'//diag(ilist), (/ 'lev' /), 'A','/m',&
220 3072 : 'Aerosol extinction 1020 nm, day only', flag_xyfill=.true.)
221 : call addfld ('ABSORBdn'//diag(ilist), (/ 'lev' /), 'A','/m',&
222 3072 : 'Aerosol absorption, day only', flag_xyfill=.true.)
223 : call addfld ('AODVISdn'//diag(ilist), horiz_only, 'A',' ', &
224 1536 : 'Aerosol optical depth 550 nm', flag_xyfill=.true.)
225 : call addfld ('AODVISstdn'//diag(ilist), horiz_only, 'A',' ', &
226 1536 : 'Stratospheric aerosol optical depth 550 nm, day only', flag_xyfill=.true.)
227 : call addfld ('AODNIRstdn'//diag(ilist), horiz_only, 'A',' ', &
228 1536 : 'Stratospheric aerosol optical depth 1020 nm, day only', flag_xyfill=.true.)
229 : call addfld ('AODUVstdn'//diag(ilist), horiz_only, 'A',' ', &
230 1536 : 'Stratospheric aerosol optical depth 350 nm, day only', flag_xyfill=.true.)
231 : call addfld ('AODUVdn'//diag(ilist), horiz_only, 'A',' ', &
232 1536 : 'Aerosol optical depth 350 nm, day only', flag_xyfill=.true.)
233 : call addfld ('AODNIRdn'//diag(ilist), horiz_only, 'A',' ', &
234 1536 : 'Aerosol optical depth 1020 nm, day only', flag_xyfill=.true.)
235 : call addfld ('AODABSdn'//diag(ilist), horiz_only, 'A',' ', &
236 1536 : 'Aerosol absorption optical depth 550 nm, day only', flag_xyfill=.true.)
237 : call addfld ('AODxASYMdn'//diag(ilist), horiz_only, 'A',' ', &
238 1536 : 'Aerosol optical depth 550 * asymmetry factor, day only', flag_xyfill=.true.)
239 : call addfld ('EXTxASYMdn'//diag(ilist), (/ 'lev' /), 'A',' ', &
240 3072 : 'extinction 550 nm * asymmetry factor, day only', flag_xyfill=.true.)
241 : call addfld ('AODTOTdn'//diag(ilist), horiz_only, 'A','1',&
242 1536 : 'Aerosol optical depth summed over all sw wavelengths, day only')
243 :
244 1536 : if (lw10um_indx>0) then
245 : call addfld('AODABSLW'//diag(ilist), (/ 'lev' /), 'A','/m',&
246 3072 : 'Aerosol long-wave absorption optical depth at 10 microns')
247 : end if
248 : call addfld ('TOTABSLW'//diag(ilist), (/ 'lev' /), 'A',' ', &
249 3072 : 'LW Aero total abs')
250 :
251 1536 : if (ilist>0 .and. history_aero_optics) then
252 0 : call add_default ('EXTINCT'//diag(ilist), 1, ' ')
253 0 : call add_default ('ABSORB'//diag(ilist), 1, ' ')
254 0 : call add_default ('AODVIS'//diag(ilist), 1, ' ')
255 0 : call add_default ('AODVISst'//diag(ilist), 1, ' ')
256 0 : call add_default ('AODABS'//diag(ilist), 1, ' ')
257 : end if
258 :
259 : end if
260 : end do
261 :
262 1536 : if (num_aero_models>0) then
263 :
264 6144 : allocate(burden_fields(num_aero_models), stat=istat)
265 1536 : if (istat/=0) then
266 0 : call endrun(prefix//'array allocation error: burden_fields')
267 : end if
268 6144 : allocate(aodbin_fields(num_aero_models), stat=istat)
269 1536 : if (istat/=0) then
270 0 : call endrun(prefix//'array allocation error: aodbin_fields')
271 : end if
272 6144 : allocate(aoddust_fields(num_aero_models), stat=istat)
273 1536 : if (istat/=0) then
274 0 : call endrun(prefix//'array allocation error: aoddust_fields')
275 : end if
276 :
277 6144 : allocate(burdendn_fields(num_aero_models), stat=istat)
278 1536 : if (istat/=0) then
279 0 : call endrun(prefix//'array allocation error: burdendn_fields')
280 : end if
281 6144 : allocate(aodbindn_fields(num_aero_models), stat=istat)
282 1536 : if (istat/=0) then
283 0 : call endrun(prefix//'array allocation error: aodbindn_fields')
284 : end if
285 6144 : allocate(aoddustdn_fields(num_aero_models), stat=istat)
286 1536 : if (istat/=0) then
287 0 : call endrun(prefix//'array allocation error: aoddustdn_fields')
288 : end if
289 :
290 3072 : do n = 1,num_aero_models
291 :
292 4608 : allocate(burden_fields(n)%name(aero_props(n)%obj%nbins()), stat=istat)
293 1536 : if (istat/=0) then
294 0 : call endrun(prefix//'array allocation error: burden_fields(n)%name')
295 : end if
296 4608 : allocate(aodbin_fields(n)%name(aero_props(n)%obj%nbins()), stat=istat)
297 1536 : if (istat/=0) then
298 0 : call endrun(prefix//'array allocation error: aodbin_fields(n)%name')
299 : end if
300 4608 : allocate(aoddust_fields(n)%name(aero_props(n)%obj%nbins()), stat=istat)
301 1536 : if (istat/=0) then
302 0 : call endrun(prefix//'array allocation error: aoddust_fields(n)%name')
303 : end if
304 :
305 4608 : allocate(burdendn_fields(n)%name(aero_props(n)%obj%nbins()), stat=istat)
306 1536 : if (istat/=0) then
307 0 : call endrun(prefix//'array allocation error: burdendn_fields(n)%name')
308 : end if
309 4608 : allocate(aodbindn_fields(n)%name(aero_props(n)%obj%nbins()), stat=istat)
310 1536 : if (istat/=0) then
311 0 : call endrun(prefix//'array allocation error: aodbindn_fields(n)%name')
312 : end if
313 4608 : allocate(aoddustdn_fields(n)%name(aero_props(n)%obj%nbins()), stat=istat)
314 1536 : if (istat/=0) then
315 0 : call endrun(prefix//'array allocation error: aoddustdn_fields(n)%name')
316 : end if
317 :
318 9216 : do m = 1, aero_props(n)%obj%nbins()
319 :
320 6144 : write(fldname,'(a,i2.2)') 'BURDEN', m
321 6144 : burden_fields(n)%name(m) = fldname
322 6144 : write(lngname,'(a,i2.2)') 'Aerosol burden bin ', m
323 6144 : call addfld (fldname, horiz_only, 'A', 'kg/m2', lngname, flag_xyfill=.true.)
324 6144 : if (history_aero_optics) then
325 0 : call add_default (fldname, 1, ' ')
326 : end if
327 :
328 6144 : fldname = 'AOD_'//trim(aero_props(n)%obj%bin_name(0,m))
329 6144 : aodbin_fields(n)%name(m) = fldname
330 6144 : lngname = 'Aerosol optical depth, day only, 550 nm, '//trim(aero_props(n)%obj%bin_name(0,m))
331 6144 : call addfld (aodbin_fields(n)%name(m), horiz_only, 'A', ' ', lngname, flag_xyfill=.true.)
332 6144 : if (history_aero_optics) then
333 0 : call add_default (fldname, 1, ' ')
334 : end if
335 :
336 6144 : write(fldname,'(a,i2.2)') 'AODDUST', m
337 6144 : aoddust_fields(n)%name(m) = fldname
338 6144 : write(lngname,'(a,i2,a)') 'Aerosol optical depth, day only, 550 nm mode ',m,' from dust'
339 6144 : call addfld (aoddust_fields(n)%name(m), horiz_only, 'A', ' ', lngname, flag_xyfill=.true.)
340 6144 : if (history_aero_optics) then
341 0 : call add_default (fldname, 1, ' ')
342 : end if
343 :
344 6144 : write(fldname,'(a,i2.2)') 'BURDENdn', m
345 6144 : burdendn_fields(n)%name(m) = fldname
346 6144 : write(lngname,'(a,i2)') 'Aerosol burden, day night, bin ', m
347 6144 : call addfld (burdendn_fields(n)%name(m), horiz_only, 'A', 'kg/m2', lngname, flag_xyfill=.true.)
348 6144 : if (history_aero_optics) then
349 0 : call add_default (fldname, 1, ' ')
350 : end if
351 :
352 6144 : fldname = 'AODdn_'//trim(aero_props(n)%obj%bin_name(0,m))
353 6144 : aodbindn_fields(n)%name(m) = fldname
354 6144 : lngname = 'Aerosol optical depth 550 nm, day night, '//trim(aero_props(n)%obj%bin_name(0,m))
355 6144 : call addfld (aodbindn_fields(n)%name(m), horiz_only, 'A', ' ', lngname, flag_xyfill=.true.)
356 6144 : if (history_aero_optics) then
357 0 : call add_default (fldname, 1, ' ')
358 : end if
359 :
360 6144 : write(fldname,'(a,i2.2)') 'AODdnDUST', m
361 6144 : aoddustdn_fields(n)%name(m) = fldname
362 6144 : write(lngname,'(a,i2,a)') 'Aerosol optical depth 550 nm, day night, bin ',m,' from dust'
363 6144 : call addfld (aoddustdn_fields(n)%name(m), horiz_only, 'A', ' ', lngname, flag_xyfill=.true.)
364 7680 : if (history_aero_optics) then
365 0 : call add_default (fldname, 1, ' ')
366 : end if
367 :
368 : end do
369 :
370 : end do
371 :
372 : end if
373 :
374 : call addfld ('AODDUST', horiz_only, 'A',' ', 'Aerosol optical depth 550 nm from dust, day only', &
375 1536 : flag_xyfill=.true.)
376 : call addfld ('AODSO4', horiz_only, 'A',' ', 'Aerosol optical depth 550 nm from SO4, day only', &
377 1536 : flag_xyfill=.true.)
378 : call addfld ('AODPOM', horiz_only, 'A',' ', 'Aerosol optical depth 550 nm from POM, day only', &
379 1536 : flag_xyfill=.true.)
380 : call addfld ('AODSOA', horiz_only, 'A',' ', 'Aerosol optical depth 550 nm from SOA, day only', &
381 1536 : flag_xyfill=.true.)
382 : call addfld ('AODBC', horiz_only, 'A',' ', 'Aerosol optical depth 550 nm from BC, day only', &
383 1536 : flag_xyfill=.true.)
384 : call addfld ('AODSS', horiz_only, 'A',' ', 'Aerosol optical depth 550 nm from seasalt, day only', &
385 1536 : flag_xyfill=.true.)
386 : call addfld ('AODABSBC', horiz_only, 'A',' ', 'Aerosol absorption optical depth 550 nm from BC, day only',&
387 1536 : flag_xyfill=.true.)
388 : call addfld ('BURDENDUST', horiz_only, 'A','kg/m2', 'Dust aerosol burden, day only' , &
389 1536 : flag_xyfill=.true.)
390 : call addfld ('BURDENSO4', horiz_only, 'A','kg/m2', 'Sulfate aerosol burden, day only' , &
391 1536 : flag_xyfill=.true.)
392 : call addfld ('BURDENPOM', horiz_only, 'A','kg/m2', 'POM aerosol burden, day only' , &
393 1536 : flag_xyfill=.true.)
394 : call addfld ('BURDENSOA', horiz_only, 'A','kg/m2', 'SOA aerosol burden, day only' , &
395 1536 : flag_xyfill=.true.)
396 : call addfld ('BURDENBC', horiz_only, 'A','kg/m2', 'Black carbon aerosol burden, day only', &
397 1536 : flag_xyfill=.true.)
398 : call addfld ('BURDENSEASALT', horiz_only, 'A','kg/m2', 'Seasalt aerosol burden, day only' , &
399 1536 : flag_xyfill=.true.)
400 : call addfld ('SSAVIS', horiz_only, 'A',' ', 'Aerosol single-scatter albedo, day only', &
401 1536 : flag_xyfill=.true.)
402 :
403 : call addfld ('AODDUSTdn', horiz_only, 'A',' ', 'Aerosol optical depth 550 nm from dust, day night', &
404 1536 : flag_xyfill=.true.)
405 : call addfld ('AODSO4dn', horiz_only, 'A',' ', 'Aerosol optical depth 550 nm from SO4, day night', &
406 1536 : flag_xyfill=.true.)
407 : call addfld ('AODPOMdn', horiz_only, 'A',' ', 'Aerosol optical depth 550 nm from POM, day night', &
408 1536 : flag_xyfill=.true.)
409 : call addfld ('AODSOAdn', horiz_only, 'A',' ', 'Aerosol optical depth 550 nm from SOA, day night', &
410 1536 : flag_xyfill=.true.)
411 : call addfld ('AODBCdn', horiz_only, 'A',' ', 'Aerosol optical depth 550 nm from BC, day night', &
412 1536 : flag_xyfill=.true.)
413 : call addfld ('AODSSdn', horiz_only, 'A',' ', 'Aerosol optical depth 550 nm from seasalt, day night', &
414 1536 : flag_xyfill=.true.)
415 : call addfld ('AODABSBCdn', horiz_only, 'A',' ', 'Aerosol absorption optical depth 550 nm from BC, day night',&
416 1536 : flag_xyfill=.true.)
417 : call addfld ('BURDENDUSTdn', horiz_only, 'A','kg/m2', 'Dust aerosol burden, day night' , &
418 1536 : flag_xyfill=.true.)
419 : call addfld ('BURDENSO4dn', horiz_only, 'A','kg/m2', 'Sulfate aerosol burden, day night' , &
420 1536 : flag_xyfill=.true.)
421 : call addfld ('BURDENPOMdn', horiz_only, 'A','kg/m2', 'POM aerosol burden, day night' , &
422 1536 : flag_xyfill=.true.)
423 : call addfld ('BURDENSOAdn', horiz_only, 'A','kg/m2', 'SOA aerosol burden, day night' , &
424 1536 : flag_xyfill=.true.)
425 : call addfld ('BURDENBCdn', horiz_only, 'A','kg/m2', 'Black carbon aerosol burden, day night', &
426 1536 : flag_xyfill=.true.)
427 : call addfld ('BURDENSEASALTdn', horiz_only, 'A','kg/m2', 'Seasalt aerosol burden, day night' , &
428 1536 : flag_xyfill=.true.)
429 : call addfld ('SSAVISdn', horiz_only, 'A',' ', 'Aerosol single-scatter albedo, day night', &
430 1536 : flag_xyfill=.true.)
431 :
432 1536 : if (history_amwg) then
433 1536 : call add_default ('AODDUST01' , 1, ' ')
434 1536 : call add_default ('AODDUST03' , 1, ' ')
435 1536 : call add_default ('AODDUST' , 1, ' ')
436 1536 : call add_default ('AODVIS' , 1, ' ')
437 : end if
438 :
439 1536 : if (history_dust) then
440 0 : call add_default ('AODDUST01' , 1, ' ')
441 0 : call add_default ('AODDUST02' , 1, ' ')
442 0 : call add_default ('AODDUST03' , 1, ' ')
443 : end if
444 :
445 1536 : if (history_aero_optics) then
446 0 : call add_default ('AODDUST01' , 1, ' ')
447 0 : call add_default ('AODDUST03' , 1, ' ')
448 0 : call add_default ('ABSORB' , 1, ' ')
449 0 : call add_default ('AODVIS' , 1, ' ')
450 0 : call add_default ('AODUV' , 1, ' ')
451 0 : call add_default ('AODNIR' , 1, ' ')
452 0 : call add_default ('AODABS' , 1, ' ')
453 0 : call add_default ('AODABSBC' , 1, ' ')
454 0 : call add_default ('AODDUST' , 1, ' ')
455 0 : call add_default ('AODSO4' , 1, ' ')
456 0 : call add_default ('AODPOM' , 1, ' ')
457 0 : call add_default ('AODSOA' , 1, ' ')
458 0 : call add_default ('AODBC' , 1, ' ')
459 0 : call add_default ('AODSS' , 1, ' ')
460 0 : call add_default ('BURDEN01' , 1, ' ')
461 0 : call add_default ('BURDEN02' , 1, ' ')
462 0 : call add_default ('BURDEN03' , 1, ' ')
463 0 : call add_default ('BURDENDUST' , 1, ' ')
464 0 : call add_default ('BURDENSO4' , 1, ' ')
465 0 : call add_default ('BURDENPOM' , 1, ' ')
466 0 : call add_default ('BURDENSOA' , 1, ' ')
467 0 : call add_default ('BURDENBC' , 1, ' ')
468 0 : call add_default ('BURDENSEASALT', 1, ' ')
469 0 : call add_default ('SSAVIS' , 1, ' ')
470 0 : call add_default ('EXTINCT' , 1, ' ')
471 0 : call add_default ('AODxASYM' , 1, ' ')
472 0 : call add_default ('EXTxASYM' , 1, ' ')
473 :
474 0 : call add_default ('AODdnDUST01' , 1, ' ')
475 0 : call add_default ('AODdnDUST03' , 1, ' ')
476 0 : call add_default ('ABSORBdn' , 1, ' ')
477 0 : call add_default ('AODVISdn' , 1, ' ')
478 0 : call add_default ('AODUVdn' , 1, ' ')
479 0 : call add_default ('AODNIRdn' , 1, ' ')
480 0 : call add_default ('AODABSdn' , 1, ' ')
481 0 : call add_default ('AODABSBCdn' , 1, ' ')
482 0 : call add_default ('AODDUSTdn' , 1, ' ')
483 0 : call add_default ('AODSO4dn' , 1, ' ')
484 0 : call add_default ('AODPOMdn' , 1, ' ')
485 0 : call add_default ('AODSOAdn' , 1, ' ')
486 0 : call add_default ('AODBCdn' , 1, ' ')
487 0 : call add_default ('AODSSdn' , 1, ' ')
488 0 : call add_default ('BURDENdn01' , 1, ' ')
489 0 : call add_default ('BURDENdn02' , 1, ' ')
490 0 : call add_default ('BURDENdn03' , 1, ' ')
491 0 : call add_default ('BURDENDUSTdn' , 1, ' ')
492 0 : call add_default ('BURDENSO4dn' , 1, ' ')
493 0 : call add_default ('BURDENPOMdn' , 1, ' ')
494 0 : call add_default ('BURDENSOAdn' , 1, ' ')
495 0 : call add_default ('BURDENBCdn' , 1, ' ')
496 0 : call add_default ('BURDENSEASALTdn', 1, ' ')
497 0 : call add_default ('SSAVISdn' , 1, ' ')
498 0 : call add_default ('EXTINCTdn' , 1, ' ')
499 0 : call add_default ('AODxASYMdn' , 1, ' ')
500 0 : call add_default ('EXTxASYMdn' , 1, ' ')
501 : end if
502 :
503 1536 : end subroutine aerosol_optics_cam_init
504 :
505 : !===============================================================================
506 0 : subroutine aerosol_optics_cam_final
507 :
508 : integer :: iaermod
509 :
510 0 : do iaermod = 1,num_aero_models
511 0 : if (associated(aero_props(iaermod)%obj)) then
512 0 : deallocate(aero_props(iaermod)%obj)
513 0 : nullify(aero_props(iaermod)%obj)
514 : end if
515 : end do
516 :
517 0 : if (allocated(aero_props)) then
518 0 : deallocate(aero_props)
519 : endif
520 :
521 1536 : end subroutine aerosol_optics_cam_final
522 :
523 : !===============================================================================
524 30960 : subroutine aerosol_optics_cam_sw(list_idx, state, pbuf, nnite, idxnite, tauxar, wa, ga, fa)
525 :
526 : ! calculates aerosol sw radiative properties
527 :
528 : integer, intent(in) :: list_idx ! index of the climate or a diagnostic list
529 : type(physics_state), intent(in), target :: state ! state variables
530 :
531 : type(physics_buffer_desc), pointer :: pbuf(:)
532 : integer, intent(in) :: nnite ! number of night columns
533 : integer, intent(in) :: idxnite(nnite) ! local column indices of night columns
534 :
535 : real(r8), intent(inout) :: tauxar(pcols,0:pver,nswbands) ! layer extinction optical depth
536 : real(r8), intent(inout) :: wa(pcols,0:pver,nswbands) ! layer single-scatter albedo
537 : real(r8), intent(inout) :: ga(pcols,0:pver,nswbands) ! asymmetry factor
538 : real(r8), intent(inout) :: fa(pcols,0:pver,nswbands) ! forward scattered fraction
539 :
540 : character(len=*), parameter :: prefix = 'aerosol_optics_cam_sw: '
541 :
542 : integer :: ibin, nbins
543 : integer :: iwav, ilev
544 : integer :: icol, istat
545 : integer :: lchnk, ncol
546 :
547 30960 : type(aero_state_t), allocatable :: aero_state(:) ! array of aerosol state objects to allow for
548 : ! multiple aerosol representations in the same sim
549 : ! such as MAM and CARMA
550 :
551 : class(aerosol_optics), pointer :: aero_optics
552 :
553 : real(r8) :: dopaer(pcols)
554 : real(r8) :: mass(pcols,pver)
555 : real(r8) :: air_density(pcols,pver)
556 :
557 30960 : real(r8), allocatable :: pext(:) ! parameterized specific extinction (m2/kg)
558 30960 : real(r8), allocatable :: pabs(:) ! parameterized specific absorption (m2/kg)
559 30960 : real(r8), allocatable :: palb(:) ! parameterized single scattering albedo
560 30960 : real(r8), allocatable :: pasm(:) ! parameterized asymmetry factor
561 :
562 : character(len=ot_length) :: opticstype
563 : integer :: iaermod
564 :
565 : real(r8) :: aodvis(pcols) ! extinction optical depth in vis
566 : real(r8) :: aoduv(pcols) ! extinction optical depth in uv
567 : real(r8) :: aodnir(pcols) ! extinction optical depth in nir
568 : real(r8) :: absorb(pcols,pver)
569 : real(r8) :: aodabs(pcols) ! absorption optical depth
570 :
571 : real(r8) :: aodabsbc(pcols) ! absorption optical depth of BC
572 :
573 : real(r8) :: aodtot(pcols)
574 :
575 : real(r8) :: extinct(pcols,pver)
576 : real(r8) :: extinctnir(pcols,pver)
577 : real(r8) :: extinctuv(pcols,pver)
578 :
579 : real(r8) :: asymvis(pcols) ! asymmetry factor * optical depth
580 : real(r8) :: asymext(pcols,pver) ! asymmetry factor * extinction
581 :
582 : real(r8) :: wetvol(pcols,pver)
583 : real(r8) :: watervol(pcols,pver)
584 :
585 : real(r8) :: vol(pcols)
586 : real(r8) :: dustvol(pcols)
587 :
588 : real(r8) :: scatdust(pcols)
589 : real(r8) :: absdust(pcols)
590 : real(r8) :: dustaodbin(pcols)
591 :
592 : real(r8) :: scatbc(pcols)
593 : real(r8) :: absbc(pcols)
594 :
595 : real(r8) :: scatpom(pcols)
596 : real(r8) :: abspom(pcols)
597 :
598 : real(r8) :: scatsslt(pcols)
599 : real(r8) :: abssslt(pcols)
600 :
601 : real(r8) :: scatsoa(pcols)
602 : real(r8) :: abssoa(pcols)
603 :
604 : real(r8) :: scatsulf(pcols)
605 : real(r8) :: abssulf(pcols)
606 :
607 : real(r8) :: burden(pcols)
608 : real(r8) :: burdendust(pcols), burdenso4(pcols), burdenbc(pcols), &
609 : burdenpom(pcols), burdensoa(pcols), burdenseasalt(pcols)
610 :
611 : real(r8) :: hygrodust(pcols), hygrosulf(pcols), hygrobc(pcols), &
612 : hygropom(pcols), hygrosoa(pcols), hygrosslt(pcols)
613 :
614 : real(r8) :: aodbin(pcols)
615 :
616 30960 : complex(r8), pointer :: specrefindex(:) ! species refractive index
617 :
618 : class(aerosol_state), pointer :: aerostate
619 : class(aerosol_properties), pointer :: aeroprops
620 :
621 : real(r8) :: specdens
622 : character(len=32) :: spectype ! species type
623 30960 : real(r8), pointer :: specmmr(:,:)
624 : real(r8) :: hygro_aer !
625 :
626 : real(r8) :: scath2o, absh2o, sumscat, sumabs, sumhygro
627 :
628 : real(r8) :: aodc ! aod of component
629 :
630 : ! total species AOD
631 : real(r8) :: dustaod(pcols), sulfaod(pcols), bcaod(pcols), &
632 : pomaod(pcols), soaaod(pcols), ssltaod(pcols)
633 :
634 : real(r8) :: aodvisst(pcols) ! stratospheric extinction optical depth
635 : real(r8) :: aoduvst(pcols) ! stratospheric extinction optical depth in uv
636 : real(r8) :: aodnirst(pcols) ! stratospheric extinction optical depth in nir
637 : real(r8) :: ssavis(pcols)
638 : integer :: troplev(pcols)
639 :
640 30960 : nullify(aero_optics)
641 :
642 30960 : call tropopause_findChemTrop(state, troplev)
643 :
644 30960 : lchnk = state%lchnk
645 30960 : ncol = state%ncol
646 :
647 48108240 : mass(:ncol,:) = state%pdeldry(:ncol,:)*rga
648 48108240 : air_density(:ncol,:) = state%pmid(:ncol,:)/(rair*state%t(:ncol,:))
649 :
650 30960 : aodvis = 0._r8
651 30960 : aodnir = 0._r8
652 30960 : aoduv = 0._r8
653 30960 : aodabs = 0._r8
654 30960 : absorb = 0._r8
655 30960 : aodtot = 0._r8
656 30960 : tauxar = 0._r8
657 30960 : extinct = 0._r8
658 30960 : extinctnir = 0._r8
659 30960 : extinctuv = 0._r8
660 30960 : asymvis = 0.0_r8
661 30960 : asymext = 0.0_r8
662 30960 : ssavis = 0.0_r8
663 30960 : aodvisst = 0.0_r8
664 30960 : aoduvst = 0.0_r8
665 30960 : aodnirst = 0.0_r8
666 :
667 30960 : burdendust = 0.0_r8
668 30960 : burdenso4 = 0.0_r8
669 30960 : burdenbc = 0.0_r8
670 30960 : burdenpom = 0.0_r8
671 30960 : burdensoa = 0.0_r8
672 30960 : burdenseasalt = 0.0_r8
673 :
674 30960 : aodabsbc = 0.0_r8
675 30960 : dustaod = 0.0_r8
676 30960 : sulfaod = 0.0_r8
677 30960 : pomaod = 0.0_r8
678 30960 : soaaod = 0.0_r8
679 30960 : bcaod = 0.0_r8
680 30960 : ssltaod = 0.0_r8
681 :
682 30960 : if (num_aero_models<1) return
683 :
684 123840 : allocate(aero_state(num_aero_models), stat=istat)
685 30960 : if (istat/=0) then
686 0 : call endrun(prefix//'array allocation error: aero_state')
687 : end if
688 :
689 30960 : iaermod = 0
690 30960 : if (modal_active) then
691 30960 : iaermod = iaermod+1
692 30960 : aero_state(iaermod)%obj => modal_aerosol_state( state, pbuf )
693 : end if
694 :
695 92880 : allocate(pext(ncol), stat=istat)
696 30960 : if (istat/=0) then
697 0 : call endrun(prefix//'array allocation error: pext')
698 : end if
699 92880 : allocate(pabs(ncol), stat=istat)
700 30960 : if (istat/=0) then
701 0 : call endrun(prefix//'array allocation error: pabs')
702 : end if
703 92880 : allocate(palb(ncol), stat=istat)
704 30960 : if (istat/=0) then
705 0 : call endrun(prefix//'array allocation error: palb')
706 : end if
707 92880 : allocate(pasm(ncol), stat=istat)
708 30960 : if (istat/=0) then
709 0 : call endrun(prefix//'array allocation error: pasm')
710 : end if
711 :
712 61920 : aeromodel: do iaermod = 1,num_aero_models
713 :
714 30960 : aeroprops => aero_props(iaermod)%obj
715 30960 : aerostate => aero_state(iaermod)%obj
716 :
717 30960 : nbins=aeroprops%nbins(list_idx)
718 :
719 185760 : binloop: do ibin = 1, nbins
720 :
721 123840 : dustaodbin(:) = 0._r8
722 123840 : burden(:) = 0._r8
723 123840 : aodbin(:) = 0.0_r8
724 :
725 123840 : call aeroprops%optics_params(list_idx, ibin, opticstype=opticstype)
726 :
727 247680 : select case (trim(opticstype))
728 : case('modal') ! refractive method
729 : aero_optics=>refractive_aerosol_optics(aeroprops, aerostate, list_idx, ibin, &
730 123840 : ncol, pver, nswbands, nlwbands, crefwsw, crefwlw)
731 : case default
732 247680 : call endrun(prefix//'optics method not recognized')
733 : end select
734 :
735 123840 : if (associated(aero_optics)) then
736 :
737 192432960 : wetvol(:ncol,:pver) = aerostate%wet_volume(aeroprops, list_idx, ibin, ncol, pver)
738 192432960 : watervol(:ncol,:pver) = aerostate%water_volume(aeroprops, list_idx, ibin, ncol, pver)
739 :
740 1857600 : wavelength: do iwav = 1, nswbands
741 :
742 163097280 : vertical: do ilev = 1, pver
743 :
744 161239680 : call aero_optics%sw_props(ncol, ilev, iwav, pext, pabs, palb, pasm )
745 :
746 161239680 : call init_diags
747 :
748 2694061440 : column: do icol = 1,ncol
749 2531088000 : dopaer(icol) = pext(icol)*mass(icol,ilev)
750 2531088000 : tauxar(icol,ilev,iwav) = tauxar(icol,ilev,iwav) + dopaer(icol)
751 2531088000 : wa(icol,ilev,iwav) = wa(icol,ilev,iwav) + dopaer(icol)*palb(icol)
752 2531088000 : ga(icol,ilev,iwav) = ga(icol,ilev,iwav) + dopaer(icol)*palb(icol)*pasm(icol)
753 2531088000 : fa(icol,ilev,iwav) = fa(icol,ilev,iwav) + dopaer(icol)*palb(icol)*pasm(icol)*pasm(icol)
754 :
755 2692327680 : call update_diags
756 :
757 : end do column
758 :
759 : end do vertical
760 : end do wavelength
761 :
762 : else
763 0 : call endrun(prefix//'aero_optics object pointer not associated')
764 : end if
765 :
766 123840 : deallocate(aero_optics)
767 : nullify(aero_optics)
768 :
769 154800 : call output_bin_diags
770 :
771 : end do binloop
772 : end do aeromodel
773 :
774 30960 : call output_tot_diags
775 :
776 30960 : deallocate(pext)
777 30960 : deallocate(pabs)
778 30960 : deallocate(palb)
779 30960 : deallocate(pasm)
780 :
781 61920 : do iaermod = 1,num_aero_models
782 61920 : deallocate(aero_state(iaermod)%obj)
783 61920 : nullify(aero_state(iaermod)%obj)
784 : end do
785 :
786 30960 : deallocate(aero_state)
787 :
788 : contains
789 :
790 : !===============================================================================
791 161239680 : subroutine init_diags
792 2853567360 : dustvol(:ncol) = 0._r8
793 2692327680 : scatdust(:ncol) = 0._r8
794 2692327680 : absdust(:ncol) = 0._r8
795 2692327680 : hygrodust(:ncol) = 0._r8
796 2692327680 : scatsulf(:ncol) = 0._r8
797 2692327680 : abssulf(:ncol) = 0._r8
798 2692327680 : hygrosulf(:ncol) = 0._r8
799 2692327680 : scatbc(:ncol) = 0._r8
800 2692327680 : absbc(:ncol) = 0._r8
801 2692327680 : hygrobc(:ncol) = 0._r8
802 2692327680 : scatpom(:ncol) = 0._r8
803 2692327680 : abspom(:ncol) = 0._r8
804 2692327680 : hygropom(:ncol) = 0._r8
805 2692327680 : scatsoa(:ncol) = 0._r8
806 2692327680 : abssoa(:ncol) = 0._r8
807 2692327680 : hygrosoa(:ncol) = 0._r8
808 2692327680 : scatsslt(:ncol) = 0._r8
809 2692327680 : abssslt(:ncol) = 0._r8
810 2692327680 : hygrosslt(:ncol) = 0._r8
811 161239680 : end subroutine init_diags
812 :
813 : !===============================================================================
814 2531088000 : subroutine update_diags
815 :
816 : integer :: ispec
817 :
818 2531088000 : if (iwav==idx_uv_diag) then
819 180792000 : aoduv(icol) = aoduv(icol) + dopaer(icol)
820 180792000 : extinctuv(icol,ilev) = extinctuv(icol,ilev) + dopaer(icol)*air_density(icol,ilev)/mass(icol,ilev)
821 180792000 : if (ilev<=troplev(icol)) then
822 101232476 : aoduvst(icol) = aoduvst(icol) + dopaer(icol)
823 : end if
824 :
825 2350296000 : else if (iwav==idx_sw_diag) then ! vis
826 180792000 : aodvis(icol) = aodvis(icol) + dopaer(icol)
827 180792000 : aodabs(icol) = aodabs(icol) + pabs(icol)*mass(icol,ilev)
828 180792000 : extinct(icol,ilev) = extinct(icol,ilev) + dopaer(icol)*air_density(icol,ilev)/mass(icol,ilev)
829 180792000 : absorb(icol,ilev) = absorb(icol,ilev) + pabs(icol)*air_density(icol,ilev)
830 180792000 : ssavis(icol) = ssavis(icol) + dopaer(icol)*palb(icol)
831 180792000 : asymvis(icol) = asymvis(icol) + dopaer(icol)*pasm(icol)
832 180792000 : asymext(icol,ilev) = asymext(icol,ilev) + dopaer(icol)*pasm(icol)*air_density(icol,ilev)/mass(icol,ilev)
833 :
834 180792000 : aodbin(icol) = aodbin(icol) + dopaer(icol)
835 :
836 180792000 : if (ilev<=troplev(icol)) then
837 101232476 : aodvisst(icol) = aodvisst(icol) + dopaer(icol)
838 : end if
839 :
840 : ! loop over species ...
841 :
842 858762000 : do ispec = 1, aeroprops%nspecies(list_idx,ibin)
843 : call aeroprops%get(ibin, ispec, list_ndx=list_idx, density=specdens, &
844 677970000 : spectype=spectype, refindex_sw=specrefindex, hygro=hygro_aer)
845 677970000 : call aerostate%get_ambient_mmr(list_idx, ispec, ibin, specmmr)
846 :
847 677970000 : burden(icol) = burden(icol) + specmmr(icol,ilev)*mass(icol,ilev)
848 :
849 677970000 : vol(icol) = specmmr(icol,ilev)/specdens
850 :
851 1536732000 : select case ( trim(spectype) )
852 : case('dust')
853 135594000 : dustvol(icol) = vol(icol)
854 135594000 : burdendust(icol) = burdendust(icol) + specmmr(icol,ilev)*mass(icol,ilev)
855 135594000 : scatdust(icol) = vol(icol) * specrefindex(iwav)%re
856 135594000 : absdust(icol) =-vol(icol) * specrefindex(iwav)%im
857 135594000 : hygrodust(icol)= vol(icol)*hygro_aer
858 : case('black-c')
859 90396000 : burdenbc(icol) = burdenbc(icol) + specmmr(icol,ilev)*mass(icol,ilev)
860 90396000 : scatbc(icol) = vol(icol) * specrefindex(iwav)%re
861 90396000 : absbc(icol) =-vol(icol) * specrefindex(iwav)%im
862 90396000 : hygrobc(icol)= vol(icol)*hygro_aer
863 : case('sulfate')
864 135594000 : burdenso4(icol) = burdenso4(icol) + specmmr(icol,ilev)*mass(icol,ilev)
865 135594000 : scatsulf(icol) = vol(icol) * specrefindex(iwav)%re
866 135594000 : abssulf(icol) =-vol(icol) * specrefindex(iwav)%im
867 135594000 : hygrosulf(icol)= vol(icol)*hygro_aer
868 : case('p-organic')
869 90396000 : burdenpom(icol) = burdenpom(icol) + specmmr(icol,ilev)*mass(icol,ilev)
870 90396000 : scatpom(icol) = vol(icol) * specrefindex(iwav)%re
871 90396000 : abspom(icol) =-vol(icol) * specrefindex(iwav)%im
872 90396000 : hygropom(icol)= vol(icol)*hygro_aer
873 : case('s-organic')
874 90396000 : burdensoa(icol) = burdensoa(icol) + specmmr(icol,ilev)*mass(icol,ilev)
875 90396000 : scatsoa(icol) = vol(icol) * specrefindex(iwav)%re
876 90396000 : abssoa(icol) = -vol(icol) * specrefindex(iwav)%im
877 90396000 : hygrosoa(icol)= vol(icol)*hygro_aer
878 : case('seasalt')
879 135594000 : burdenseasalt(icol) = burdenseasalt(icol) + specmmr(icol,ilev)*mass(icol,ilev)
880 135594000 : scatsslt(icol) = vol(icol) * specrefindex(iwav)%re
881 135594000 : abssslt(icol) = -vol(icol) * specrefindex(iwav)%im
882 1491534000 : hygrosslt(icol)= vol(icol)*hygro_aer
883 : end select
884 : end do
885 :
886 180792000 : if (wetvol(icol,ilev)>1.e-40_r8 .and. vol(icol)>0._r8) then
887 :
888 180792000 : dustaodbin(icol) = dustaodbin(icol) + dopaer(icol)*dustvol(icol)/wetvol(icol,ilev)
889 :
890 : ! partition optical depth into contributions from each constituent
891 : ! assume contribution is proportional to refractive index X volume
892 :
893 180792000 : scath2o = watervol(icol,ilev)*crefwsw(iwav)%re
894 180792000 : absh2o = -watervol(icol,ilev)*crefwsw(iwav)%im
895 : sumscat = scatsulf(icol) + scatpom(icol) + scatsoa(icol) + scatbc(icol) + &
896 180792000 : scatdust(icol) + scatsslt(icol) + scath2o
897 : sumabs = abssulf(icol) + abspom(icol) + abssoa(icol) + absbc(icol) + &
898 180792000 : absdust(icol) + abssslt(icol) + absh2o
899 : sumhygro = hygrosulf(icol) + hygropom(icol) + hygrosoa(icol) + hygrobc(icol) + &
900 180792000 : hygrodust(icol) + hygrosslt(icol)
901 :
902 180792000 : scatdust(icol) = (scatdust(icol) + scath2o*hygrodust(icol)/sumhygro)/sumscat
903 180792000 : absdust(icol) = (absdust(icol) + absh2o*hygrodust(icol)/sumhygro)/sumabs
904 :
905 180792000 : scatsulf(icol) = (scatsulf(icol) + scath2o*hygrosulf(icol)/sumhygro)/sumscat
906 180792000 : abssulf(icol) = (abssulf(icol) + absh2o*hygrosulf(icol)/sumhygro)/sumabs
907 :
908 180792000 : scatpom(icol) = (scatpom(icol) + scath2o*hygropom(icol)/sumhygro)/sumscat
909 180792000 : abspom(icol) = (abspom(icol) + absh2o*hygropom(icol)/sumhygro)/sumabs
910 :
911 180792000 : scatsoa(icol) = (scatsoa(icol) + scath2o*hygrosoa(icol)/sumhygro)/sumscat
912 180792000 : abssoa(icol) = (abssoa(icol) + absh2o*hygrosoa(icol)/sumhygro)/sumabs
913 :
914 180792000 : scatbc(icol)= (scatbc(icol) + scath2o*hygrobc(icol)/sumhygro)/sumscat
915 180792000 : absbc(icol) = (absbc(icol) + absh2o*hygrobc(icol)/sumhygro)/sumabs
916 :
917 180792000 : scatsslt(icol) = (scatsslt(icol) + scath2o*hygrosslt(icol)/sumhygro)/sumscat
918 180792000 : abssslt(icol) = (abssslt(icol) + absh2o*hygrosslt(icol)/sumhygro)/sumabs
919 :
920 :
921 180792000 : aodabsbc(icol) = aodabsbc(icol) + absbc(icol)*dopaer(icol)*(1.0_r8-palb(icol))
922 :
923 :
924 :
925 180792000 : aodc = (absdust(icol)*(1.0_r8 - palb(icol)) + palb(icol)*scatdust(icol))*dopaer(icol)
926 180792000 : dustaod(icol) = dustaod(icol) + aodc
927 :
928 180792000 : aodc = (abssulf(icol)*(1.0_r8 - palb(icol)) + palb(icol)*scatsulf(icol))*dopaer(icol)
929 180792000 : sulfaod(icol) = sulfaod(icol) + aodc
930 :
931 180792000 : aodc = (abspom(icol)*(1.0_r8 - palb(icol)) + palb(icol)*scatpom(icol))*dopaer(icol)
932 180792000 : pomaod(icol) = pomaod(icol) + aodc
933 :
934 180792000 : aodc = (abssoa(icol)*(1.0_r8 - palb(icol)) + palb(icol)*scatsoa(icol))*dopaer(icol)
935 180792000 : soaaod(icol) = soaaod(icol) + aodc
936 :
937 180792000 : aodc = (absbc(icol)*(1.0_r8 - palb(icol)) + palb(icol)*scatbc(icol))*dopaer(icol)
938 180792000 : bcaod(icol) = bcaod(icol) + aodc
939 :
940 180792000 : aodc = (abssslt(icol)*(1.0_r8 - palb(icol)) + palb(icol)*scatsslt(icol))*dopaer(icol)
941 180792000 : ssltaod(icol) = ssltaod(icol) + aodc
942 :
943 : end if
944 2169504000 : else if (iwav==idx_nir_diag) then
945 180792000 : aodnir(icol) = aodnir(icol) + dopaer(icol)
946 180792000 : extinctnir(icol,ilev) = extinctnir(icol,ilev) + dopaer(icol)*air_density(icol,ilev)/mass(icol,ilev)
947 :
948 180792000 : if (ilev<=troplev(icol)) then
949 101232476 : aodnirst(icol) = aodnirst(icol) + dopaer(icol)
950 : end if
951 :
952 : end if
953 :
954 2531088000 : aodtot(icol) = aodtot(icol) + dopaer(icol)
955 :
956 2531088000 : end subroutine update_diags
957 :
958 : !===============================================================================
959 123840 : subroutine output_bin_diags
960 :
961 : integer :: icol
962 :
963 123840 : if (list_idx == 0) then
964 :
965 123840 : call outfld(burdendn_fields(iaermod)%name(ibin), burden, pcols, lchnk)
966 123840 : call outfld(aoddustdn_fields(iaermod)%name(ibin), dustaodbin, pcols, lchnk)
967 123840 : call outfld(aodbindn_fields(iaermod)%name(ibin), aodbin, pcols, lchnk)
968 :
969 1095840 : do icol = 1, nnite
970 972000 : burden(idxnite(icol)) = fillvalue
971 972000 : aodbin(idxnite(icol)) = fillvalue
972 1095840 : dustaodbin(idxnite(icol)) = fillvalue
973 : end do
974 :
975 123840 : call outfld(burden_fields(iaermod)%name(ibin), burden, pcols, lchnk)
976 123840 : call outfld(aoddust_fields(iaermod)%name(ibin), dustaodbin, pcols, lchnk)
977 123840 : call outfld(aodbin_fields(iaermod)%name(ibin), aodbin, pcols, lchnk)
978 :
979 : endif
980 :
981 123840 : end subroutine output_bin_diags
982 :
983 : !===============================================================================
984 30960 : subroutine output_tot_diags
985 :
986 : integer :: icol
987 :
988 30960 : call outfld('AODUVdn'//diag(list_idx), aoduv, pcols, lchnk)
989 30960 : call outfld('AODVISdn'//diag(list_idx), aodvis, pcols, lchnk)
990 30960 : call outfld('AODABSdn'//diag(list_idx), aodabs, pcols, lchnk)
991 30960 : call outfld('AODNIRdn'//diag(list_idx), aodnir, pcols, lchnk)
992 30960 : call outfld('AODTOTdn'//diag(list_idx), aodtot, pcols, lchnk)
993 30960 : call outfld('EXTINCTUVdn'//diag(list_idx), extinctuv, pcols, lchnk)
994 30960 : call outfld('EXTINCTNIRdn'//diag(list_idx), extinctnir, pcols, lchnk)
995 30960 : call outfld('EXTINCTdn'//diag(list_idx), extinct, pcols, lchnk)
996 30960 : call outfld('ABSORBdn'//diag(list_idx), absorb, pcols, lchnk)
997 30960 : call outfld('EXTxASYMdn'//diag(list_idx), asymext, pcols, lchnk)
998 30960 : call outfld('AODxASYMdn'//diag(list_idx), asymvis, pcols, lchnk)
999 :
1000 30960 : call outfld('AODVISstdn'//diag(list_idx), aodvisst,pcols, lchnk)
1001 30960 : call outfld('AODUVstdn'//diag(list_idx), aoduvst, pcols, lchnk)
1002 30960 : call outfld('AODNIRstdn'//diag(list_idx), aodnirst,pcols, lchnk)
1003 :
1004 273960 : do icol = 1, nnite
1005 243000 : aodvis(idxnite(icol)) = fillvalue
1006 243000 : aodnir(idxnite(icol)) = fillvalue
1007 243000 : aoduv(idxnite(icol)) = fillvalue
1008 243000 : aodabs(idxnite(icol)) = fillvalue
1009 243000 : aodtot(idxnite(icol)) = fillvalue
1010 22842000 : extinct(idxnite(icol),:) = fillvalue
1011 22842000 : extinctnir(idxnite(icol),:) = fillvalue
1012 22842000 : extinctuv(idxnite(icol),:) = fillvalue
1013 22842000 : absorb(idxnite(icol),:) = fillvalue
1014 22842000 : asymext(idxnite(icol),:) = fillvalue
1015 243000 : asymvis(idxnite(icol)) = fillvalue
1016 243000 : aodabs(idxnite(icol)) = fillvalue
1017 243000 : aodvisst(idxnite(icol)) = fillvalue
1018 243000 : aoduvst(idxnite(icol)) = fillvalue
1019 273960 : aodnirst(idxnite(icol)) = fillvalue
1020 : end do
1021 :
1022 30960 : call outfld('AODUV'//diag(list_idx), aoduv, pcols, lchnk)
1023 30960 : call outfld('AODVIS'//diag(list_idx), aodvis, pcols, lchnk)
1024 30960 : call outfld('AODABS'//diag(list_idx), aodabs, pcols, lchnk)
1025 30960 : call outfld('AODNIR'//diag(list_idx), aodnir, pcols, lchnk)
1026 30960 : call outfld('AODTOT'//diag(list_idx), aodtot, pcols, lchnk)
1027 30960 : call outfld('EXTINCTUV'//diag(list_idx), extinctuv, pcols, lchnk)
1028 30960 : call outfld('EXTINCTNIR'//diag(list_idx), extinctnir, pcols, lchnk)
1029 30960 : call outfld('EXTINCT'//diag(list_idx), extinct, pcols, lchnk)
1030 30960 : call outfld('ABSORB'//diag(list_idx), absorb, pcols, lchnk)
1031 30960 : call outfld('EXTxASYM'//diag(list_idx), asymext, pcols, lchnk)
1032 30960 : call outfld('AODxASYM'//diag(list_idx), asymvis, pcols, lchnk)
1033 30960 : call outfld('AODVISst'//diag(list_idx), aodvisst,pcols, lchnk)
1034 30960 : call outfld('AODUVst'//diag(list_idx), aoduvst, pcols, lchnk)
1035 30960 : call outfld('AODNIRst'//diag(list_idx), aodnirst,pcols, lchnk)
1036 :
1037 : ! These diagnostics are output only for climate list
1038 30960 : if (list_idx == 0) then
1039 516960 : do icol = 1, ncol
1040 516960 : if (aodvis(icol) > 1.e-10_r8) then
1041 486000 : ssavis(icol) = ssavis(icol)/aodvis(icol)
1042 : else
1043 0 : ssavis(icol) = 0.925_r8
1044 : endif
1045 : end do
1046 30960 : call outfld('SSAVISdn', ssavis, pcols, lchnk)
1047 :
1048 30960 : call outfld('BURDENDUSTdn', burdendust, pcols, lchnk)
1049 30960 : call outfld('BURDENSO4dn' , burdenso4, pcols, lchnk)
1050 30960 : call outfld('BURDENPOMdn' , burdenpom, pcols, lchnk)
1051 30960 : call outfld('BURDENSOAdn' , burdensoa, pcols, lchnk)
1052 30960 : call outfld('BURDENBCdn' , burdenbc, pcols, lchnk)
1053 30960 : call outfld('BURDENSEASALTdn', burdenseasalt, pcols, lchnk)
1054 :
1055 30960 : call outfld('AODABSBCdn', aodabsbc, pcols, lchnk)
1056 :
1057 30960 : call outfld('AODDUSTdn', dustaod, pcols, lchnk)
1058 30960 : call outfld('AODSO4dn', sulfaod, pcols, lchnk)
1059 30960 : call outfld('AODPOMdn', pomaod, pcols, lchnk)
1060 30960 : call outfld('AODSOAdn', soaaod, pcols, lchnk)
1061 30960 : call outfld('AODBCdn', bcaod, pcols, lchnk)
1062 30960 : call outfld('AODSSdn', ssltaod, pcols, lchnk)
1063 :
1064 :
1065 273960 : do icol = 1, nnite
1066 :
1067 243000 : ssavis(idxnite(icol)) = fillvalue
1068 243000 : asymvis(idxnite(icol)) = fillvalue
1069 :
1070 243000 : burdendust(idxnite(icol)) = fillvalue
1071 243000 : burdenso4(idxnite(icol)) = fillvalue
1072 243000 : burdenpom(idxnite(icol)) = fillvalue
1073 243000 : burdensoa(idxnite(icol)) = fillvalue
1074 243000 : burdenbc(idxnite(icol)) = fillvalue
1075 243000 : burdenseasalt(idxnite(icol)) = fillvalue
1076 243000 : aodabsbc(idxnite(icol)) = fillvalue
1077 :
1078 243000 : dustaod(idxnite(icol)) = fillvalue
1079 243000 : sulfaod(idxnite(icol)) = fillvalue
1080 243000 : pomaod(idxnite(icol)) = fillvalue
1081 243000 : soaaod(idxnite(icol)) = fillvalue
1082 243000 : bcaod(idxnite(icol)) = fillvalue
1083 273960 : ssltaod(idxnite(icol)) = fillvalue
1084 :
1085 : end do
1086 :
1087 30960 : call outfld('SSAVIS', ssavis, pcols, lchnk)
1088 30960 : call outfld('AODxASYM', asymvis, pcols, lchnk)
1089 30960 : call outfld('BURDENDUST', burdendust, pcols, lchnk)
1090 30960 : call outfld('BURDENSO4' , burdenso4, pcols, lchnk)
1091 30960 : call outfld('BURDENPOM' , burdenpom, pcols, lchnk)
1092 30960 : call outfld('BURDENSOA' , burdensoa, pcols, lchnk)
1093 30960 : call outfld('BURDENBC' , burdenbc, pcols, lchnk)
1094 30960 : call outfld('BURDENSEASALT', burdenseasalt, pcols, lchnk)
1095 30960 : call outfld('AODABSBC', aodabsbc, pcols, lchnk)
1096 30960 : call outfld('AODDUST', dustaod, pcols, lchnk)
1097 30960 : call outfld('AODSO4', sulfaod, pcols, lchnk)
1098 30960 : call outfld('AODPOM', pomaod, pcols, lchnk)
1099 30960 : call outfld('AODSOA', soaaod, pcols, lchnk)
1100 30960 : call outfld('AODBC', bcaod, pcols, lchnk)
1101 30960 : call outfld('AODSS', ssltaod, pcols, lchnk)
1102 :
1103 : end if
1104 :
1105 30960 : end subroutine output_tot_diags
1106 :
1107 : end subroutine aerosol_optics_cam_sw
1108 :
1109 : !===============================================================================
1110 30960 : subroutine aerosol_optics_cam_lw(list_idx, state, pbuf, tauxar)
1111 :
1112 : ! calculates aerosol lw radiative properties
1113 :
1114 : integer, intent(in) :: list_idx ! index of the climate or a diagnostic list
1115 : type(physics_state), intent(in), target :: state ! state variables
1116 :
1117 : type(physics_buffer_desc), pointer :: pbuf(:)
1118 :
1119 : real(r8), intent(inout) :: tauxar(pcols,pver,nlwbands) ! layer absorption optical depth
1120 :
1121 :
1122 : real(r8) :: dopaer(pcols)
1123 : real(r8) :: mass(pcols,pver)
1124 :
1125 : character(len=*), parameter :: prefix = 'aerosol_optics_cam_lw: '
1126 :
1127 : integer :: ibin, nbins
1128 : integer :: iwav, ilev
1129 : integer :: ncol, icol, istat
1130 :
1131 30960 : type(aero_state_t), allocatable :: aero_state(:) ! array of aerosol state objects to allow for
1132 : ! multiple aerosol representations in the same sim
1133 : ! such as MAM and CARMA
1134 :
1135 : class(aerosol_optics), pointer :: aero_optics
1136 : class(aerosol_state), pointer :: aerostate
1137 : class(aerosol_properties), pointer :: aeroprops
1138 :
1139 30960 : real(r8), allocatable :: pabs(:)
1140 :
1141 : character(len=32) :: opticstype
1142 : integer :: iaermod
1143 :
1144 : real(r8) :: lwabs(pcols,pver)
1145 30960 : lwabs = 0._r8
1146 30960 : tauxar = 0._r8
1147 :
1148 30960 : nullify(aero_optics)
1149 :
1150 123840 : allocate(aero_state(num_aero_models), stat=istat)
1151 30960 : if (istat/=0) then
1152 0 : call endrun(prefix//'array allocation error: aero_state')
1153 : end if
1154 :
1155 30960 : iaermod = 0
1156 30960 : if (modal_active) then
1157 30960 : iaermod = iaermod+1
1158 30960 : aero_state(iaermod)%obj => modal_aerosol_state( state, pbuf )
1159 : end if
1160 :
1161 30960 : ncol = state%ncol
1162 :
1163 48108240 : mass(:ncol,:) = state%pdeldry(:ncol,:)*rga
1164 :
1165 92880 : allocate(pabs(ncol), stat=istat)
1166 30960 : if (istat/=0) then
1167 0 : call endrun(prefix//'array allocation error: pabs')
1168 : end if
1169 :
1170 61920 : aeromodel: do iaermod = 1,num_aero_models
1171 :
1172 30960 : aeroprops => aero_props(iaermod)%obj
1173 30960 : aerostate => aero_state(iaermod)%obj
1174 :
1175 30960 : nbins=aero_props(iaermod)%obj%nbins(list_idx)
1176 :
1177 185760 : binloop: do ibin = 1, nbins
1178 :
1179 123840 : call aeroprops%optics_params(list_idx, ibin, opticstype=opticstype)
1180 :
1181 247680 : select case (trim(opticstype))
1182 : case('modal') ! refractive method
1183 : aero_optics=>refractive_aerosol_optics(aeroprops, aerostate, list_idx, ibin, &
1184 123840 : ncol, pver, nswbands, nlwbands, crefwsw, crefwlw)
1185 : case default
1186 247680 : call endrun(prefix//'optics method not recognized')
1187 : end select
1188 :
1189 123840 : if (associated(aero_optics)) then
1190 :
1191 2105280 : wavelength: do iwav = 1, nlwbands
1192 :
1193 186379200 : vertical: do ilev = 1, pver
1194 184273920 : call aero_optics%lw_props(ncol, ilev, iwav, pabs )
1195 :
1196 3078927360 : column: do icol = 1, ncol
1197 2892672000 : dopaer(icol) = pabs(icol)*mass(icol,ilev)
1198 2892672000 : tauxar(icol,ilev,iwav) = tauxar(icol,ilev,iwav) + dopaer(icol)
1199 3076945920 : lwabs(icol,ilev) = lwabs(icol,ilev) + pabs(icol)
1200 : end do column
1201 :
1202 : end do vertical
1203 :
1204 : end do wavelength
1205 :
1206 : else
1207 0 : call endrun(prefix//'aero_optics object pointer not associated')
1208 : end if
1209 :
1210 247680 : deallocate(aero_optics)
1211 30960 : nullify(aero_optics)
1212 :
1213 : end do binloop
1214 : end do aeromodel
1215 :
1216 30960 : call outfld('TOTABSLW'//diag(list_idx), lwabs(:,:), pcols, state%lchnk)
1217 :
1218 30960 : if (lw10um_indx>0) then
1219 30960 : call outfld('AODABSLW'//diag(list_idx), tauxar(:,:,lw10um_indx), pcols, state%lchnk)
1220 : end if
1221 :
1222 30960 : deallocate(pabs)
1223 :
1224 61920 : do iaermod = 1,num_aero_models
1225 30960 : deallocate(aero_state(iaermod)%obj)
1226 61920 : nullify(aero_state(iaermod)%obj)
1227 : end do
1228 :
1229 30960 : deallocate(aero_state)
1230 :
1231 30960 : end subroutine aerosol_optics_cam_lw
1232 :
1233 : !===============================================================================
1234 : ! Private routines
1235 : !===============================================================================
1236 :
1237 1536 : subroutine read_water_refindex(infilename)
1238 : use cam_pio_utils, only: cam_pio_openfile
1239 : use pio, only: file_desc_t, var_desc_t, pio_inq_dimlen, pio_inq_dimid, pio_inq_varid, &
1240 : pio_get_var, PIO_NOWRITE, pio_closefile, pio_noerr
1241 :
1242 :
1243 : ! read water refractive index file and set module data
1244 :
1245 : character*(*), intent(in) :: infilename ! modal optics filename
1246 :
1247 : ! Local variables
1248 :
1249 : integer :: i, ierr
1250 : type(file_desc_t) :: ncid ! pio file handle
1251 : integer :: did ! dimension ids
1252 : integer :: dimlen ! dimension lengths
1253 : type(var_desc_t) :: vid ! variable ids
1254 : real(r8) :: refrwsw(nswbands), refiwsw(nswbands) ! real, imaginary ref index for water visible
1255 : real(r8) :: refrwlw(nlwbands), refiwlw(nlwbands) ! real, imaginary ref index for water infrared
1256 :
1257 : character(len=*), parameter :: prefix = 'read_water_refindex: '
1258 : !----------------------------------------------------------------------------
1259 :
1260 : ! open file
1261 1536 : call cam_pio_openfile(ncid, infilename, PIO_NOWRITE)
1262 :
1263 : ! inquire dimensions. Check that file values match parameter values.
1264 :
1265 1536 : ierr = pio_inq_dimid(ncid, 'lw_band', did)
1266 1536 : if (ierr /= pio_noerr ) then
1267 0 : call endrun(prefix//'pio_inq_dimid lw_band')
1268 : end if
1269 1536 : ierr = pio_inq_dimlen(ncid, did, dimlen)
1270 1536 : if (ierr /= pio_noerr ) then
1271 0 : call endrun(prefix//'pio_inq_dimlen lw_band')
1272 : end if
1273 1536 : if (dimlen /= nlwbands) then
1274 0 : write(iulog,*) 'lw_band len=', dimlen, ' from ', infilename, ' ne nlwbands=', nlwbands
1275 0 : call endrun(prefix//'bad lw_band value')
1276 : endif
1277 :
1278 1536 : ierr = pio_inq_dimid(ncid, 'sw_band', did)
1279 1536 : if (ierr /= pio_noerr ) then
1280 0 : call endrun(prefix//'pio_inq_dimid sw_band')
1281 : end if
1282 1536 : ierr = pio_inq_dimlen(ncid, did, dimlen)
1283 1536 : if (ierr /= pio_noerr ) then
1284 0 : call endrun(prefix//'pio_inq_dimlen sw_band')
1285 : end if
1286 1536 : if (dimlen /= nswbands) then
1287 0 : write(iulog,*) 'sw_band len=', dimlen, ' from ', infilename, ' ne nswbands=', nswbands
1288 0 : call endrun(prefix//'bad sw_band value')
1289 : endif
1290 :
1291 : ! read variables
1292 1536 : ierr = pio_inq_varid(ncid, 'refindex_real_water_sw', vid)
1293 1536 : if (ierr /= pio_noerr ) then
1294 0 : call endrun(prefix//'pio_inq_varid refindex_real_water_sw')
1295 : end if
1296 1536 : ierr = pio_get_var(ncid, vid, refrwsw)
1297 1536 : if (ierr /= pio_noerr ) then
1298 0 : call endrun(prefix//'pio_get_var refrwsw')
1299 : end if
1300 :
1301 1536 : ierr = pio_inq_varid(ncid, 'refindex_im_water_sw', vid)
1302 1536 : if (ierr /= pio_noerr ) then
1303 0 : call endrun(prefix//'pio_inq_varid refindex_im_water_sw')
1304 : end if
1305 1536 : ierr = pio_get_var(ncid, vid, refiwsw)
1306 1536 : if (ierr /= pio_noerr ) then
1307 0 : call endrun(prefix//'pio_get_var refiwsw')
1308 : end if
1309 :
1310 1536 : ierr = pio_inq_varid(ncid, 'refindex_real_water_lw', vid)
1311 1536 : if (ierr /= pio_noerr ) then
1312 0 : call endrun(prefix//'pio_inq_varid refindex_real_water_lw')
1313 : end if
1314 1536 : ierr = pio_get_var(ncid, vid, refrwlw)
1315 1536 : if (ierr /= pio_noerr ) then
1316 0 : call endrun(prefix//'pio_get_var refrwlw')
1317 : end if
1318 :
1319 1536 : ierr = pio_inq_varid(ncid, 'refindex_im_water_lw', vid)
1320 1536 : if (ierr /= pio_noerr ) then
1321 0 : call endrun(prefix//'pio_inq_varid refindex_im_water_lw')
1322 : end if
1323 1536 : ierr = pio_get_var(ncid, vid, refiwlw)
1324 1536 : if (ierr /= pio_noerr ) then
1325 0 : call endrun(prefix//'pio_get_var refiwlw')
1326 : end if
1327 :
1328 : ! set complex representation of refractive indices as module data
1329 23040 : do i = 1, nswbands
1330 23040 : crefwsw(i) = cmplx(refrwsw(i), abs(refiwsw(i)), kind=r8)
1331 : end do
1332 26112 : do i = 1, nlwbands
1333 26112 : crefwlw(i) = cmplx(refrwlw(i), abs(refiwlw(i)), kind=r8)
1334 : end do
1335 :
1336 1536 : call pio_closefile(ncid)
1337 :
1338 3072 : end subroutine read_water_refindex
1339 :
1340 1536 : end module aerosol_optics_cam
|
