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