| 1 | ! radiation_interface.F90 - Monochromatic gas/cloud optics for testing |
|---|
| 2 | |
|---|
| 3 | ! (C) Copyright 2014- 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 | |
|---|
| 15 | ! Modifications |
|---|
| 16 | ! 2017-04-11 R. Hogan Receive "surface" dummy argument |
|---|
| 17 | ! 2017-09-13 R. Hogan Revert |
|---|
| 18 | ! 2018-08-29 R. Hogan Particulate single-scattering albedo / asymmetry from namelist |
|---|
| 19 | |
|---|
| 20 | module radiation_monochromatic |
|---|
| 21 | |
|---|
| 22 | implicit none |
|---|
| 23 | |
|---|
| 24 | public :: setup_gas_optics, gas_optics, set_gas_units, & |
|---|
| 25 | & setup_cloud_optics, cloud_optics, & |
|---|
| 26 | & setup_aerosol_optics, add_aerosol_optics |
|---|
| 27 | |
|---|
| 28 | contains |
|---|
| 29 | |
|---|
| 30 | ! Provides elemental function "delta_eddington" |
|---|
| 31 | #include "radiation_delta_eddington.h" |
|---|
| 32 | |
|---|
| 33 | !--------------------------------------------------------------------- |
|---|
| 34 | ! Setup the arrays in the config object corresponding to the |
|---|
| 35 | ! monochromatic gas optics model. The directory argument is not |
|---|
| 36 | ! used, since no look-up tables need to be loaded. |
|---|
| 37 | subroutine setup_gas_optics(config, directory) |
|---|
| 38 | |
|---|
| 39 | use radiation_config, only : config_type |
|---|
| 40 | |
|---|
| 41 | type(config_type), intent(inout) :: config |
|---|
| 42 | character(len=*), intent(in) :: directory |
|---|
| 43 | |
|---|
| 44 | ! In the monochromatic model we have simply one band and g-point |
|---|
| 45 | ! in both the longwave and shortwave parts of the spectrum |
|---|
| 46 | config%n_g_sw = 1 |
|---|
| 47 | config%n_g_lw = 1 |
|---|
| 48 | config%n_bands_sw = 1 |
|---|
| 49 | config%n_bands_lw = 1 |
|---|
| 50 | |
|---|
| 51 | ! Allocate arrays |
|---|
| 52 | allocate(config%i_band_from_g_sw (config%n_g_sw)) |
|---|
| 53 | allocate(config%i_band_from_g_lw (config%n_g_lw)) |
|---|
| 54 | allocate(config%i_band_from_reordered_g_sw(config%n_g_sw)) |
|---|
| 55 | allocate(config%i_band_from_reordered_g_lw(config%n_g_lw)) |
|---|
| 56 | allocate(config%i_g_from_reordered_g_sw(config%n_g_sw)) |
|---|
| 57 | allocate(config%i_g_from_reordered_g_lw(config%n_g_lw)) |
|---|
| 58 | |
|---|
| 59 | ! Indices are trivial... |
|---|
| 60 | config%i_band_from_g_sw = 1 |
|---|
| 61 | config%i_band_from_g_lw = 1 |
|---|
| 62 | config%i_g_from_reordered_g_sw = 1 |
|---|
| 63 | config%i_g_from_reordered_g_lw = 1 |
|---|
| 64 | config%i_band_from_reordered_g_sw = 1 |
|---|
| 65 | config%i_band_from_reordered_g_lw = 1 |
|---|
| 66 | |
|---|
| 67 | end subroutine setup_gas_optics |
|---|
| 68 | |
|---|
| 69 | |
|---|
| 70 | !--------------------------------------------------------------------- |
|---|
| 71 | ! Dummy routine for scaling gas mixing ratios |
|---|
| 72 | subroutine set_gas_units(gas) |
|---|
| 73 | |
|---|
| 74 | use radiation_gas, only : gas_type |
|---|
| 75 | type(gas_type), intent(inout) :: gas |
|---|
| 76 | |
|---|
| 77 | end subroutine set_gas_units |
|---|
| 78 | |
|---|
| 79 | |
|---|
| 80 | !--------------------------------------------------------------------- |
|---|
| 81 | ! Dummy setup routine for cloud optics: in fact, no setup is |
|---|
| 82 | ! required for monochromatic case |
|---|
| 83 | subroutine setup_cloud_optics(config) |
|---|
| 84 | |
|---|
| 85 | use radiation_config, only : config_type |
|---|
| 86 | type(config_type), intent(inout) :: config |
|---|
| 87 | |
|---|
| 88 | end subroutine setup_cloud_optics |
|---|
| 89 | |
|---|
| 90 | |
|---|
| 91 | !--------------------------------------------------------------------- |
|---|
| 92 | ! Dummy subroutine since no aerosols are represented in |
|---|
| 93 | ! monochromatic case |
|---|
| 94 | subroutine setup_aerosol_optics(config) |
|---|
| 95 | |
|---|
| 96 | use radiation_config, only : config_type |
|---|
| 97 | type(config_type), intent(inout) :: config |
|---|
| 98 | |
|---|
| 99 | end subroutine setup_aerosol_optics |
|---|
| 100 | |
|---|
| 101 | |
|---|
| 102 | !--------------------------------------------------------------------- |
|---|
| 103 | ! Compute gas optical depths, shortwave scattering, Planck function |
|---|
| 104 | ! and incoming shortwave radiation at top-of-atmosphere |
|---|
| 105 | subroutine gas_optics(ncol,nlev,istartcol,iendcol, & |
|---|
| 106 | config, single_level, thermodynamics, gas, lw_albedo, & |
|---|
| 107 | od_lw, od_sw, ssa_sw, planck_hl, lw_emission, & |
|---|
| 108 | incoming_sw) |
|---|
| 109 | |
|---|
| 110 | use parkind1, only : jprb |
|---|
| 111 | use radiation_config, only : config_type |
|---|
| 112 | use radiation_thermodynamics, only : thermodynamics_type |
|---|
| 113 | use radiation_single_level, only : single_level_type |
|---|
| 114 | use radiation_gas, only : gas_type |
|---|
| 115 | use radiation_constants, only : Pi, StefanBoltzmann |
|---|
| 116 | |
|---|
| 117 | ! Inputs |
|---|
| 118 | integer, intent(in) :: ncol ! number of columns |
|---|
| 119 | integer, intent(in) :: nlev ! number of model levels |
|---|
| 120 | integer, intent(in) :: istartcol, iendcol ! range of columns to process |
|---|
| 121 | type(config_type), intent(in) :: config |
|---|
| 122 | type(single_level_type), intent(in) :: single_level |
|---|
| 123 | type(thermodynamics_type),intent(in) :: thermodynamics |
|---|
| 124 | type(gas_type), intent(in) :: gas |
|---|
| 125 | |
|---|
| 126 | ! Longwave albedo of the surface |
|---|
| 127 | real(jprb), dimension(config%n_g_lw,istartcol:iendcol), & |
|---|
| 128 | & intent(in) :: lw_albedo |
|---|
| 129 | |
|---|
| 130 | ! Outputs |
|---|
| 131 | |
|---|
| 132 | ! Gaseous layer optical depth in longwave and shortwave, and |
|---|
| 133 | ! shortwave single scattering albedo (i.e. fraction of extinction |
|---|
| 134 | ! due to Rayleigh scattering) at each g-point |
|---|
| 135 | real(jprb), dimension(config%n_g_lw,nlev,istartcol:iendcol), intent(out) :: & |
|---|
| 136 | & od_lw |
|---|
| 137 | real(jprb), dimension(config%n_g_sw,nlev,istartcol:iendcol), intent(out) :: & |
|---|
| 138 | & od_sw, ssa_sw |
|---|
| 139 | |
|---|
| 140 | ! The Planck function (emitted flux from a black body) at half |
|---|
| 141 | ! levels and at the surface at each longwave g-point |
|---|
| 142 | real(jprb), dimension(config%n_g_lw,nlev+1,istartcol:iendcol), intent(out) :: & |
|---|
| 143 | & planck_hl |
|---|
| 144 | real(jprb), dimension(config%n_g_lw,istartcol:iendcol), intent(out) :: & |
|---|
| 145 | & lw_emission |
|---|
| 146 | |
|---|
| 147 | ! The incoming shortwave flux into a plane perpendicular to the |
|---|
| 148 | ! incoming radiation at top-of-atmosphere in each of the shortwave |
|---|
| 149 | ! g-points |
|---|
| 150 | real(jprb), dimension(config%n_g_sw,istartcol:iendcol), intent(out) :: & |
|---|
| 151 | & incoming_sw |
|---|
| 152 | |
|---|
| 153 | ! Ratio of the optical depth of the entire atmospheric column that |
|---|
| 154 | ! is in the current layer |
|---|
| 155 | real(jprb), dimension(istartcol:iendcol) :: extinction_fraction |
|---|
| 156 | |
|---|
| 157 | ! In the monochromatic model, the absorption by the atmosphere is |
|---|
| 158 | ! assumed proportional to the mass in each layer, so is defined in |
|---|
| 159 | ! terms of a total zenith optical depth and then distributed with |
|---|
| 160 | ! height according to the pressure. |
|---|
| 161 | !real(jprb), parameter :: total_od_sw = 0.10536_jprb ! Transmittance 0.9 |
|---|
| 162 | !real(jprb), parameter :: total_od_lw = 1.6094_jprb ! Transmittance 0.2 |
|---|
| 163 | |
|---|
| 164 | integer :: jlev |
|---|
| 165 | |
|---|
| 166 | DO jlev = 1,nlev |
|---|
| 167 | ! The fraction of the total optical depth in the current layer |
|---|
| 168 | ! is proportional to the fraction of the mass of the atmosphere |
|---|
| 169 | ! in the current layer, computed from pressure assuming |
|---|
| 170 | ! hydrostatic balance |
|---|
| 171 | extinction_fraction = & |
|---|
| 172 | & (thermodynamics%pressure_hl(istartcol:iendcol,jlev+1) & |
|---|
| 173 | & -thermodynamics%pressure_hl(istartcol:iendcol,jlev)) & |
|---|
| 174 | & /thermodynamics%pressure_hl(istartcol:iendcol,nlev) |
|---|
| 175 | |
|---|
| 176 | ! Assign longwave and shortwave optical depths |
|---|
| 177 | od_lw(1,jlev,:) = config%mono_lw_total_od*extinction_fraction |
|---|
| 178 | od_sw(1,jlev,:) = config%mono_sw_total_od*extinction_fraction |
|---|
| 179 | end do |
|---|
| 180 | |
|---|
| 181 | ! Shortwave single-scattering albedo is almost entirely Rayleigh |
|---|
| 182 | ! scattering |
|---|
| 183 | ssa_sw = 0.999999_jprb |
|---|
| 184 | |
|---|
| 185 | ! Entire shortwave spectrum represented in one band |
|---|
| 186 | incoming_sw(1,:) = single_level%solar_irradiance |
|---|
| 187 | |
|---|
| 188 | if (single_level%is_simple_surface) then |
|---|
| 189 | if (config%mono_lw_wavelength <= 0.0_jprb) then |
|---|
| 190 | ! Entire longwave spectrum represented in one band |
|---|
| 191 | lw_emission(1,istartcol:iendcol) & |
|---|
| 192 | & = StefanBoltzmann * single_level%skin_temperature(istartcol:iendcol)**4 & |
|---|
| 193 | & * single_level%lw_emissivity(istartcol:iendcol,1) |
|---|
| 194 | DO jlev = 1,nlev+1 |
|---|
| 195 | planck_hl(1,jlev,istartcol:iendcol) = StefanBoltzmann * thermodynamics%temperature_hl(istartcol:iendcol,jlev)**4 |
|---|
| 196 | end do |
|---|
| 197 | else |
|---|
| 198 | ! Single wavelength: multiply by pi to convert W sr-1 m-3 to W m-3 |
|---|
| 199 | lw_emission(1,istartcol:iendcol) = Pi*planck_function(config%mono_lw_wavelength, & |
|---|
| 200 | & single_level%skin_temperature(istartcol:iendcol)) & |
|---|
| 201 | & * single_level%lw_emissivity(istartcol:iendcol,1) |
|---|
| 202 | DO jlev = 1,nlev+1 |
|---|
| 203 | planck_hl(1,jlev,istartcol:iendcol) = Pi*planck_function(config%mono_lw_wavelength, & |
|---|
| 204 | & thermodynamics%temperature_hl(istartcol:iendcol,jlev)) |
|---|
| 205 | end do |
|---|
| 206 | end if |
|---|
| 207 | else |
|---|
| 208 | lw_emission = transpose(single_level%lw_emission) |
|---|
| 209 | end if |
|---|
| 210 | |
|---|
| 211 | end subroutine gas_optics |
|---|
| 212 | |
|---|
| 213 | |
|---|
| 214 | !--------------------------------------------------------------------- |
|---|
| 215 | ! Compute cloud optical depth, single-scattering albedo and |
|---|
| 216 | ! g factor in the longwave and shortwave |
|---|
| 217 | subroutine cloud_optics(nlev,istartcol,iendcol, & |
|---|
| 218 | & config, thermodynamics, cloud, & |
|---|
| 219 | & od_lw_cloud, ssa_lw_cloud, g_lw_cloud, & |
|---|
| 220 | & od_sw_cloud, ssa_sw_cloud, g_sw_cloud) |
|---|
| 221 | |
|---|
| 222 | use parkind1, only : jprb |
|---|
| 223 | use radiation_config, only : config_type |
|---|
| 224 | use radiation_thermodynamics, only : thermodynamics_type |
|---|
| 225 | use radiation_cloud, only : cloud_type |
|---|
| 226 | use radiation_constants, only : AccelDueToGravity, & |
|---|
| 227 | & DensityLiquidWater, DensitySolidIce |
|---|
| 228 | |
|---|
| 229 | ! Inputs |
|---|
| 230 | integer, intent(in) :: nlev ! number of model levels |
|---|
| 231 | integer, intent(in) :: istartcol, iendcol ! range of columns to process |
|---|
| 232 | type(config_type), intent(in) :: config |
|---|
| 233 | type(thermodynamics_type),intent(in) :: thermodynamics |
|---|
| 234 | type(cloud_type), intent(in) :: cloud |
|---|
| 235 | |
|---|
| 236 | ! Outputs |
|---|
| 237 | |
|---|
| 238 | ! Layer optical depth, single scattering albedo and g factor of |
|---|
| 239 | ! clouds in each longwave band, where the latter two |
|---|
| 240 | ! variables are only defined if cloud longwave scattering is |
|---|
| 241 | ! enabled (otherwise both are treated as zero). |
|---|
| 242 | real(jprb), dimension(config%n_bands_lw,nlev,istartcol:iendcol), intent(out) :: & |
|---|
| 243 | & od_lw_cloud |
|---|
| 244 | real(jprb), dimension(config%n_bands_lw_if_scattering,nlev,istartcol:iendcol), & |
|---|
| 245 | & intent(out) :: ssa_lw_cloud, g_lw_cloud |
|---|
| 246 | |
|---|
| 247 | ! Layer optical depth, single scattering albedo and g factor of |
|---|
| 248 | ! clouds in each shortwave band |
|---|
| 249 | real(jprb), dimension(config%n_g_sw,nlev,istartcol:iendcol), intent(out) :: & |
|---|
| 250 | & od_sw_cloud, ssa_sw_cloud, g_sw_cloud |
|---|
| 251 | |
|---|
| 252 | ! In-cloud liquid and ice water path in a layer, in kg m-2 |
|---|
| 253 | real(jprb), dimension(nlev,istartcol:iendcol) :: lwp_kg_m2, iwp_kg_m2 |
|---|
| 254 | |
|---|
| 255 | integer :: jlev, jcol |
|---|
| 256 | |
|---|
| 257 | ! Factor to convert from gridbox-mean mass mixing ratio to |
|---|
| 258 | ! in-cloud water path |
|---|
| 259 | real(jprb) :: factor |
|---|
| 260 | |
|---|
| 261 | ! Convert cloud mixing ratio into liquid and ice water path |
|---|
| 262 | ! in each layer |
|---|
| 263 | DO jlev = 1, nlev |
|---|
| 264 | DO jcol = istartcol, iendcol |
|---|
| 265 | ! Factor to convert from gridbox-mean mass mixing ratio to |
|---|
| 266 | ! in-cloud water path involves the pressure difference in |
|---|
| 267 | ! Pa, acceleration due to gravity and cloud fraction |
|---|
| 268 | ! adjusted to avoid division by zero. |
|---|
| 269 | factor = ( thermodynamics%pressure_hl(jcol,jlev+1) & |
|---|
| 270 | & -thermodynamics%pressure_hl(jcol,jlev ) ) & |
|---|
| 271 | & / (AccelDueToGravity & |
|---|
| 272 | & * max(epsilon(1.0_jprb), cloud%fraction(jcol,jlev))) |
|---|
| 273 | lwp_kg_m2(jlev,jcol) = factor * cloud%q_liq(jcol,jlev) |
|---|
| 274 | iwp_kg_m2(jlev,jcol) = factor * cloud%q_ice(jcol,jlev) |
|---|
| 275 | end do |
|---|
| 276 | end do |
|---|
| 277 | |
|---|
| 278 | ! Geometric optics approximation: particles treated as much larger |
|---|
| 279 | ! than the wavelength in both shortwave and longwave |
|---|
| 280 | od_sw_cloud(1,:,:) & |
|---|
| 281 | & = (3.0_jprb/(2.0_jprb*DensityLiquidWater)) & |
|---|
| 282 | & * lwp_kg_m2 / transpose(cloud%re_liq(istartcol:iendcol,:)) & |
|---|
| 283 | & + (3.0_jprb / (2.0_jprb * DensitySolidIce)) & |
|---|
| 284 | & * iwp_kg_m2 / transpose(cloud%re_ice(istartcol:iendcol,:)) |
|---|
| 285 | od_lw_cloud(1,:,:) = lwp_kg_m2 * 137.22_jprb & |
|---|
| 286 | & + (3.0_jprb / (2.0_jprb * DensitySolidIce)) & |
|---|
| 287 | & * iwp_kg_m2 / transpose(cloud%re_ice(istartcol:iendcol,:)) |
|---|
| 288 | |
|---|
| 289 | if (config%iverbose >= 4) then |
|---|
| 290 | DO jcol = istartcol,iendcol |
|---|
| 291 | write(*,'(a,i0,a,f7.3,a,f7.3)') 'Profile ', jcol, ': shortwave optical depth = ', & |
|---|
| 292 | & sum(od_sw_cloud(1,:,jcol)*cloud%fraction(jcol,:)), & |
|---|
| 293 | & ', longwave optical depth = ', & |
|---|
| 294 | & sum(od_lw_cloud(1,:,jcol)*cloud%fraction(jcol,:)) |
|---|
| 295 | ! PRINT *, 'LWP = ', sum(lwp_kg_m2(:,istartcol)*cloud%fraction(istartcol,:)) |
|---|
| 296 | end do |
|---|
| 297 | end if |
|---|
| 298 | |
|---|
| 299 | ssa_sw_cloud = config%mono_sw_single_scattering_albedo |
|---|
| 300 | g_sw_cloud = config%mono_sw_asymmetry_factor |
|---|
| 301 | |
|---|
| 302 | ! In-place delta-Eddington scaling |
|---|
| 303 | call delta_eddington(od_sw_cloud, ssa_sw_cloud, g_sw_cloud) |
|---|
| 304 | |
|---|
| 305 | if (config%do_lw_cloud_scattering) then |
|---|
| 306 | ssa_lw_cloud = config%mono_lw_single_scattering_albedo |
|---|
| 307 | g_lw_cloud = config%mono_lw_asymmetry_factor |
|---|
| 308 | ! In-place delta-Eddington scaling |
|---|
| 309 | call delta_eddington(od_lw_cloud, ssa_lw_cloud, g_lw_cloud) |
|---|
| 310 | end if |
|---|
| 311 | |
|---|
| 312 | end subroutine cloud_optics |
|---|
| 313 | |
|---|
| 314 | |
|---|
| 315 | !--------------------------------------------------------------------- |
|---|
| 316 | ! Dummy subroutine since no aerosols are represented in |
|---|
| 317 | ! monochromatic case |
|---|
| 318 | subroutine add_aerosol_optics(nlev,istartcol,iendcol, & |
|---|
| 319 | & config, thermodynamics, gas, aerosol, & |
|---|
| 320 | & od_lw, ssa_lw, g_lw, od_sw, ssa_sw, g_sw) |
|---|
| 321 | |
|---|
| 322 | use parkind1, only : jprb |
|---|
| 323 | |
|---|
| 324 | use radiation_config, only : config_type |
|---|
| 325 | use radiation_thermodynamics, only : thermodynamics_type |
|---|
| 326 | use radiation_gas, only : gas_type |
|---|
| 327 | use radiation_aerosol, only : aerosol_type |
|---|
| 328 | |
|---|
| 329 | integer, intent(in) :: nlev ! number of model levels |
|---|
| 330 | integer, intent(in) :: istartcol, iendcol ! range of columns to process |
|---|
| 331 | type(config_type), intent(in), target :: config |
|---|
| 332 | type(thermodynamics_type),intent(in) :: thermodynamics |
|---|
| 333 | type(gas_type), intent(in) :: gas |
|---|
| 334 | type(aerosol_type), intent(in) :: aerosol |
|---|
| 335 | ! Optical depth, single scattering albedo and asymmetry factor of |
|---|
| 336 | ! the atmosphere (gases on input, gases and aerosols on output) |
|---|
| 337 | ! for each g point. Note that longwave ssa and asymmetry and |
|---|
| 338 | ! shortwave asymmetry are all zero for gases, so are not yet |
|---|
| 339 | ! defined on input and are therefore intent(out). |
|---|
| 340 | real(jprb), dimension(config%n_g_lw,nlev,istartcol:iendcol), intent(inout) :: od_lw |
|---|
| 341 | real(jprb), dimension(config%n_g_lw_if_scattering,nlev,istartcol:iendcol), & |
|---|
| 342 | & intent(out) :: ssa_lw, g_lw |
|---|
| 343 | real(jprb), dimension(config%n_g_sw,nlev,istartcol:iendcol), intent(inout) & |
|---|
| 344 | & :: od_sw, ssa_sw |
|---|
| 345 | real(jprb), dimension(config%n_g_sw,nlev,istartcol:iendcol), intent(out) :: g_sw |
|---|
| 346 | |
|---|
| 347 | g_sw(:,:,istartcol:iendcol) = 0.0_jprb |
|---|
| 348 | |
|---|
| 349 | if (config%do_lw_aerosol_scattering) then |
|---|
| 350 | ssa_lw(:,:,istartcol:iendcol) = 0.0_jprb |
|---|
| 351 | g_lw(:,:,istartcol:iendcol) = 0.0_jprb |
|---|
| 352 | end if |
|---|
| 353 | |
|---|
| 354 | end subroutine add_aerosol_optics |
|---|
| 355 | |
|---|
| 356 | !--------------------------------------------------------------------- |
|---|
| 357 | ! Planck function in terms of wavelength |
|---|
| 358 | elemental function planck_function(wavelength, temperature) |
|---|
| 359 | |
|---|
| 360 | use parkind1, only : jprb |
|---|
| 361 | |
|---|
| 362 | use radiation_constants, only : BoltzmannConstant, PlanckConstant, & |
|---|
| 363 | & SpeedOfLight |
|---|
| 364 | |
|---|
| 365 | real(jprb), intent(in) :: wavelength ! metres |
|---|
| 366 | real(jprb), intent(in) :: temperature ! Kelvin |
|---|
| 367 | |
|---|
| 368 | ! Output in W sr-1 m-3 |
|---|
| 369 | real(jprb) :: planck_function |
|---|
| 370 | |
|---|
| 371 | if (temperature > 0.0_jprb) then |
|---|
| 372 | planck_function = 2.0_jprb * PlanckConstant * SpeedOfLight**2 & |
|---|
| 373 | & / (wavelength**5 & |
|---|
| 374 | & * (exp(PlanckConstant*SpeedOfLight & |
|---|
| 375 | & /(wavelength*BoltzmannConstant*temperature)) - 1.0_jprb)) |
|---|
| 376 | else |
|---|
| 377 | planck_function = 0.0_jprb |
|---|
| 378 | end if |
|---|
| 379 | |
|---|
| 380 | end function planck_function |
|---|
| 381 | |
|---|
| 382 | end module radiation_monochromatic |
|---|
