Line data Source code
1 : module phys_prop
2 :
3 : ! Properties of aerosols that are used by radiation and other parameterizations.
4 :
5 : ! *****N.B.*****
6 : ! This module is a utility used by the rad_constituents module. The properties stored
7 : ! here are meant to be accessed via that module. This module knows nothing about how
8 : ! this data is associated with the constituents that are radiatively active or those that
9 : ! are being used for diagnostic calculations. That is the responsibility of the
10 : ! rad_constituents module.
11 :
12 : use shr_kind_mod, only: r8 => shr_kind_r8
13 : use spmd_utils, only: masterproc
14 : use radconstants, only: nlwbands, nswbands, idx_sw_diag
15 : use ioFileMod, only: getfil
16 : use cam_pio_utils, only: cam_pio_openfile
17 : use pio, only: file_desc_t, var_desc_t, pio_get_var, pio_inq_varid, &
18 : pio_inq_dimlen, pio_inq_dimid , pio_nowrite, pio_closefile, &
19 : pio_seterrorhandling, PIO_BCAST_ERROR, PIO_INTERNAL_ERROR, PIO_NOERR
20 :
21 : use cam_logfile, only: iulog
22 : use cam_abortutils, only: endrun
23 :
24 : implicit none
25 : private
26 : save
27 :
28 : integer, parameter, public :: ot_length = 32
29 :
30 : public :: &
31 : physprop_accum_unique_files, &! Make a list of the unique set of files that contain properties
32 : ! This is an initialization step that must be done before calling physprop_init
33 : physprop_init, &! Initialization -- read the input datasets
34 : physprop_get_id, &! Return ID used to access the property data from the input files
35 : physprop_get ! Return data for specified ID
36 :
37 : ! Data from one input dataset is stored in a structure of type(physprop_type).
38 : type :: physprop_type
39 : character(len=256) :: sourcefile ! Absolute pathname of data file.
40 : character(len=ot_length) :: opticsmethod ! one of {hygro,nonhygro}
41 :
42 : ! for hygroscopic species of externally mixed aerosols
43 : real(r8), pointer :: sw_hygro_ext(:,:)
44 : real(r8), pointer :: sw_hygro_ssa(:,:)
45 : real(r8), pointer :: sw_hygro_asm(:,:)
46 : real(r8), pointer :: lw_hygro_abs(:,:)
47 :
48 : ! for nonhygroscopic species of externally mixed aerosols
49 : real(r8), pointer :: sw_nonhygro_ext(:)
50 : real(r8), pointer :: sw_nonhygro_ssa(:)
51 : real(r8), pointer :: sw_nonhygro_asm(:)
52 : real(r8), pointer :: sw_nonhygro_scat(:)
53 : real(r8), pointer :: sw_nonhygro_ascat(:)
54 : real(r8), pointer :: lw_abs(:)
55 :
56 : ! complex refractive index
57 : complex(r8), pointer :: refindex_aer_sw(:)
58 : complex(r8), pointer :: refindex_aer_lw(:)
59 :
60 : ! for radius-dependent mass-specific quantities
61 : real(r8), pointer :: r_sw_ext(:,:)
62 : real(r8), pointer :: r_sw_scat(:,:)
63 : real(r8), pointer :: r_sw_ascat(:,:)
64 : real(r8), pointer :: r_lw_abs(:,:)
65 : real(r8), pointer :: mu(:)
66 :
67 : ! for modal optics
68 : real(r8), pointer :: extpsw(:,:,:,:) ! specific extinction
69 : real(r8), pointer :: abspsw(:,:,:,:) ! specific absorption
70 : real(r8), pointer :: asmpsw(:,:,:,:) ! asymmetry factor
71 : real(r8), pointer :: absplw(:,:,:,:) ! specific absorption
72 : real(r8), pointer :: refrtabsw(:,:) ! table of real refractive indices for aerosols visible
73 : real(r8), pointer :: refitabsw(:,:) ! table of imag refractive indices for aerosols visible
74 : real(r8), pointer :: refrtablw(:,:) ! table of real refractive indices for aerosols infrared
75 : real(r8), pointer :: refitablw(:,:) ! table of imag refractive indices for aerosols infrared
76 :
77 : ! microphysics parameters.
78 : character(len=32) :: aername ! for output of number concentration
79 : real(r8) :: density_aer ! density of aerosol (kg/m3)
80 : real(r8) :: hygro_aer ! hygroscopicity of aerosol
81 : real(r8) :: dryrad_aer ! number mode radius (m) of aerosol size distribution
82 : real(r8) :: dispersion_aer ! geometric standard deviation of aerosol size distribution
83 : real(r8) :: num_to_mass_aer ! ratio of number concentration to mass concentration (#/kg)
84 : ! *** Is this actually (kg/#) ???
85 : ! mode parameters
86 : integer :: ncoef ! number of Chebyshev coefficients
87 : integer :: prefr ! dimension in table of real refractive indices
88 : integer :: prefi ! dimension in table of imag refractive indices
89 : real(r8) :: sigmag ! geometric standard deviation of the number distribution for aerosol mode
90 : real(r8) :: dgnum ! geometric dry mean diameter of the number distribution for aerosol mode
91 : real(r8) :: dgnumlo ! lower limit of dgnum
92 : real(r8) :: dgnumhi ! upper limit of dgnum
93 : real(r8) :: rhcrystal ! crystalization relative humidity for mode
94 : real(r8) :: rhdeliques ! deliquescence relative humidity for mode
95 :
96 : endtype physprop_type
97 :
98 : ! This module stores data in an array of physprop_type structures. The way this data
99 : ! is accessed outside the module is via a physprop ID, which is an index into the array.
100 : integer :: numphysprops = 0 ! an incremental total across ALL clim and diag constituents
101 : type (physprop_type), pointer :: physprop(:)
102 :
103 : ! Temporary storage location for filenames in namelist, and construction of dynamic index
104 : ! to properties. The unique filenames specified in the namelist are the identifiers of
105 : ! the properties. Searching the uniquefilenames array provides the index into the physprop
106 : ! array.
107 : character(len=256), allocatable :: uniquefilenames(:)
108 :
109 : ! Number of evenly spaced intervals in rh used in this module and in the aer_rad_props module
110 : ! for calculations of aerosol hygroscopic growth.
111 : integer, parameter, public :: nrh = 1000
112 :
113 : !================================================================================================
114 : contains
115 : !================================================================================================
116 :
117 1536 : subroutine physprop_accum_unique_files(radname, type)
118 :
119 : ! Count number of aerosols in input radname array. Aerosols are identified
120 : ! as strings with a ".nc" suffix.
121 : ! Construct a cumulative list of unique filenames containing physical property data.
122 :
123 : character(len=*), intent(in) :: radname(:)
124 : character(len=1), intent(in) :: type(:)
125 :
126 : integer :: ncnst, i
127 : character(len=*), parameter :: subname = 'physprop_accum_unique_files'
128 : !------------------------------------------------------------------------------------
129 :
130 : ! Initial guess for number of files we need.
131 1536 : if (.not. allocated(uniquefilenames)) allocate(uniquefilenames(50))
132 :
133 1536 : ncnst = ubound(radname, 1)
134 :
135 13824 : do i = 1, ncnst
136 :
137 : ! check if radname is either a bulk aerosol or a mode
138 13824 : if (type(i) == 'A' .or. type(i) == 'M') then
139 :
140 : ! check if this filename has been used by another aerosol. If not
141 : ! then add it to the list of unique names.
142 0 : if (physprop_get_id(radname(i)) < 0) then
143 0 : numphysprops = numphysprops + 1
144 0 : if (numphysprops > size(uniquefilenames)) then
145 0 : call double_capacity(uniquefilenames)
146 : end if
147 0 : uniquefilenames(numphysprops) = trim(radname(i))
148 : endif
149 :
150 : endif
151 : enddo
152 :
153 : contains
154 :
155 : ! Simple routine to re-allocate an array with twice the size, but with
156 : ! the inital values being preserved.
157 0 : subroutine double_capacity(array)
158 : character(len=256), intent(inout), allocatable :: array(:)
159 :
160 0 : character(len=256), allocatable :: tmp(:)
161 : integer :: ierr
162 :
163 0 : allocate(tmp(size(array)*2), stat=ierr)
164 0 : if ( ierr /= 0 ) then
165 0 : call endrun('physprop_accum_unique_files: Allocation error.')
166 : end if
167 :
168 0 : tmp(:size(array)) = array
169 :
170 0 : deallocate(array, stat=ierr)
171 : if ( ierr /= 0 ) then
172 0 : call endrun('physprop_accum_unique_files: Deallocation error.')
173 : end if
174 :
175 0 : call move_alloc(tmp, array)
176 :
177 0 : end subroutine double_capacity
178 :
179 : end subroutine physprop_accum_unique_files
180 :
181 : !================================================================================================
182 :
183 1536 : subroutine physprop_init()
184 :
185 : ! Read properties from the aerosol data files.
186 :
187 : ! ***N.B.*** The calls to physprop_accum_unique_files must be made before calling
188 : ! this init routine. physprop_accum_unique_files is responsible for building
189 : ! the list of files to be read here.
190 :
191 : ! Local variables
192 : integer :: fileindex
193 : type(file_desc_t) :: nc_id ! index to netcdf file
194 : character(len=256) :: locfn ! path to actual file used
195 : character(len=32) :: aername_str ! string read from netCDF file -- may contain trailing
196 : ! nulls which aren't dealt with by trim()
197 :
198 : integer :: ierr ! error codes from mpi
199 :
200 : !------------------------------------------------------------------------------------
201 :
202 3072 : allocate(physprop(numphysprops))
203 :
204 1536 : do fileindex = 1, numphysprops
205 0 : nullify(physprop(fileindex)%sw_hygro_ext)
206 0 : nullify(physprop(fileindex)%sw_hygro_ssa)
207 0 : nullify(physprop(fileindex)%sw_hygro_asm)
208 0 : nullify(physprop(fileindex)%lw_hygro_abs)
209 :
210 0 : nullify(physprop(fileindex)%sw_nonhygro_ext)
211 0 : nullify(physprop(fileindex)%sw_nonhygro_ssa)
212 0 : nullify(physprop(fileindex)%sw_nonhygro_asm)
213 0 : nullify(physprop(fileindex)%sw_nonhygro_scat)
214 0 : nullify(physprop(fileindex)%sw_nonhygro_ascat)
215 0 : nullify(physprop(fileindex)%lw_abs)
216 :
217 0 : nullify(physprop(fileindex)%refindex_aer_sw)
218 0 : nullify(physprop(fileindex)%refindex_aer_lw)
219 :
220 0 : nullify(physprop(fileindex)%r_sw_ext)
221 0 : nullify(physprop(fileindex)%r_sw_scat)
222 0 : nullify(physprop(fileindex)%r_sw_ascat)
223 0 : nullify(physprop(fileindex)%r_lw_abs)
224 0 : nullify(physprop(fileindex)%mu)
225 :
226 0 : nullify(physprop(fileindex)%extpsw)
227 0 : nullify(physprop(fileindex)%abspsw)
228 0 : nullify(physprop(fileindex)%asmpsw)
229 0 : nullify(physprop(fileindex)%absplw)
230 0 : nullify(physprop(fileindex)%refrtabsw)
231 0 : nullify(physprop(fileindex)%refitabsw)
232 0 : nullify(physprop(fileindex)%refrtablw)
233 0 : nullify(physprop(fileindex)%refitablw)
234 :
235 0 : call getfil(uniquefilenames(fileindex), locfn, 0)
236 0 : physprop(fileindex)%sourcefile = locfn
237 :
238 : ! Open the physprop file
239 0 : call cam_pio_openfile(nc_id, locfn, PIO_NOWRITE)
240 :
241 0 : call aerosol_optics_init(physprop(fileindex), nc_id)
242 :
243 : ! Close the physprop file
244 1536 : call pio_closefile(nc_id)
245 :
246 : end do
247 1536 : end subroutine physprop_init
248 :
249 : !================================================================================================
250 :
251 0 : integer function physprop_get_id(filename)
252 :
253 : ! Look for filename in the global list of unique filenames (module data uniquefilenames).
254 : ! If found, return it's index in the list. Otherwise return -1.
255 :
256 : character(len=*), intent(in) :: filename
257 : integer iphysprop
258 :
259 0 : physprop_get_id = -1
260 0 : do iphysprop = 1, numphysprops
261 0 : if(trim(uniquefilenames(iphysprop)) == trim(filename) ) then
262 0 : physprop_get_id = iphysprop
263 0 : return
264 : endif
265 : enddo
266 :
267 : end function physprop_get_id
268 :
269 : !================================================================================================
270 :
271 0 : subroutine physprop_get(id, sourcefile, opticstype, &
272 : sw_hygro_ext, sw_hygro_ssa, sw_hygro_asm, lw_hygro_abs, &
273 : sw_nonhygro_ext, sw_nonhygro_ssa, sw_nonhygro_asm, &
274 : sw_nonhygro_scat, sw_nonhygro_ascat, lw_abs, &
275 : refindex_aer_sw, refindex_aer_lw, &
276 : r_sw_ext, r_sw_scat, r_sw_ascat, r_lw_abs, mu, &
277 : extpsw, abspsw, asmpsw, absplw, refrtabsw, &
278 : refitabsw, refrtablw, refitablw, &
279 0 : aername, density_aer, hygro_aer, dryrad_aer, dispersion_aer, &
280 : num_to_mass_aer, ncoef, prefr, prefi, sigmag, &
281 : dgnum, dgnumlo, dgnumhi, rhcrystal, rhdeliques)
282 :
283 : ! Return requested properties for specified ID.
284 :
285 : ! Arguments
286 : integer, intent(in) :: id
287 : character(len=256), optional, intent(out) :: sourcefile ! Absolute pathname of data file.
288 : character(len=ot_length), optional, intent(out) :: opticstype
289 : real(r8), optional, pointer :: sw_hygro_ext(:,:)
290 : real(r8), optional, pointer :: sw_hygro_ssa(:,:)
291 : real(r8), optional, pointer :: sw_hygro_asm(:,:)
292 : real(r8), optional, pointer :: lw_hygro_abs(:,:)
293 : real(r8), optional, pointer :: sw_nonhygro_ext(:)
294 : real(r8), optional, pointer :: sw_nonhygro_ssa(:)
295 : real(r8), optional, pointer :: sw_nonhygro_asm(:)
296 : real(r8), optional, pointer :: sw_nonhygro_scat(:)
297 : real(r8), optional, pointer :: sw_nonhygro_ascat(:)
298 : real(r8), optional, pointer :: lw_abs(:)
299 : complex(r8), optional, pointer :: refindex_aer_sw(:)
300 : complex(r8), optional, pointer :: refindex_aer_lw(:)
301 : real(r8), optional, pointer :: r_sw_ext(:,:)
302 : real(r8), optional, pointer :: r_sw_scat(:,:)
303 : real(r8), optional, pointer :: r_sw_ascat(:,:)
304 : real(r8), optional, pointer :: r_lw_abs(:,:)
305 : real(r8), optional, pointer :: mu(:)
306 : real(r8), optional, pointer :: extpsw(:,:,:,:)
307 : real(r8), optional, pointer :: abspsw(:,:,:,:)
308 : real(r8), optional, pointer :: asmpsw(:,:,:,:)
309 : real(r8), optional, pointer :: absplw(:,:,:,:)
310 : real(r8), optional, pointer :: refrtabsw(:,:)
311 : real(r8), optional, pointer :: refitabsw(:,:)
312 : real(r8), optional, pointer :: refrtablw(:,:)
313 : real(r8), optional, pointer :: refitablw(:,:)
314 : character(len=20), optional, intent(out) :: aername
315 : real(r8), optional, intent(out) :: density_aer
316 : real(r8), optional, intent(out) :: hygro_aer
317 : real(r8), optional, intent(out) :: dryrad_aer
318 : real(r8), optional, intent(out) :: dispersion_aer
319 : real(r8), optional, intent(out) :: num_to_mass_aer
320 : integer, optional, intent(out) :: ncoef
321 : integer, optional, intent(out) :: prefr
322 : integer, optional, intent(out) :: prefi
323 : real(r8), optional, intent(out) :: sigmag
324 : real(r8), optional, intent(out) :: dgnum
325 : real(r8), optional, intent(out) :: dgnumlo
326 : real(r8), optional, intent(out) :: dgnumhi
327 : real(r8), optional, intent(out) :: rhcrystal
328 : real(r8), optional, intent(out) :: rhdeliques
329 :
330 : ! Local variables
331 : character(len=*), parameter :: subname = 'physprop_get'
332 : !------------------------------------------------------------------------------------
333 :
334 0 : if (id <= 0 .or. id > numphysprops) then
335 0 : write(iulog,*) subname//': illegal ID value: ', id
336 0 : call endrun('physprop_get: ID out of range')
337 : end if
338 :
339 0 : if (present(sourcefile)) sourcefile = physprop(id)%sourcefile
340 0 : if (present(opticstype)) opticstype = physprop(id)%opticsmethod
341 0 : if (present(sw_hygro_ext)) sw_hygro_ext => physprop(id)%sw_hygro_ext
342 0 : if (present(sw_hygro_ssa)) sw_hygro_ssa => physprop(id)%sw_hygro_ssa
343 0 : if (present(sw_hygro_asm)) sw_hygro_asm => physprop(id)%sw_hygro_asm
344 0 : if (present(lw_hygro_abs)) lw_hygro_abs => physprop(id)%lw_hygro_abs
345 0 : if (present(sw_nonhygro_ext)) sw_nonhygro_ext => physprop(id)%sw_nonhygro_ext
346 0 : if (present(sw_nonhygro_ssa)) sw_nonhygro_ssa => physprop(id)%sw_nonhygro_ssa
347 0 : if (present(sw_nonhygro_asm)) sw_nonhygro_asm => physprop(id)%sw_nonhygro_asm
348 0 : if (present(sw_nonhygro_scat)) sw_nonhygro_scat => physprop(id)%sw_nonhygro_scat
349 0 : if (present(sw_nonhygro_ascat)) sw_nonhygro_ascat => physprop(id)%sw_nonhygro_ascat
350 0 : if (present(lw_abs)) lw_abs => physprop(id)%lw_abs
351 :
352 0 : if (present(refindex_aer_sw)) refindex_aer_sw => physprop(id)%refindex_aer_sw
353 0 : if (present(refindex_aer_lw)) refindex_aer_lw => physprop(id)%refindex_aer_lw
354 :
355 0 : if (present(r_sw_ext)) r_sw_ext => physprop(id)%r_sw_ext
356 0 : if (present(r_sw_scat)) r_sw_scat => physprop(id)%r_sw_scat
357 0 : if (present(r_sw_ascat)) r_sw_ascat => physprop(id)%r_sw_ascat
358 0 : if (present(r_lw_abs)) r_lw_abs => physprop(id)%r_lw_abs
359 0 : if (present(mu)) mu => physprop(id)%mu
360 :
361 0 : if (present(extpsw)) extpsw => physprop(id)%extpsw
362 0 : if (present(abspsw)) abspsw => physprop(id)%abspsw
363 0 : if (present(asmpsw)) asmpsw => physprop(id)%asmpsw
364 0 : if (present(absplw)) absplw => physprop(id)%absplw
365 0 : if (present(refrtabsw)) refrtabsw => physprop(id)%refrtabsw
366 0 : if (present(refitabsw)) refitabsw => physprop(id)%refitabsw
367 0 : if (present(refrtablw)) refrtablw => physprop(id)%refrtablw
368 0 : if (present(refitablw)) refitablw => physprop(id)%refitablw
369 :
370 0 : if (present(aername)) aername = physprop(id)%aername
371 0 : if (present(density_aer)) density_aer = physprop(id)%density_aer
372 0 : if (present(hygro_aer)) hygro_aer = physprop(id)%hygro_aer
373 0 : if (present(dryrad_aer)) dryrad_aer = physprop(id)%dryrad_aer
374 0 : if (present(dispersion_aer)) dispersion_aer = physprop(id)%dispersion_aer
375 0 : if (present(num_to_mass_aer)) num_to_mass_aer = physprop(id)%num_to_mass_aer
376 :
377 0 : if (present(ncoef)) ncoef = physprop(id)%ncoef
378 0 : if (present(prefr)) prefr = physprop(id)%prefr
379 0 : if (present(prefi)) prefi = physprop(id)%prefi
380 0 : if (present(sigmag)) sigmag = physprop(id)%sigmag
381 0 : if (present(dgnum)) dgnum = physprop(id)%dgnum
382 0 : if (present(dgnumlo)) dgnumlo = physprop(id)%dgnumlo
383 0 : if (present(dgnumhi)) dgnumhi = physprop(id)%dgnumhi
384 0 : if (present(rhcrystal)) rhcrystal = physprop(id)%rhcrystal
385 0 : if (present(rhdeliques)) rhdeliques = physprop(id)%rhdeliques
386 :
387 0 : end subroutine physprop_get
388 :
389 : !================================================================================================
390 : ! Private methods
391 : !================================================================================================
392 :
393 0 : subroutine aerosol_optics_init(phys_prop, nc_id)
394 :
395 : ! Determine the opticstype, then call the
396 : ! appropriate routine to read the data.
397 :
398 : type(physprop_type), intent(inout) :: phys_prop ! data after interp onto cam rh mesh
399 : type(file_desc_t), intent(inout) :: nc_id ! indentifier for netcdf file
400 :
401 : integer :: opticslength_id, opticslength
402 : type(var_desc_t) :: op_type_id
403 : integer :: ierr ! mpi error codes
404 : character(len=ot_length) :: opticstype_str ! string read from netCDF file -- may contain trailing
405 : ! nulls which aren't dealt with by trim()
406 : !------------------------------------------------------------------------------------
407 :
408 0 : ierr = pio_inq_dimid(nc_id, 'opticsmethod_len', opticslength_id)
409 0 : ierr = pio_inq_dimlen(nc_id, opticslength_id, opticslength)
410 0 : if ( opticslength .gt. ot_length ) then
411 0 : call endrun(" optics type length in "//phys_prop%sourcefile//" excedes maximum length of 32")
412 : endif
413 0 : ierr = pio_inq_varid(nc_id, 'opticsmethod', op_type_id)
414 0 : ierr = pio_get_var(nc_id, op_type_id,phys_prop%opticsmethod )
415 :
416 0 : select case (phys_prop%opticsmethod)
417 : case ('zero')
418 0 : call zero_optics_init(phys_prop, nc_id)
419 :
420 : case ('hygro')
421 0 : call hygro_optics_init(phys_prop, nc_id)
422 :
423 : case ('hygroscopic')
424 0 : call hygroscopic_optics_init(phys_prop, nc_id)
425 :
426 : case ('nonhygro')
427 0 : call nonhygro_optics_init(phys_prop, nc_id)
428 :
429 : case ('insoluble')
430 0 : call insoluble_optics_init(phys_prop, nc_id)
431 :
432 : case ('volcanic_radius','volcanic_radius1','volcanic_radius2','volcanic_radius3')
433 0 : call volcanic_radius_optics_init(phys_prop, nc_id)
434 :
435 : case ('volcanic')
436 0 : call volcanic_optics_init(phys_prop, nc_id)
437 :
438 : case ('modal')
439 0 : call modal_optics_init(phys_prop, nc_id)
440 :
441 : ! other types of optics can be added here
442 :
443 : case default
444 : call endrun('aerosol_optics_init: unsupported optics type '//&
445 0 : trim(phys_prop%opticsmethod)//' in file '//phys_prop%sourcefile)
446 : end select
447 :
448 0 : end subroutine aerosol_optics_init
449 :
450 : !================================================================================================
451 :
452 0 : subroutine hygro_optics_init(phys_prop, nc_id)
453 :
454 : ! Read optics data of type 'hygro' and interpolate it to CAM's rh mesh.
455 :
456 : type (physprop_type), intent(inout) :: phys_prop ! data after interp onto cam rh mesh
457 : type (file_desc_t), intent(inout) :: nc_id ! indentifier for netcdf file
458 :
459 : ! Local variables
460 : integer :: ierr ! error flag
461 :
462 : integer :: rh_idx_id, lw_band_id, sw_band_id
463 : integer :: kbnd, krh
464 : integer :: rh_id, sw_ext_id, sw_ssa_id, sw_asm_id, lw_ext_id
465 : integer :: nbnd, swbands
466 :
467 : ! temp data from hygroscopic file before interpolation onto cam-rh-mesh
468 : integer :: nfilerh ! number of rh values in file
469 0 : real(r8), allocatable, dimension(:) :: frh
470 0 : real(r8), allocatable, dimension(:,:) :: fsw_ext
471 0 : real(r8), allocatable, dimension(:,:) :: fsw_ssa
472 0 : real(r8), allocatable, dimension(:,:) :: fsw_asm
473 :
474 : real(r8) :: rh ! real rh value on cam rh mesh (indexvalue)
475 : !------------------------------------------------------------------------------------
476 :
477 0 : allocate(phys_prop%sw_hygro_ext(nrh,nswbands))
478 0 : allocate(phys_prop%sw_hygro_ssa(nrh,nswbands))
479 0 : allocate(phys_prop%sw_hygro_asm(nrh,nswbands))
480 0 : allocate(phys_prop%lw_abs(nlwbands))
481 :
482 0 : ierr = pio_inq_dimid(nc_id, 'rh_idx', rh_idx_id)
483 :
484 0 : ierr = pio_inq_dimlen(nc_id, rh_idx_id, nfilerh)
485 :
486 0 : ierr = pio_inq_dimid(nc_id, 'lw_band', lw_band_id)
487 :
488 0 : ierr = pio_inq_dimid(nc_id, 'sw_band', sw_band_id)
489 :
490 0 : ierr = pio_inq_dimlen(nc_id, lw_band_id, nbnd)
491 :
492 0 : if (nbnd .ne. nlwbands) call endrun(phys_prop%sourcefile// &
493 0 : ' has the wrong number of lwbands')
494 :
495 0 : ierr = pio_inq_dimlen(nc_id, sw_band_id, swbands)
496 :
497 0 : if(swbands .ne. nswbands) call endrun(phys_prop%sourcefile// &
498 0 : ' has the wrong number of sw bands')
499 :
500 0 : ierr = pio_inq_varid(nc_id, 'rh', rh_id)
501 :
502 0 : ierr = pio_inq_varid(nc_id, 'ext_sw', sw_ext_id)
503 :
504 0 : ierr = pio_inq_varid(nc_id, 'ssa_sw', sw_ssa_id)
505 :
506 0 : ierr = pio_inq_varid(nc_id, 'asm_sw', sw_asm_id)
507 :
508 0 : ierr = pio_inq_varid(nc_id, 'abs_lw', lw_ext_id)
509 :
510 : ! specific optical properties on file's rh mesh
511 0 : allocate(fsw_ext(nfilerh,nswbands))
512 0 : allocate(fsw_asm(nfilerh,nswbands))
513 0 : allocate(fsw_ssa(nfilerh,nswbands))
514 0 : allocate(frh(nfilerh))
515 :
516 0 : ierr = pio_get_var(nc_id, rh_id, frh)
517 :
518 0 : ierr = pio_get_var(nc_id, sw_ext_id, fsw_ext)
519 :
520 0 : ierr = pio_get_var(nc_id, sw_ssa_id, fsw_ssa)
521 :
522 0 : ierr = pio_get_var(nc_id, sw_asm_id, fsw_asm)
523 :
524 0 : ierr = pio_get_var(nc_id, lw_ext_id, phys_prop%lw_abs)
525 :
526 : ! interpolate onto cam's rh mesh
527 0 : do kbnd = 1,nswbands
528 0 : do krh = 1, nrh
529 0 : rh = 1.0_r8 / nrh * (krh - 1)
530 0 : phys_prop%sw_hygro_ext(krh,kbnd) = &
531 0 : exp_interpol( frh, fsw_ext(:,kbnd) / fsw_ext(1,kbnd), rh ) &
532 0 : * fsw_ext(1, kbnd)
533 0 : phys_prop%sw_hygro_ssa(krh,kbnd) = &
534 0 : lin_interpol( frh, fsw_ssa(:,kbnd) / fsw_ssa(1,kbnd), rh ) &
535 0 : * fsw_ssa(1, kbnd)
536 0 : phys_prop%sw_hygro_asm(krh,kbnd) = &
537 0 : lin_interpol( frh, fsw_asm(:,kbnd) / fsw_asm(1,kbnd), rh ) &
538 0 : * fsw_asm(1, kbnd)
539 : enddo
540 : enddo
541 :
542 0 : deallocate (fsw_ext, fsw_asm, fsw_ssa, frh)
543 :
544 : ! read refractive index data if available
545 0 : call refindex_aer_init(phys_prop, nc_id)
546 :
547 : ! read bulk aero props
548 0 : call bulk_props_init(phys_prop, nc_id)
549 :
550 0 : end subroutine hygro_optics_init
551 :
552 : !================================================================================================
553 :
554 0 : subroutine zero_optics_init(phys_prop, nc_id)
555 :
556 : ! Read optics data of type 'nonhygro'
557 :
558 : type (physprop_type), intent(inout) :: phys_prop ! storage for file data
559 : type (file_desc_t), intent(inout) :: nc_id ! indentifier for netcdf file
560 :
561 : ! Local variables
562 : integer :: lw_band_id, sw_band_id
563 : integer :: sw_ext_id, sw_ssa_id, sw_asm_id, lw_ext_id
564 : integer :: swbands, nbnd
565 : integer :: ierr ! error flag
566 : !------------------------------------------------------------------------------------
567 :
568 : ! perhaps this doesn't even need allocated.
569 0 : allocate (phys_prop%sw_nonhygro_ext(nswbands))
570 0 : allocate (phys_prop%sw_nonhygro_ssa(nswbands))
571 0 : allocate (phys_prop%sw_nonhygro_asm(nswbands))
572 0 : allocate (phys_prop%lw_abs(nlwbands))
573 :
574 0 : phys_prop%sw_nonhygro_ext = 0._r8
575 0 : phys_prop%sw_nonhygro_ssa = 0._r8
576 0 : phys_prop%sw_nonhygro_asm = 0._r8
577 0 : phys_prop%lw_abs = 0._r8
578 :
579 0 : end subroutine zero_optics_init
580 :
581 : !================================================================================================
582 :
583 0 : subroutine insoluble_optics_init(phys_prop, nc_id)
584 :
585 : ! Read optics data of type 'nonhygro'
586 :
587 : type (physprop_type), intent(inout) :: phys_prop ! storage for file data
588 : type (file_desc_t), intent(inout) :: nc_id ! indentifier for netcdf file
589 :
590 : ! Local variables
591 : integer :: lw_band_id, sw_band_id
592 : integer :: sw_ext_id, sw_ssa_id, sw_asm_id, lw_ext_id
593 : integer :: swbands, nbnd
594 : integer :: ierr ! error flag
595 : integer :: start(2), count(2)
596 : !------------------------------------------------------------------------------------
597 :
598 0 : allocate (phys_prop%sw_nonhygro_ext(nswbands))
599 0 : allocate (phys_prop%sw_nonhygro_ssa(nswbands))
600 0 : allocate (phys_prop%sw_nonhygro_asm(nswbands))
601 0 : allocate (phys_prop%lw_abs(nlwbands))
602 :
603 0 : ierr = pio_inq_dimid(nc_id, 'lw_band', lw_band_id)
604 :
605 0 : ierr = pio_inq_dimid(nc_id, 'sw_band', sw_band_id)
606 :
607 0 : ierr = pio_inq_dimlen(nc_id, lw_band_id, nbnd)
608 :
609 0 : if (nbnd .ne. nlwbands) call endrun(phys_prop%sourcefile// &
610 0 : ' has the wrong number of lwbands')
611 :
612 0 : ierr = pio_inq_dimlen(nc_id, sw_band_id, swbands)
613 :
614 0 : if (swbands .ne. nswbands) call endrun(phys_prop%sourcefile// &
615 0 : ' has the wrong number of sw bands')
616 :
617 : ! read file data
618 0 : ierr = pio_inq_varid(nc_id, 'ext_sw', sw_ext_id)
619 0 : ierr = pio_inq_varid(nc_id, 'ssa_sw', sw_ssa_id)
620 0 : ierr = pio_inq_varid(nc_id, 'asm_sw', sw_asm_id)
621 0 : ierr = pio_inq_varid(nc_id, 'abs_lw', lw_ext_id)
622 :
623 0 : start = 1
624 0 : count=(/1,swbands/)
625 :
626 0 : ierr = pio_get_var(nc_id, sw_ext_id, start, count, phys_prop%sw_nonhygro_ext)
627 0 : ierr = pio_get_var(nc_id, sw_ssa_id, start, count, phys_prop%sw_nonhygro_ssa)
628 0 : ierr = pio_get_var(nc_id, sw_asm_id, start, count, phys_prop%sw_nonhygro_asm)
629 0 : count = (/1,nbnd/)
630 0 : ierr = pio_get_var(nc_id, lw_ext_id, start, count, phys_prop%lw_abs)
631 :
632 : ! read refractive index data if available
633 0 : call refindex_aer_init(phys_prop, nc_id)
634 :
635 : ! read bulk aero props
636 0 : call bulk_props_init(phys_prop, nc_id)
637 :
638 0 : end subroutine insoluble_optics_init
639 :
640 : !================================================================================================
641 :
642 0 : subroutine volcanic_radius_optics_init(phys_prop, nc_id)
643 :
644 : ! Read optics data of type 'volcanic_radius'
645 :
646 : type (physprop_type), intent(inout) :: phys_prop ! storage for file data
647 : type (file_desc_t), intent(inout) :: nc_id ! indentifier for netcdf file
648 :
649 : ! Local variables
650 : integer :: lw_band_id, sw_band_id, mu_id, mu_did
651 : integer :: sw_ext_id, sw_scat_id, sw_ascat_id, lw_abs_id
652 : integer :: swbands, nbnd, n_mu_samples
653 : integer :: ierr ! error flag
654 : !------------------------------------------------------------------------------------
655 :
656 0 : ierr = pio_inq_dimid(nc_id, 'mu_samples', mu_did)
657 0 : ierr = pio_inq_dimlen(nc_id, mu_did, n_mu_samples)
658 :
659 0 : allocate (phys_prop%r_sw_ext(nswbands,n_mu_samples))
660 0 : allocate (phys_prop%r_sw_scat(nswbands,n_mu_samples))
661 0 : allocate (phys_prop%r_sw_ascat(nswbands,n_mu_samples))
662 0 : allocate (phys_prop%r_lw_abs(nlwbands,n_mu_samples))
663 0 : allocate (phys_prop%mu(n_mu_samples))
664 :
665 0 : ierr = pio_inq_dimid(nc_id, 'lw_band', lw_band_id)
666 :
667 0 : ierr = pio_inq_dimid(nc_id, 'sw_band', sw_band_id)
668 :
669 0 : ierr = pio_inq_dimlen(nc_id, lw_band_id, nbnd)
670 :
671 0 : if (nbnd .ne. nlwbands) call endrun(phys_prop%sourcefile// &
672 0 : ' has the wrong number of lwbands')
673 :
674 0 : ierr = pio_inq_dimlen(nc_id, sw_band_id, swbands)
675 :
676 0 : if (swbands .ne. nswbands) call endrun(phys_prop%sourcefile// &
677 0 : ' has the wrong number of sw bands')
678 :
679 : ! read file data
680 0 : ierr = pio_inq_varid(nc_id, 'bext_sw', sw_ext_id)
681 0 : ierr = pio_inq_varid(nc_id, 'bsca_sw', sw_scat_id)
682 0 : ierr = pio_inq_varid(nc_id, 'basc_sw', sw_ascat_id)
683 0 : ierr = pio_inq_varid(nc_id, 'babs_lw', lw_abs_id)
684 0 : ierr = pio_inq_varid(nc_id, 'mu_samples', mu_id)
685 :
686 0 : ierr = pio_get_var(nc_id, sw_ext_id, phys_prop%r_sw_ext)
687 0 : ierr = pio_get_var(nc_id, sw_scat_id, phys_prop%r_sw_scat)
688 0 : ierr = pio_get_var(nc_id, sw_ascat_id, phys_prop%r_sw_ascat)
689 0 : ierr = pio_get_var(nc_id, lw_abs_id, phys_prop%r_lw_abs)
690 0 : ierr = pio_get_var(nc_id, mu_id, phys_prop%mu)
691 :
692 : ! read bulk aero props
693 0 : call bulk_props_init(phys_prop, nc_id)
694 :
695 0 : end subroutine volcanic_radius_optics_init
696 :
697 : !================================================================================================
698 :
699 0 : subroutine volcanic_optics_init(phys_prop, nc_id)
700 :
701 : ! Read optics data of type 'volcanic'
702 :
703 : type (physprop_type), intent(inout) :: phys_prop ! storage for file data
704 : type (file_desc_t) , intent(inout) :: nc_id ! indentifier for netcdf file
705 :
706 : ! Local variables
707 : integer :: lw_band_id, sw_band_id
708 : integer :: sw_ext_id, sw_scat_id, sw_ascat_id, lw_abs_id
709 : integer :: swbands, nbnd
710 : integer :: ierr ! error flag
711 : !------------------------------------------------------------------------------------
712 :
713 0 : allocate (phys_prop%sw_nonhygro_ext(nswbands))
714 0 : allocate (phys_prop%sw_nonhygro_scat(nswbands))
715 0 : allocate (phys_prop%sw_nonhygro_ascat(nswbands))
716 0 : allocate (phys_prop%lw_abs(nlwbands))
717 :
718 0 : ierr = pio_inq_dimid(nc_id, 'lw_band', lw_band_id)
719 0 : ierr = pio_inq_dimid(nc_id, 'sw_band', sw_band_id)
720 :
721 0 : ierr = pio_inq_dimlen(nc_id, lw_band_id, nbnd)
722 :
723 0 : if (nbnd .ne. nlwbands) call endrun(phys_prop%sourcefile// &
724 0 : ' has the wrong number of lwbands')
725 :
726 0 : ierr = pio_inq_dimlen(nc_id, sw_band_id, swbands)
727 0 : if(masterproc) write(iulog,*) 'swbands',swbands
728 :
729 0 : if (swbands .ne. nswbands) call endrun(phys_prop%sourcefile// &
730 0 : ' has the wrong number of sw bands')
731 :
732 : ! read file data
733 0 : ierr = pio_inq_varid(nc_id, 'bext_sw', sw_ext_id)
734 0 : ierr = pio_inq_varid(nc_id, 'bsca_sw', sw_scat_id)
735 0 : ierr = pio_inq_varid(nc_id, 'basc_sw', sw_ascat_id)
736 0 : ierr = pio_inq_varid(nc_id, 'babs_lw', lw_abs_id)
737 :
738 0 : ierr = pio_get_var(nc_id, sw_ext_id, phys_prop%sw_nonhygro_ext)
739 0 : ierr = pio_get_var(nc_id, sw_scat_id, phys_prop%sw_nonhygro_scat)
740 0 : ierr = pio_get_var(nc_id, sw_ascat_id, phys_prop%sw_nonhygro_ascat)
741 0 : ierr = pio_get_var(nc_id, lw_abs_id, phys_prop%lw_abs)
742 :
743 : ! read bulk aero props
744 0 : call bulk_props_init(phys_prop, nc_id)
745 :
746 0 : end subroutine volcanic_optics_init
747 :
748 : !================================================================================================
749 :
750 0 : subroutine hygroscopic_optics_init(phys_prop, nc_id)
751 :
752 : ! Read optics data of type 'hygroscopic' and interpolate it to CAM's rh mesh.
753 :
754 : type (physprop_type), intent(inout) :: phys_prop ! data after interp onto cam rh mesh
755 : type (file_desc_T), intent(inout) :: nc_id ! indentifier for netcdf file
756 :
757 : ! Local variables
758 : integer :: ierr ! error flag
759 :
760 : integer :: rh_idx_id, lw_band_id, sw_band_id
761 : integer :: kbnd, krh
762 : integer :: rh_id, sw_ext_id, sw_ssa_id, sw_asm_id, lw_ext_id
763 : integer :: nbnd, swbands
764 :
765 : ! temp data from hygroscopic file before interpolation onto cam-rh-mesh
766 : integer :: nfilerh ! number of rh values in file
767 0 : real(r8), allocatable, dimension(:) :: frh
768 0 : real(r8), allocatable, dimension(:,:) :: fsw_ext
769 0 : real(r8), allocatable, dimension(:,:) :: fsw_ssa
770 0 : real(r8), allocatable, dimension(:,:) :: fsw_asm
771 0 : real(r8), allocatable, dimension(:,:) :: flw_abs
772 :
773 : real(r8) :: rh ! real rh value on cam rh mesh (indexvalue)
774 : character(len=*), parameter :: sub = 'hygroscopic_optics_init'
775 : !------------------------------------------------------------------------------------
776 :
777 0 : allocate(phys_prop%sw_hygro_ext(nrh,nswbands))
778 0 : allocate(phys_prop%sw_hygro_ssa(nrh,nswbands))
779 0 : allocate(phys_prop%sw_hygro_asm(nrh,nswbands))
780 0 : allocate(phys_prop%lw_hygro_abs(nrh,nlwbands))
781 :
782 0 : ierr = pio_inq_dimid(nc_id, 'rh_idx', rh_idx_id)
783 0 : ierr = pio_inq_dimlen(nc_id, rh_idx_id, nfilerh)
784 :
785 0 : ierr = pio_inq_dimid(nc_id, 'lw_band', lw_band_id)
786 0 : ierr = pio_inq_dimlen(nc_id, lw_band_id, nbnd)
787 0 : if (nbnd .ne. nlwbands) call endrun(phys_prop%sourcefile// &
788 0 : ' has the wrong number of lwbands')
789 :
790 0 : ierr = pio_inq_dimid(nc_id, 'sw_band', sw_band_id)
791 0 : ierr = pio_inq_dimlen(nc_id, sw_band_id, swbands)
792 0 : if(swbands .ne. nswbands) call endrun(phys_prop%sourcefile// &
793 0 : ' has the wrong number of sw bands')
794 :
795 0 : ierr = pio_inq_varid(nc_id, 'rh', rh_id)
796 0 : ierr = pio_inq_varid(nc_id, 'ext_sw', sw_ext_id)
797 0 : ierr = pio_inq_varid(nc_id, 'ssa_sw', sw_ssa_id)
798 0 : ierr = pio_inq_varid(nc_id, 'asm_sw', sw_asm_id)
799 0 : ierr = pio_inq_varid(nc_id, 'abs_lw', lw_ext_id)
800 :
801 : ! specific optical properties on file's rh mesh
802 0 : allocate(fsw_ext(nfilerh,nswbands))
803 0 : allocate(fsw_asm(nfilerh,nswbands))
804 0 : allocate(fsw_ssa(nfilerh,nswbands))
805 0 : allocate(flw_abs(nfilerh,nlwbands))
806 0 : allocate(frh(nfilerh))
807 :
808 0 : ierr = pio_get_var(nc_id, rh_id, frh)
809 0 : ierr = pio_get_var(nc_id, sw_ext_id, fsw_ext)
810 0 : ierr = pio_get_var(nc_id, sw_ssa_id, fsw_ssa)
811 0 : ierr = pio_get_var(nc_id, sw_asm_id, fsw_asm)
812 0 : ierr = pio_get_var(nc_id, lw_ext_id, flw_abs)
813 :
814 : ! interpolate onto cam's rh mesh
815 0 : do kbnd = 1,nswbands
816 0 : do krh = 1, nrh
817 0 : rh = 1.0_r8 / nrh * (krh - 1)
818 0 : phys_prop%sw_hygro_ext(krh,kbnd) = &
819 0 : exp_interpol( frh, fsw_ext(:,kbnd) / fsw_ext(1,kbnd), rh ) &
820 0 : * fsw_ext(1, kbnd)
821 0 : phys_prop%sw_hygro_ssa(krh,kbnd) = &
822 0 : lin_interpol( frh, fsw_ssa(:,kbnd) / fsw_ssa(1,kbnd), rh ) &
823 0 : * fsw_ssa(1, kbnd)
824 0 : phys_prop%sw_hygro_asm(krh,kbnd) = &
825 0 : lin_interpol( frh, fsw_asm(:,kbnd) / fsw_asm(1,kbnd), rh ) &
826 0 : * fsw_asm(1, kbnd)
827 : enddo
828 : enddo
829 0 : do kbnd = 1,nlwbands
830 0 : do krh = 1, nrh
831 0 : rh = 1.0_r8 / nrh * (krh - 1)
832 0 : phys_prop%lw_hygro_abs(krh,kbnd) = &
833 0 : exp_interpol( frh, flw_abs(:,kbnd) / flw_abs(1,kbnd), rh ) &
834 0 : * flw_abs(1, kbnd)
835 : enddo
836 : enddo
837 :
838 0 : deallocate (fsw_ext, fsw_asm, fsw_ssa, flw_abs, frh)
839 :
840 : ! read refractive index data if available
841 0 : call refindex_aer_init(phys_prop, nc_id)
842 :
843 : ! read bulk aero props
844 0 : call bulk_props_init(phys_prop, nc_id)
845 :
846 0 : end subroutine hygroscopic_optics_init
847 :
848 : !================================================================================================
849 :
850 0 : subroutine nonhygro_optics_init(phys_prop, nc_id)
851 :
852 : ! Read optics data of type 'nonhygro'
853 :
854 : type (physprop_type), intent(inout) :: phys_prop ! storage for file data
855 : type (file_desc_t) , intent(inout) :: nc_id ! indentifier for netcdf file
856 :
857 : ! Local variables
858 : integer :: lw_band_id, sw_band_id
859 : integer :: sw_ext_id, sw_ssa_id, sw_asm_id, lw_ext_id
860 : integer :: swbands, nbnd
861 : integer :: ierr ! error flag
862 : !------------------------------------------------------------------------------------
863 :
864 0 : allocate (phys_prop%sw_nonhygro_ext(nswbands))
865 0 : allocate (phys_prop%sw_nonhygro_ssa(nswbands))
866 0 : allocate (phys_prop%sw_nonhygro_asm(nswbands))
867 0 : allocate (phys_prop%lw_abs(nlwbands))
868 :
869 0 : ierr = pio_inq_dimid(nc_id, 'lw_band', lw_band_id)
870 0 : ierr = pio_inq_dimid(nc_id, 'sw_band', sw_band_id)
871 :
872 0 : ierr = pio_inq_dimlen(nc_id, lw_band_id, nbnd)
873 :
874 0 : if (nbnd .ne. nlwbands) call endrun(phys_prop%sourcefile// &
875 0 : ' has the wrong number of lwbands')
876 :
877 0 : ierr = pio_inq_dimlen(nc_id, sw_band_id, swbands)
878 :
879 0 : if (swbands .ne. nswbands) call endrun(phys_prop%sourcefile// &
880 0 : ' has the wrong number of sw bands')
881 :
882 : ! read file data
883 0 : ierr = pio_inq_varid(nc_id, 'ext_sw', sw_ext_id)
884 0 : ierr = pio_inq_varid(nc_id, 'ssa_sw', sw_ssa_id)
885 0 : ierr = pio_inq_varid(nc_id, 'asm_sw', sw_asm_id)
886 0 : ierr = pio_inq_varid(nc_id, 'abs_lw', lw_ext_id)
887 :
888 0 : ierr = pio_get_var(nc_id, sw_ext_id, phys_prop%sw_nonhygro_ext)
889 0 : ierr = pio_get_var(nc_id, sw_ssa_id, phys_prop%sw_nonhygro_ssa)
890 0 : ierr = pio_get_var(nc_id, sw_asm_id, phys_prop%sw_nonhygro_asm)
891 0 : ierr = pio_get_var(nc_id, lw_ext_id, phys_prop%lw_abs)
892 :
893 : ! read refractive index data if available
894 0 : call refindex_aer_init(phys_prop, nc_id)
895 :
896 : ! read bulk aero props
897 0 : call bulk_props_init(phys_prop, nc_id)
898 :
899 0 : end subroutine nonhygro_optics_init
900 :
901 : !================================================================================================
902 :
903 0 : subroutine refindex_aer_init(phys_prop, nc_id)
904 :
905 : ! Read refractive indices of aerosol
906 :
907 : type (physprop_type), intent(inout) :: phys_prop ! storage for file data
908 : type (file_desc_T), intent(inout) :: nc_id ! indentifier for netcdf file
909 :
910 : ! Local variables
911 : integer :: i
912 : integer :: istat1, istat2, istat3 ! status flags
913 : integer :: vid_real, vid_im ! variable ids
914 0 : real(r8), pointer :: ref_real(:), ref_im(:) ! tmp storage for components of complex index
915 : character(len=*), parameter :: subname = 'refindex_aer_init'
916 : !------------------------------------------------------------------------------------
917 :
918 : ! assume that the dimensions lw_band and sw_band have already been checked
919 : ! by the calling subroutine
920 :
921 : ! Check that the variables are present before allocating storage and reading.
922 : ! Since we're setting complex data values, both the real and imaginary parts must
923 : ! be present or neither will be read.
924 :
925 : ! set PIO to return control to the caller when variable not found
926 0 : call pio_seterrorhandling(nc_id, PIO_BCAST_ERROR)
927 :
928 0 : istat1 = pio_inq_varid(nc_id, 'refindex_real_aer_sw', vid_real)
929 0 : istat2 = pio_inq_varid(nc_id, 'refindex_im_aer_sw', vid_im)
930 :
931 0 : if (istat1 == PIO_NOERR .and. istat2 == PIO_NOERR) then
932 :
933 0 : allocate(ref_real(nswbands), ref_im(nswbands))
934 :
935 0 : istat3 = pio_get_var(nc_id, vid_real, ref_real)
936 0 : if (istat3 /= PIO_NOERR) then
937 0 : call endrun(subname//': ERROR reading refindex_real_aer_sw')
938 : end if
939 :
940 0 : istat3 = pio_get_var(nc_id, vid_im, ref_im)
941 0 : if (istat3 /= PIO_NOERR) then
942 0 : call endrun(subname//': ERROR reading refindex_im_aer_sw')
943 : end if
944 :
945 : ! successfully read refindex data -- set complex values in physprop object
946 0 : allocate(phys_prop%refindex_aer_sw(nswbands))
947 0 : do i = 1, nswbands
948 0 : phys_prop%refindex_aer_sw(i) = cmplx(ref_real(i), abs(ref_im(i)),&
949 0 : kind=r8)
950 : end do
951 :
952 0 : deallocate(ref_real, ref_im)
953 :
954 : end if
955 :
956 0 : istat1 = pio_inq_varid(nc_id, 'refindex_real_aer_lw', vid_real)
957 0 : istat2 = pio_inq_varid(nc_id, 'refindex_im_aer_lw', vid_im)
958 :
959 0 : if (istat1 == PIO_NOERR .and. istat2 == PIO_NOERR) then
960 :
961 0 : allocate(ref_real(nlwbands), ref_im(nlwbands))
962 :
963 0 : istat3 = pio_get_var(nc_id, vid_real, ref_real)
964 0 : if (istat3 /= PIO_NOERR) then
965 0 : call endrun(subname//': ERROR reading refindex_real_aer_lw')
966 : end if
967 :
968 0 : istat3 = pio_get_var(nc_id, vid_im, ref_im)
969 0 : if (istat3 /= PIO_NOERR) then
970 0 : call endrun(subname//': ERROR reading refindex_im_aer_lw')
971 : end if
972 :
973 : ! successfully read refindex data -- set complex value in physprop object
974 0 : allocate(phys_prop%refindex_aer_lw(nlwbands))
975 0 : do i = 1, nlwbands
976 0 : phys_prop%refindex_aer_lw(i) = cmplx(ref_real(i), abs(ref_im(i)),&
977 0 : kind=r8)
978 : end do
979 :
980 0 : deallocate(ref_real, ref_im)
981 :
982 : end if
983 :
984 : ! reset PIO to handle errors internally
985 0 : call pio_seterrorhandling(nc_id, PIO_INTERNAL_ERROR)
986 :
987 0 : end subroutine refindex_aer_init
988 :
989 : !================================================================================================
990 :
991 0 : subroutine modal_optics_init(props, ncid)
992 :
993 : ! Read optics data for modal aerosols
994 :
995 : type (physprop_type), intent(inout) :: props ! storage for file data
996 : type (file_desc_T), intent(inout) :: ncid ! indentifier for netcdf file
997 :
998 : ! Local variables
999 : integer :: ierr
1000 : integer :: did
1001 : integer :: ival
1002 : type(var_desc_t) :: vid
1003 0 : real(r8), pointer :: rval(:,:,:,:,:) ! temp array used to eliminate a singleton dimension
1004 :
1005 : character(len=*), parameter :: subname = 'modal_optics_init'
1006 : !------------------------------------------------------------------------------------
1007 :
1008 : ! Check dimensions for number of lw and sw bands
1009 :
1010 0 : ierr = pio_inq_dimid(ncid, 'lw_band', did)
1011 0 : ierr = pio_inq_dimlen(ncid, did, ival)
1012 0 : if (ival .ne. nlwbands) call endrun(subname//':'//props%sourcefile// &
1013 0 : ' has the wrong number of lw bands')
1014 :
1015 0 : ierr = pio_inq_dimid(ncid, 'sw_band', did)
1016 0 : ierr = pio_inq_dimlen(ncid, did, ival)
1017 0 : if (ival .ne. nswbands) call endrun(subname//':'//props%sourcefile// &
1018 0 : ' has the wrong number of sw bands')
1019 :
1020 : ! Get other dimensions
1021 0 : ierr = pio_inq_dimid(ncid, 'coef_number', did)
1022 0 : ierr = pio_inq_dimlen(ncid, did, props%ncoef)
1023 :
1024 0 : ierr = pio_inq_dimid(ncid, 'refindex_real', did)
1025 0 : ierr = pio_inq_dimlen(ncid, did, props%prefr)
1026 :
1027 0 : ierr = pio_inq_dimid(ncid, 'refindex_im', did)
1028 0 : ierr = pio_inq_dimlen(ncid, did, props%prefi)
1029 :
1030 : ! Allocate arrays
1031 : allocate( &
1032 : props%extpsw(props%ncoef,props%prefr,props%prefi,nswbands), &
1033 : props%abspsw(props%ncoef,props%prefr,props%prefi,nswbands), &
1034 : props%asmpsw(props%ncoef,props%prefr,props%prefi,nswbands), &
1035 : props%absplw(props%ncoef,props%prefr,props%prefi,nlwbands), &
1036 : props%refrtabsw(props%prefr,nswbands), &
1037 : props%refitabsw(props%prefi,nswbands), &
1038 : props%refrtablw(props%prefr,nlwbands), &
1039 0 : props%refitablw(props%prefi,nlwbands) )
1040 :
1041 :
1042 : ! allocate temp to remove the mode dimension from the sw variables
1043 0 : allocate(rval(props%ncoef,props%prefr,props%prefi,1,nswbands))
1044 :
1045 0 : ierr = pio_inq_varid(ncid, 'extpsw', vid)
1046 0 : ierr = pio_get_var(ncid, vid, rval)
1047 0 : props%extpsw = rval(:,:,:,1,:)
1048 :
1049 0 : ierr = pio_inq_varid(ncid, 'abspsw', vid)
1050 0 : ierr = pio_get_var(ncid, vid, rval)
1051 0 : props%abspsw = rval(:,:,:,1,:)
1052 :
1053 0 : ierr = pio_inq_varid(ncid, 'asmpsw', vid)
1054 0 : ierr = pio_get_var(ncid, vid, rval)
1055 0 : props%asmpsw = rval(:,:,:,1,:)
1056 :
1057 0 : deallocate(rval)
1058 :
1059 : ! allocate temp to remove the mode dimension from the lw variables
1060 0 : allocate(rval(props%ncoef,props%prefr,props%prefi,1,nlwbands))
1061 :
1062 0 : ierr = pio_inq_varid(ncid, 'absplw', vid)
1063 0 : ierr = pio_get_var(ncid, vid, rval)
1064 0 : props%absplw = rval(:,:,:,1,:)
1065 :
1066 0 : deallocate(rval)
1067 :
1068 0 : ierr = pio_inq_varid(ncid, 'refindex_real_sw', vid)
1069 0 : ierr = pio_get_var(ncid, vid, props%refrtabsw)
1070 :
1071 0 : ierr = pio_inq_varid(ncid, 'refindex_im_sw', vid)
1072 0 : ierr = pio_get_var(ncid, vid, props%refitabsw)
1073 :
1074 0 : ierr = pio_inq_varid(ncid, 'refindex_real_lw', vid)
1075 0 : ierr = pio_get_var(ncid, vid, props%refrtablw)
1076 :
1077 0 : ierr = pio_inq_varid(ncid, 'refindex_im_lw', vid)
1078 0 : ierr = pio_get_var(ncid, vid, props%refitablw)
1079 :
1080 0 : ierr = pio_inq_varid(ncid, 'sigmag', vid)
1081 0 : ierr = pio_get_var(ncid, vid, props%sigmag)
1082 :
1083 0 : ierr = pio_inq_varid(ncid, 'dgnum', vid)
1084 0 : ierr = pio_get_var(ncid, vid, props%dgnum)
1085 :
1086 0 : ierr = pio_inq_varid(ncid, 'dgnumlo', vid)
1087 0 : ierr = pio_get_var(ncid, vid, props%dgnumlo)
1088 :
1089 0 : ierr = pio_inq_varid(ncid, 'dgnumhi', vid)
1090 0 : ierr = pio_get_var(ncid, vid, props%dgnumhi)
1091 :
1092 0 : ierr = pio_inq_varid(ncid, 'rhcrystal', vid)
1093 0 : ierr = pio_get_var(ncid, vid, props%rhcrystal)
1094 :
1095 0 : ierr = pio_inq_varid(ncid, 'rhdeliques', vid)
1096 0 : ierr = pio_get_var(ncid, vid, props%rhdeliques)
1097 :
1098 0 : end subroutine modal_optics_init
1099 :
1100 : !================================================================================================
1101 :
1102 0 : subroutine bulk_props_init(physprop, nc_id)
1103 :
1104 : ! Read props for bulk aerosols
1105 :
1106 : type (physprop_type), intent(inout) :: physprop ! storage for file data
1107 : type (file_desc_T), intent(inout) :: nc_id ! indentifier for netcdf file
1108 :
1109 : ! Local variables
1110 : integer :: ierr
1111 :
1112 : type(var_desc_T) :: vid
1113 :
1114 : logical :: debug = .false.
1115 :
1116 : character(len=*), parameter :: subname = 'bulk_props_init'
1117 : !------------------------------------------------------------------------------------
1118 :
1119 : ! read microphys
1120 0 : ierr = pio_inq_varid(nc_id, 'name', vid)
1121 0 : ierr = pio_get_var(nc_id, vid, physprop%aername)
1122 :
1123 : ! use GLC function to remove trailing nulls and blanks.
1124 : ! physprop%aername = aername_str(:GLC(aername_str))
1125 :
1126 0 : ierr = pio_inq_varid(nc_id, 'density', vid)
1127 0 : ierr = pio_get_var(nc_id, vid, physprop%density_aer)
1128 :
1129 0 : ierr = pio_inq_varid(nc_id, 'sigma_logr', vid)
1130 0 : ierr = pio_get_var(nc_id, vid, physprop%dispersion_aer)
1131 :
1132 0 : ierr = pio_inq_varid(nc_id, 'dryrad', vid)
1133 0 : ierr = pio_get_var(nc_id, vid, physprop%dryrad_aer)
1134 :
1135 0 : ierr = pio_inq_varid(nc_id, 'hygroscopicity', vid)
1136 0 : ierr = pio_get_var(nc_id, vid, physprop%hygro_aer)
1137 :
1138 0 : ierr = pio_inq_varid(nc_id, 'num_to_mass_ratio', vid)
1139 0 : ierr = pio_get_var(nc_id, vid, physprop%num_to_mass_aer)
1140 :
1141 : ! Output select data to log file
1142 0 : if (debug .and. masterproc .and. idx_sw_diag > 0) then
1143 0 : if (trim(physprop%aername) == 'SULFATE') then
1144 0 : write(iulog, '(2x, a)') '_______ hygroscopic growth in visible band _______'
1145 0 : call aer_optics_log_rh('SO4', physprop%sw_hygro_ext(:,idx_sw_diag), &
1146 0 : physprop%sw_hygro_ssa(:,idx_sw_diag), physprop%sw_hygro_asm(:,idx_sw_diag))
1147 : end if
1148 0 : write(iulog, *) subname//': finished for ', trim(physprop%aername)
1149 : end if
1150 :
1151 0 : end subroutine bulk_props_init
1152 :
1153 : !================================================================================================
1154 :
1155 0 : function exp_interpol(x, f, y) result(g)
1156 : ! Purpose:
1157 : ! interpolates f(x) to point y
1158 : ! assuming f(x) = f(x0) exp a(x - x0)
1159 : ! where a = ( ln f(x1) - ln f(x0) ) / (x1 - x0)
1160 : ! x0 <= x <= x1
1161 : ! assumes x is monotonically increasing
1162 : ! Author: D. Fillmore
1163 :
1164 : implicit none
1165 :
1166 : real(r8), intent(in), dimension(:) :: x ! grid points
1167 : real(r8), intent(in), dimension(:) :: f ! grid function values
1168 : real(r8), intent(in) :: y ! interpolation point
1169 : real(r8) :: g ! interpolated function value
1170 :
1171 : integer :: k ! interpolation point index
1172 : integer :: n ! length of x
1173 : real(r8) :: a
1174 :
1175 0 : n = size(x)
1176 :
1177 : ! find k such that x(k) < y =< x(k+1)
1178 : ! set k = 1 if y <= x(1) and k = n-1 if y > x(n)
1179 :
1180 0 : if (y <= x(1)) then
1181 : k = 1
1182 0 : else if (y >= x(n)) then
1183 0 : k = n - 1
1184 : else
1185 : k = 1
1186 0 : do while (y > x(k+1) .and. k < n)
1187 0 : k = k + 1
1188 : end do
1189 : end if
1190 :
1191 : ! interpolate
1192 0 : a = ( log( f(k+1) / f(k) ) ) / ( x(k+1) - x(k) )
1193 0 : g = f(k) * exp( a * (y - x(k)) )
1194 : return
1195 : end function exp_interpol
1196 :
1197 : !================================================================================================
1198 :
1199 0 : function lin_interpol(x, f, y) result(g)
1200 : ! Purpose:
1201 : ! interpolates f(x) to point y
1202 : ! assuming f(x) = f(x0) + a * (x - x0)
1203 : ! where a = ( f(x1) - f(x0) ) / (x1 - x0)
1204 : ! x0 <= x <= x1
1205 : ! assumes x is monotonically increasing
1206 : ! Author: D. Fillmore
1207 :
1208 : implicit none
1209 :
1210 : real(r8), intent(in), dimension(:) :: x ! grid points
1211 : real(r8), intent(in), dimension(:) :: f ! grid function values
1212 : real(r8), intent(in) :: y ! interpolation point
1213 : real(r8) :: g ! interpolated function value
1214 :
1215 : integer :: k ! interpolation point index
1216 : integer :: n ! length of x
1217 : real(r8) :: a
1218 :
1219 0 : n = size(x)
1220 :
1221 : ! find k such that x(k) < y =< x(k+1)
1222 : ! set k = 1 if y <= x(1) and k = n-1 if y > x(n)
1223 :
1224 0 : if (y <= x(1)) then
1225 : k = 1
1226 0 : else if (y >= x(n)) then
1227 0 : k = n - 1
1228 : else
1229 : k = 1
1230 0 : do while (y > x(k+1) .and. k < n)
1231 0 : k = k + 1
1232 : end do
1233 : end if
1234 :
1235 : ! interpolate
1236 0 : a = ( f(k+1) - f(k) ) / ( x(k+1) - x(k) )
1237 0 : g = f(k) + a * (y - x(k))
1238 : return
1239 : end function lin_interpol
1240 :
1241 : !================================================================================================
1242 :
1243 : subroutine aer_optics_log(name, ext, ssa, asm)
1244 :
1245 : ! Purpose:
1246 : ! write aerosol optical constants to log file
1247 :
1248 : ! Author: D. Fillmore
1249 :
1250 : character(len=*), intent(in) :: name
1251 : real(r8), intent(in) :: ext(:)
1252 : real(r8), intent(in) :: ssa(:)
1253 : real(r8), intent(in) :: asm(:)
1254 :
1255 : integer :: kbnd, nbnd
1256 : !------------------------------------------------------------------------------------
1257 :
1258 : nbnd = ubound(ext, 1)
1259 :
1260 : write(iulog, '(2x, a)') name
1261 : write(iulog, '(2x, a, 4x, a, 4x, a, 4x, a)') 'SW band', 'ext (m^2 kg^-1)', ' ssa', ' asm'
1262 : do kbnd = 1, nbnd
1263 : write(iulog, '(2x, i7, 4x, f13.2, 4x, f4.2, 4x, f4.2)') kbnd, ext(kbnd), ssa(kbnd), asm(kbnd)
1264 : end do
1265 :
1266 : end subroutine aer_optics_log
1267 :
1268 : !================================================================================================
1269 :
1270 :
1271 0 : subroutine aer_optics_log_rh(name, ext, ssa, asm)
1272 :
1273 : ! Purpose:
1274 : ! write out aerosol optical properties
1275 : ! for a set of test rh values
1276 : ! to test hygroscopic growth interpolation
1277 :
1278 : ! Author: D. Fillmore
1279 :
1280 : character(len=*), intent(in) :: name
1281 : real(r8), intent(in) :: ext(nrh)
1282 : real(r8), intent(in) :: ssa(nrh)
1283 : real(r8), intent(in) :: asm(nrh)
1284 :
1285 : integer :: krh_test
1286 : integer, parameter :: nrh_test = 36
1287 : integer :: krh
1288 : real(r8) :: rh
1289 : real(r8) :: rh_test(nrh_test)
1290 : real(r8) :: exti
1291 : real(r8) :: ssai
1292 : real(r8) :: asmi
1293 : real(r8) :: wrh
1294 : !------------------------------------------------------------------------------------
1295 :
1296 0 : do krh_test = 1, nrh_test
1297 0 : rh_test(krh_test) = sqrt(sqrt(sqrt(sqrt(((krh_test - 1.0_r8) / (nrh_test - 1))))))
1298 : enddo
1299 0 : write(iulog, '(2x, a)') name
1300 0 : write(iulog, '(2x, a, 4x, a, 4x, a, 4x, a)') ' rh', 'ext (m^2 kg^-1)', ' ssa', ' asm'
1301 :
1302 : ! loop through test rh values
1303 0 : do krh_test = 1, nrh_test
1304 : ! find corresponding rh index
1305 0 : rh = rh_test(krh_test)
1306 0 : krh = min(floor( (rh) * nrh ) + 1, nrh - 1)
1307 0 : wrh = (rh) *nrh - krh
1308 0 : exti = ext(krh + 1) * (wrh + 1) - ext(krh) * wrh
1309 0 : ssai = ssa(krh + 1) * (wrh + 1) - ssa(krh) * wrh
1310 0 : asmi = asm(krh + 1) * (wrh + 1) - asm(krh) * wrh
1311 0 : write(iulog, '(2x, f5.3, 4x, f13.3, 4x, f5.3, 4x, f5.3)') rh_test(krh_test), exti, ssai, asmi
1312 : end do
1313 :
1314 0 : end subroutine aer_optics_log_rh
1315 :
1316 :
1317 : !================================================================================================
1318 :
1319 0 : end module phys_prop
|
