1 | ! radiation_general_cloud_optics_data.F90 - Type to store generalized cloud optical properties |
---|
2 | ! |
---|
3 | ! (C) Copyright 2019- ECMWF. |
---|
4 | ! |
---|
5 | ! This software is licensed under the terms of the Apache Licence Version 2.0 |
---|
6 | ! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. |
---|
7 | ! |
---|
8 | ! In applying this licence, ECMWF does not waive the privileges and immunities |
---|
9 | ! granted to it by virtue of its status as an intergovernmental organisation |
---|
10 | ! nor does it submit to any jurisdiction. |
---|
11 | ! |
---|
12 | ! Author: Robin Hogan |
---|
13 | ! Email: r.j.hogan@ecmwf.int |
---|
14 | ! License: see the COPYING file for details |
---|
15 | ! |
---|
16 | |
---|
17 | module radiation_general_cloud_optics_data |
---|
18 | |
---|
19 | use parkind1, only : jprb |
---|
20 | |
---|
21 | implicit none |
---|
22 | |
---|
23 | public |
---|
24 | |
---|
25 | !--------------------------------------------------------------------- |
---|
26 | ! This type holds the configuration information to compute optical |
---|
27 | ! properties for a particular type of cloud or hydrometeor in one of |
---|
28 | ! the shortwave or longwave |
---|
29 | type general_cloud_optics_type |
---|
30 | ! Band-specific (or g-point-specific) values as a look-up table |
---|
31 | ! versus effective radius dimensioned (nband,n_effective_radius) |
---|
32 | |
---|
33 | ! Extinction coefficient per unit mass (m2 kg-1) |
---|
34 | real(jprb), allocatable, dimension(:,:) :: & |
---|
35 | & mass_ext |
---|
36 | |
---|
37 | ! Single-scattering albedo and asymmetry factor (dimensionless) |
---|
38 | real(jprb), allocatable, dimension(:,:) :: & |
---|
39 | & ssa, asymmetry |
---|
40 | |
---|
41 | ! Number of effective radius coefficients, start value and |
---|
42 | ! interval in look-up table |
---|
43 | integer :: n_effective_radius = 0 |
---|
44 | real(jprb) :: effective_radius_0, d_effective_radius |
---|
45 | |
---|
46 | ! Name of cloud/precip type and scattering model |
---|
47 | ! (e.g. "mie_droplet", "fu-muskatel_ice"). These are used to |
---|
48 | ! generate the name of the data file from which the coefficients |
---|
49 | ! are read. |
---|
50 | character(len=511) :: type_name |
---|
51 | |
---|
52 | ! Do we use bands or g-points? |
---|
53 | logical :: use_bands = .false. |
---|
54 | |
---|
55 | contains |
---|
56 | procedure :: setup => setup_general_cloud_optics |
---|
57 | procedure :: add_optical_properties |
---|
58 | procedure :: save => save_general_cloud_optics_data |
---|
59 | |
---|
60 | end type general_cloud_optics_type |
---|
61 | |
---|
62 | contains |
---|
63 | |
---|
64 | ! Provides elemental function "delta_eddington" |
---|
65 | #include "radiation_delta_eddington.h" |
---|
66 | |
---|
67 | !--------------------------------------------------------------------- |
---|
68 | ! Setup cloud optics coefficients by reading them from a file |
---|
69 | subroutine setup_general_cloud_optics(this, file_name, specdef, & |
---|
70 | & use_bands, use_thick_averaging, & |
---|
71 | & weighting_temperature, & |
---|
72 | & iverbose) |
---|
73 | |
---|
74 | use yomhook, only : lhook, dr_hook, jphook |
---|
75 | use easy_netcdf, only : netcdf_file |
---|
76 | use radiation_spectral_definition, only : spectral_definition_type |
---|
77 | use radiation_io, only : nulout, nulerr, radiation_abort |
---|
78 | |
---|
79 | class(general_cloud_optics_type), intent(inout) :: this |
---|
80 | character(len=*), intent(in) :: file_name |
---|
81 | type(spectral_definition_type), intent(in) :: specdef |
---|
82 | logical, intent(in), optional :: use_bands, use_thick_averaging |
---|
83 | real(jprb), intent(in), optional :: weighting_temperature ! K |
---|
84 | integer, intent(in), optional :: iverbose |
---|
85 | |
---|
86 | ! Spectral properties read from file, dimensioned (wavenumber, |
---|
87 | ! n_effective_radius) |
---|
88 | real(jprb), dimension(:,:), allocatable :: mass_ext, & ! m2 kg-1 |
---|
89 | & ssa, asymmetry |
---|
90 | |
---|
91 | ! Reflectance of an infinitely thick cloud, needed for thick |
---|
92 | ! averaging |
---|
93 | real(jprb), dimension(:,:), allocatable :: ref_inf |
---|
94 | |
---|
95 | ! Coordinate variables from file |
---|
96 | real(jprb), dimension(:), allocatable :: wavenumber ! cm-1 |
---|
97 | real(jprb), dimension(:), allocatable :: effective_radius ! m |
---|
98 | |
---|
99 | ! Matrix mapping optical properties in the file to values per |
---|
100 | ! g-point or band, such that in the thin-averaging case, |
---|
101 | ! this%mass_ext=matmul(mapping,file%mass_ext), so mapping is |
---|
102 | ! dimensioned (ngpoint,nwav) |
---|
103 | real(jprb), dimension(:,:), allocatable :: mapping |
---|
104 | |
---|
105 | ! The NetCDF file containing the coefficients |
---|
106 | type(netcdf_file) :: file |
---|
107 | |
---|
108 | real(jprb) :: diff_spread |
---|
109 | integer :: iverb |
---|
110 | integer :: nre ! Number of effective radii |
---|
111 | integer :: nwav ! Number of wavenumbers describing cloud |
---|
112 | |
---|
113 | logical :: use_bands_local, use_thick_averaging_local |
---|
114 | |
---|
115 | real(jphook) :: hook_handle |
---|
116 | |
---|
117 | if (lhook) call dr_hook('radiation_general_cloud_optics_data:setup',0,hook_handle) |
---|
118 | |
---|
119 | ! Set local values of optional inputs |
---|
120 | if (present(iverbose)) then |
---|
121 | iverb = iverbose |
---|
122 | else |
---|
123 | iverb = 2 |
---|
124 | end if |
---|
125 | |
---|
126 | if (present(use_bands)) then |
---|
127 | use_bands_local = use_bands |
---|
128 | else |
---|
129 | use_bands_local = .false. |
---|
130 | end if |
---|
131 | |
---|
132 | if (present(use_thick_averaging)) then |
---|
133 | use_thick_averaging_local = use_thick_averaging |
---|
134 | else |
---|
135 | use_thick_averaging_local = .false. |
---|
136 | end if |
---|
137 | |
---|
138 | ! Open the scattering file and configure the way it is read |
---|
139 | call file%open(trim(file_name), iverbose=iverb) |
---|
140 | !call file%transpose_matrices() |
---|
141 | |
---|
142 | ! Read coordinate variables |
---|
143 | call file%get('wavenumber', wavenumber) |
---|
144 | call file%get('effective_radius', effective_radius) |
---|
145 | |
---|
146 | ! Read the band-specific coefficients |
---|
147 | call file%get('mass_extinction_coefficient', mass_ext) |
---|
148 | call file%get('single_scattering_albedo', ssa) |
---|
149 | call file%get('asymmetry_factor', asymmetry) |
---|
150 | |
---|
151 | ! Close scattering file |
---|
152 | call file%close() |
---|
153 | |
---|
154 | ! Check effective radius is evenly spaced |
---|
155 | nre = size(effective_radius) |
---|
156 | ! Fractional range of differences, should be near zero for evenly |
---|
157 | ! spaced data |
---|
158 | diff_spread = (maxval(effective_radius(2:nre)-effective_radius(1:nre-1)) & |
---|
159 | & -minval(effective_radius(2:nre)-effective_radius(1:nre-1))) & |
---|
160 | & / minval(abs(effective_radius(2:nre)-effective_radius(1:nre-1))) |
---|
161 | if (diff_spread > 0.01_jprb) then |
---|
162 | write(nulerr, '(a,a,a)') '*** Error: effective_radius in ', & |
---|
163 | & trim(file_name), ', is not evenly spaced to 1%' |
---|
164 | call radiation_abort('Radiation configuration error') |
---|
165 | end if |
---|
166 | |
---|
167 | ! Set up effective radius coordinate variable |
---|
168 | this%n_effective_radius = nre |
---|
169 | this%effective_radius_0 = effective_radius(1) |
---|
170 | this%d_effective_radius = effective_radius(2) - effective_radius(1) |
---|
171 | |
---|
172 | nwav = size(wavenumber) |
---|
173 | |
---|
174 | ! Define the mapping matrix |
---|
175 | call specdef%calc_mapping(wavenumber, mapping, & |
---|
176 | weighting_temperature=weighting_temperature, use_bands=use_bands) |
---|
177 | |
---|
178 | ! Thick averaging should be performed on delta-Eddington scaled |
---|
179 | ! quantities (it makes no difference to thin averaging) |
---|
180 | call delta_eddington(mass_ext, ssa, asymmetry) |
---|
181 | |
---|
182 | ! Thin averaging |
---|
183 | this%mass_ext = matmul(mapping, mass_ext) |
---|
184 | this%ssa = matmul(mapping, mass_ext*ssa) / this%mass_ext |
---|
185 | this%asymmetry = matmul(mapping, mass_ext*ssa*asymmetry) / (this%mass_ext*this%ssa) |
---|
186 | |
---|
187 | if (use_thick_averaging_local) then |
---|
188 | ! Thick averaging as described by Edwards and Slingo (1996), |
---|
189 | ! modifying only the single-scattering albedo |
---|
190 | allocate(ref_inf(nwav, nre)) |
---|
191 | |
---|
192 | ! Eqs. 18 and 17 of Edwards & Slingo (1996) |
---|
193 | ref_inf = sqrt((1.0_jprb - ssa) / (1.0_jprb - ssa*asymmetry)) |
---|
194 | ref_inf = (1.0_jprb - ref_inf) / (1.0_jprb + ref_inf) |
---|
195 | ! Here the left-hand side is actually the averaged ref_inf |
---|
196 | this%ssa = matmul(mapping, ref_inf) |
---|
197 | ! Eq. 19 of Edwards and Slingo (1996) |
---|
198 | this%ssa = 4.0_jprb * this%ssa / ((1.0_jprb + this%ssa)**2 & |
---|
199 | & - this%asymmetry * (1.0_jprb - this%ssa)**2) |
---|
200 | |
---|
201 | deallocate(ref_inf) |
---|
202 | end if |
---|
203 | |
---|
204 | deallocate(mapping) |
---|
205 | |
---|
206 | ! Revert back to unscaled quantities |
---|
207 | call revert_delta_eddington(this%mass_ext, this%ssa, this%asymmetry) |
---|
208 | |
---|
209 | if (iverb >= 2) then |
---|
210 | write(nulout,'(a,a)') ' File: ', trim(file_name) |
---|
211 | if (present(weighting_temperature)) then |
---|
212 | write(nulout,'(a,f7.1,a)') ' Weighting temperature: ', weighting_temperature, ' K' |
---|
213 | else |
---|
214 | write(nulout,'(a,f7.1,a)') ' Weighting temperature: ', specdef%reference_temperature, ' K' |
---|
215 | end if |
---|
216 | if (use_thick_averaging_local) then |
---|
217 | write(nulout,'(a)') ' SSA averaging: optically thick limit' |
---|
218 | else |
---|
219 | write(nulout,'(a)') ' SSA averaging: optically thin limit' |
---|
220 | end if |
---|
221 | if (use_bands_local) then |
---|
222 | write(nulout,'(a,i0,a)') ' Spectral discretization: ', specdef%nband, ' bands' |
---|
223 | else |
---|
224 | write(nulout,'(a,i0,a)') ' Spectral discretization: ', specdef%ng, ' g-points' |
---|
225 | end if |
---|
226 | write(nulout,'(a,i0,a,f6.1,a,f6.1,a)') ' Effective radius look-up: ', nre, ' points in range ', & |
---|
227 | & effective_radius(1)*1.0e6_jprb, '-', effective_radius(nre)*1.0e6_jprb, ' um' |
---|
228 | write(nulout,'(a,i0,a,i0,a)') ' Wavenumber range: ', int(specdef%min_wavenumber()), '-', & |
---|
229 | & int(specdef%max_wavenumber()), ' cm-1' |
---|
230 | end if |
---|
231 | |
---|
232 | if (lhook) call dr_hook('radiation_general_cloud_optics_data:setup',1,hook_handle) |
---|
233 | |
---|
234 | end subroutine setup_general_cloud_optics |
---|
235 | |
---|
236 | |
---|
237 | !--------------------------------------------------------------------- |
---|
238 | ! Add the optical properties of a particular cloud type to the |
---|
239 | ! accumulated optical properties of all cloud types |
---|
240 | subroutine add_optical_properties(this, ng, nlev, ncol, & |
---|
241 | & cloud_fraction, & |
---|
242 | & water_path, effective_radius, & |
---|
243 | & od, scat_od, scat_asymmetry) |
---|
244 | |
---|
245 | use yomhook, only : lhook, dr_hook, jphook |
---|
246 | |
---|
247 | class(general_cloud_optics_type), intent(in) :: this |
---|
248 | |
---|
249 | ! Number of g points, levels and columns |
---|
250 | integer, intent(in) :: ng, nlev, ncol |
---|
251 | |
---|
252 | ! Properties of present cloud type, dimensioned (ncol,nlev) |
---|
253 | real(jprb), intent(in) :: cloud_fraction(:,:) |
---|
254 | real(jprb), intent(in) :: water_path(:,:) ! kg m-2 |
---|
255 | real(jprb), intent(in) :: effective_radius(:,:) ! m |
---|
256 | |
---|
257 | ! Optical properties which are additive per cloud type, |
---|
258 | ! dimensioned (ng,nlev,ncol) |
---|
259 | real(jprb), intent(inout), dimension(ng,nlev,ncol) & |
---|
260 | & :: od ! Optical depth of layer |
---|
261 | real(jprb), intent(inout), dimension(ng,nlev,ncol), optional & |
---|
262 | & :: scat_od, & ! Scattering optical depth of layer |
---|
263 | & scat_asymmetry ! Scattering optical depth x asymmetry factor |
---|
264 | |
---|
265 | real(jprb) :: od_local(ng) |
---|
266 | |
---|
267 | real(jprb) :: re_index, weight1, weight2 |
---|
268 | integer :: ire |
---|
269 | |
---|
270 | integer :: jcol, jlev |
---|
271 | |
---|
272 | real(jphook) :: hook_handle |
---|
273 | |
---|
274 | if (lhook) call dr_hook('radiation_general_cloud_optics_data:add_optical_properties',0,hook_handle) |
---|
275 | |
---|
276 | if (present(scat_od)) then |
---|
277 | do jcol = 1,ncol |
---|
278 | do jlev = 1,nlev |
---|
279 | if (cloud_fraction(jcol, jlev) > 0.0_jprb) then |
---|
280 | re_index = max(1.0_jprb, min(1.0_jprb + (effective_radius(jcol,jlev)-this%effective_radius_0) & |
---|
281 | & / this%d_effective_radius, this%n_effective_radius-0.0001_jprb)) |
---|
282 | ire = int(re_index) |
---|
283 | weight2 = re_index - ire |
---|
284 | weight1 = 1.0_jprb - weight2 |
---|
285 | od_local = water_path(jcol, jlev) * (weight1*this%mass_ext(:,ire) & |
---|
286 | & +weight2*this%mass_ext(:,ire+1)) |
---|
287 | od(:,jlev,jcol) = od(:,jlev,jcol) + od_local |
---|
288 | od_local = od_local * (weight1*this%ssa(:,ire) & |
---|
289 | & +weight2*this%ssa(:,ire+1)) |
---|
290 | scat_od(:,jlev,jcol) = scat_od(:,jlev,jcol) + od_local |
---|
291 | scat_asymmetry(:,jlev,jcol) = scat_asymmetry(:,jlev,jcol) & |
---|
292 | & + od_local * (weight1*this%asymmetry(:,ire) & |
---|
293 | & +weight2*this%asymmetry(:,ire+1)) |
---|
294 | end if |
---|
295 | end do |
---|
296 | end do |
---|
297 | else |
---|
298 | ! No scattering: return the absorption optical depth |
---|
299 | do jcol = 1,ncol |
---|
300 | do jlev = 1,nlev |
---|
301 | if (water_path(jcol, jlev) > 0.0_jprb) then |
---|
302 | re_index = max(1.0_jprb, min(1.0_jprb + (effective_radius(jcol,jlev)-this%effective_radius_0) & |
---|
303 | & / this%d_effective_radius, this%n_effective_radius-0.0001_jprb)) |
---|
304 | ire = int(re_index) |
---|
305 | weight2 = re_index - ire |
---|
306 | weight1 = 1.0_jprb - weight2 |
---|
307 | od(:,jlev,jcol) = od(:,jlev,jcol) & |
---|
308 | & + water_path(jcol, jlev) * (weight1*this%mass_ext(:,ire) & |
---|
309 | & +weight2*this%mass_ext(:,ire+1)) & |
---|
310 | & * (1.0_jprb - (weight1*this%ssa(:,ire)+weight2*this%ssa(:,ire+1))) |
---|
311 | end if |
---|
312 | end do |
---|
313 | end do |
---|
314 | end if |
---|
315 | |
---|
316 | if (lhook) call dr_hook('radiation_general_cloud_optics_data:add_optical_properties',1,hook_handle) |
---|
317 | |
---|
318 | end subroutine add_optical_properties |
---|
319 | |
---|
320 | |
---|
321 | !--------------------------------------------------------------------- |
---|
322 | ! Return the Planck function (in W m-2 (cm-1)-1) for a given |
---|
323 | ! wavenumber (cm-1) and temperature (K), ensuring double precision |
---|
324 | ! for internal calculation |
---|
325 | elemental function calc_planck_function_wavenumber(wavenumber, temperature) |
---|
326 | |
---|
327 | use parkind1, only : jprb, jprd |
---|
328 | use radiation_constants, only : SpeedOfLight, BoltzmannConstant, PlanckConstant |
---|
329 | |
---|
330 | real(jprb), intent(in) :: wavenumber ! cm-1 |
---|
331 | real(jprb), intent(in) :: temperature ! K |
---|
332 | real(jprb) :: calc_planck_function_wavenumber |
---|
333 | |
---|
334 | real(jprd) :: freq ! Hz |
---|
335 | real(jprd) :: planck_fn_freq ! W m-2 Hz-1 |
---|
336 | |
---|
337 | freq = 100.0_jprd * real(SpeedOfLight,jprd) * real(wavenumber,jprd) |
---|
338 | planck_fn_freq = 2.0_jprd * real(PlanckConstant,jprd) * freq**3 & |
---|
339 | & / (real(SpeedOfLight,jprd)**2 * (exp(real(PlanckConstant,jprd)*freq & |
---|
340 | & /(real(BoltzmannConstant,jprd)*real(temperature,jprd))) - 1.0_jprd)) |
---|
341 | calc_planck_function_wavenumber = real(planck_fn_freq * 100.0_jprd * real(SpeedOfLight,jprd), jprb) |
---|
342 | |
---|
343 | end function calc_planck_function_wavenumber |
---|
344 | |
---|
345 | !--------------------------------------------------------------------- |
---|
346 | ! Save cloud optical properties in the named file |
---|
347 | subroutine save_general_cloud_optics_data(this, file_name, iverbose) |
---|
348 | |
---|
349 | use yomhook, only : lhook, dr_hook, jphook |
---|
350 | use easy_netcdf, only : netcdf_file |
---|
351 | |
---|
352 | class(general_cloud_optics_type), intent(in) :: this |
---|
353 | character(len=*), intent(in) :: file_name |
---|
354 | integer, optional, intent(in) :: iverbose |
---|
355 | |
---|
356 | ! Object for output NetCDF file |
---|
357 | type(netcdf_file) :: out_file |
---|
358 | |
---|
359 | real(jprb) :: effective_radius(this%n_effective_radius) |
---|
360 | integer :: ire |
---|
361 | |
---|
362 | real(jphook) :: hook_handle |
---|
363 | |
---|
364 | if (lhook) call dr_hook('radiation_general_cloud_optics_data:save',0,hook_handle) |
---|
365 | |
---|
366 | ! Create the file |
---|
367 | call out_file%create(trim(file_name), iverbose=iverbose) |
---|
368 | |
---|
369 | ! Define dimensions |
---|
370 | call out_file%define_dimension("band", size(this%mass_ext,1)) |
---|
371 | call out_file%define_dimension("effective_radius", this%n_effective_radius) |
---|
372 | |
---|
373 | ! Put global attributes |
---|
374 | call out_file%put_global_attributes( & |
---|
375 | & title_str="Optical properties of "//trim(this%type_name) & |
---|
376 | & //" hydrometeors using the spectral intervals of ecRad", & |
---|
377 | & source_str="ecRad offline radiation model") |
---|
378 | |
---|
379 | ! Define variables |
---|
380 | call out_file%define_variable("effective_radius", units_str="m", & |
---|
381 | & long_name="Effective radius", dim1_name="effective_radius") |
---|
382 | call out_file%define_variable("mass_extinction_coefficient", units_str="m2 kg-1", & |
---|
383 | & long_name="Mass-extinction coefficient", & |
---|
384 | & dim2_name="effective_radius", dim1_name="band") |
---|
385 | call out_file%define_variable("single_scattering_albedo", units_str="1", & |
---|
386 | & long_name="Single scattering albedo", & |
---|
387 | & dim2_name="effective_radius", dim1_name="band") |
---|
388 | call out_file%define_variable("asymmetry_factor", units_str="1", & |
---|
389 | & long_name="Asymmetry factor", & |
---|
390 | & dim2_name="effective_radius", dim1_name="band") |
---|
391 | |
---|
392 | ! Define effective radius |
---|
393 | do ire = 1,this%n_effective_radius |
---|
394 | effective_radius(ire) = this%effective_radius_0 + this%d_effective_radius*(ire-1) |
---|
395 | end do |
---|
396 | |
---|
397 | ! Write variables |
---|
398 | call out_file%put("effective_radius", effective_radius) |
---|
399 | call out_file%put("mass_extinction_coefficient", this%mass_ext) |
---|
400 | call out_file%put("single_scattering_albedo", this%ssa) |
---|
401 | call out_file%put("asymmetry_factor", this%asymmetry) |
---|
402 | |
---|
403 | call out_file%close() |
---|
404 | |
---|
405 | if (lhook) call dr_hook('radiation_general_cloud_optics_data:save',1,hook_handle) |
---|
406 | |
---|
407 | end subroutine save_general_cloud_optics_data |
---|
408 | |
---|
409 | end module radiation_general_cloud_optics_data |
---|