| 1 | ! radiation_ecckd.F90 - ecCKD generalized gas optics model |
|---|
| 2 | ! |
|---|
| 3 | ! (C) Copyright 2020- 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_ecckd |
|---|
| 18 | |
|---|
| 19 | use parkind1, only : jprb |
|---|
| 20 | use radiation_gas_constants |
|---|
| 21 | use radiation_ecckd_gas |
|---|
| 22 | use radiation_spectral_definition, only : spectral_definition_type |
|---|
| 23 | |
|---|
| 24 | implicit none |
|---|
| 25 | |
|---|
| 26 | public |
|---|
| 27 | |
|---|
| 28 | !--------------------------------------------------------------------- |
|---|
| 29 | ! This derived type contains all the data needed to describe a |
|---|
| 30 | ! correlated k-distribution gas optics model created using the ecCKD |
|---|
| 31 | ! tool |
|---|
| 32 | type ckd_model_type |
|---|
| 33 | |
|---|
| 34 | ! Gas information |
|---|
| 35 | |
|---|
| 36 | ! Number of gases |
|---|
| 37 | integer :: ngas = 0 |
|---|
| 38 | ! Array of individual-gas data objects |
|---|
| 39 | type(ckd_gas_type), allocatable :: single_gas(:) |
|---|
| 40 | ! Mapping from the "gas codes" in the radiation_gas_constants |
|---|
| 41 | ! module to an index to the single_gas array, where zero means gas |
|---|
| 42 | ! not present (or part of a composite gas) |
|---|
| 43 | integer :: i_gas_mapping(0:NMaxGases) |
|---|
| 44 | |
|---|
| 45 | ! Coordinates of main look-up table for absorption coeffts |
|---|
| 46 | |
|---|
| 47 | ! Number of pressure and temperature points |
|---|
| 48 | integer :: npress = 0 |
|---|
| 49 | integer :: ntemp = 0 |
|---|
| 50 | ! Natural logarithm of first (lowest) pressure (Pa) and increment |
|---|
| 51 | real(jprb) :: log_pressure1, d_log_pressure |
|---|
| 52 | ! First temperature profile (K), dimensioned (npress) |
|---|
| 53 | real(jprb), allocatable :: temperature1(:) |
|---|
| 54 | ! Temperature increment (K) |
|---|
| 55 | real(jprb) :: d_temperature |
|---|
| 56 | |
|---|
| 57 | ! Look-up table for Planck function |
|---|
| 58 | |
|---|
| 59 | ! Number of entries |
|---|
| 60 | integer :: nplanck = 0 |
|---|
| 61 | ! Temperature of first element of look-up table and increment (K) |
|---|
| 62 | real(jprb), allocatable :: temperature1_planck |
|---|
| 63 | real(jprb), allocatable :: d_temperature_planck |
|---|
| 64 | ! Planck function (black body flux into a horizontal plane) in W |
|---|
| 65 | ! m-2, dimensioned (ng,nplanck) |
|---|
| 66 | real(jprb), allocatable :: planck_function(:,:) |
|---|
| 67 | |
|---|
| 68 | ! Normalized solar irradiance in each g point dimensioned (ng) |
|---|
| 69 | real(jprb), allocatable :: norm_solar_irradiance(:) |
|---|
| 70 | |
|---|
| 71 | ! Rayleigh molar scattering coefficient in m2 mol-1 in each g |
|---|
| 72 | ! point |
|---|
| 73 | real(jprb), allocatable :: rayleigh_molar_scat(:) |
|---|
| 74 | |
|---|
| 75 | ! ! Spectral mapping of g points |
|---|
| 76 | |
|---|
| 77 | ! ! Number of wavenumber intervals |
|---|
| 78 | ! integer :: nwav = 0 |
|---|
| 79 | |
|---|
| 80 | ! Number of k terms / g points |
|---|
| 81 | integer :: ng = 0 |
|---|
| 82 | |
|---|
| 83 | ! Spectral definition describing bands and g points |
|---|
| 84 | type(spectral_definition_type) :: spectral_def |
|---|
| 85 | |
|---|
| 86 | ! Shortwave: true, longwave: false |
|---|
| 87 | logical :: is_sw |
|---|
| 88 | |
|---|
| 89 | contains |
|---|
| 90 | |
|---|
| 91 | procedure :: read => read_ckd_model |
|---|
| 92 | procedure :: calc_optical_depth => calc_optical_depth_ckd_model |
|---|
| 93 | procedure :: print => print_ckd_model |
|---|
| 94 | procedure :: calc_planck_function |
|---|
| 95 | procedure :: calc_incoming_sw |
|---|
| 96 | ! procedure :: deallocate => deallocate_ckd_model |
|---|
| 97 | |
|---|
| 98 | end type ckd_model_type |
|---|
| 99 | |
|---|
| 100 | |
|---|
| 101 | contains |
|---|
| 102 | |
|---|
| 103 | !--------------------------------------------------------------------- |
|---|
| 104 | ! Read a complete ecCKD gas optics model from a NetCDF file |
|---|
| 105 | ! "filename" |
|---|
| 106 | subroutine read_ckd_model(this, filename, iverbose) |
|---|
| 107 | |
|---|
| 108 | use easy_netcdf, only : netcdf_file |
|---|
| 109 | !use radiation_io, only : nulerr, radiation_abort |
|---|
| 110 | use yomhook, only : lhook, dr_hook |
|---|
| 111 | |
|---|
| 112 | class(ckd_model_type), intent(inout) :: this |
|---|
| 113 | character(len=*), intent(in) :: filename |
|---|
| 114 | integer, optional, intent(in) :: iverbose |
|---|
| 115 | |
|---|
| 116 | type(netcdf_file) :: file |
|---|
| 117 | |
|---|
| 118 | real(jprb), allocatable :: pressure_lut(:) |
|---|
| 119 | real(jprb), allocatable :: temperature_full(:,:) |
|---|
| 120 | real(jprb), allocatable :: temperature_planck(:) |
|---|
| 121 | |
|---|
| 122 | character(len=512) :: constituent_id |
|---|
| 123 | |
|---|
| 124 | integer :: iverbose_local |
|---|
| 125 | |
|---|
| 126 | ! Loop counters |
|---|
| 127 | integer :: jgas, jjgas |
|---|
| 128 | |
|---|
| 129 | integer :: istart, inext, nchar, i_gas_code |
|---|
| 130 | |
|---|
| 131 | real(jprb) :: hook_handle |
|---|
| 132 | |
|---|
| 133 | if (lhook) call dr_hook('radiation_ecckd:read_ckd_model',0,hook_handle) |
|---|
| 134 | |
|---|
| 135 | if (present(iverbose)) then |
|---|
| 136 | iverbose_local = iverbose |
|---|
| 137 | else |
|---|
| 138 | iverbose_local = 3 |
|---|
| 139 | end if |
|---|
| 140 | |
|---|
| 141 | call file%open(trim(filename), iverbose=iverbose_local) |
|---|
| 142 | |
|---|
| 143 | ! Read temperature and pressure coordinate variables |
|---|
| 144 | call file%get('pressure', pressure_lut) |
|---|
| 145 | this%log_pressure1 = log(pressure_lut(1)) |
|---|
| 146 | this%npress = size(pressure_lut) |
|---|
| 147 | this%d_log_pressure = log(pressure_lut(2)) - this%log_pressure1 |
|---|
| 148 | call file%get('temperature', temperature_full) |
|---|
| 149 | allocate(this%temperature1(this%npress)); |
|---|
| 150 | this%temperature1 = temperature_full(:,1) |
|---|
| 151 | this%d_temperature = temperature_full(1,2)-temperature_full(1,1) |
|---|
| 152 | this%ntemp = size(temperature_full,2) |
|---|
| 153 | deallocate(temperature_full) |
|---|
| 154 | |
|---|
| 155 | ! Read Planck function, or solar irradiance and Rayleigh |
|---|
| 156 | ! scattering coefficient |
|---|
| 157 | if (file%exists('solar_irradiance')) then |
|---|
| 158 | this%is_sw = .true. |
|---|
| 159 | call file%get('solar_irradiance', this%norm_solar_irradiance) |
|---|
| 160 | this%norm_solar_irradiance = this%norm_solar_irradiance & |
|---|
| 161 | & / sum(this%norm_solar_irradiance) |
|---|
| 162 | call file%get('rayleigh_molar_scattering_coeff', & |
|---|
| 163 | & this%rayleigh_molar_scat) |
|---|
| 164 | else |
|---|
| 165 | this%is_sw = .false. |
|---|
| 166 | call file%get('temperature_planck', temperature_planck) |
|---|
| 167 | this%nplanck = size(temperature_planck) |
|---|
| 168 | this%temperature1_planck = temperature_planck(1) |
|---|
| 169 | this%d_temperature_planck = temperature_planck(2) - temperature_planck(1) |
|---|
| 170 | deallocate(temperature_planck) |
|---|
| 171 | call file%get('planck_function', this%planck_function) |
|---|
| 172 | end if |
|---|
| 173 | |
|---|
| 174 | ! Read the spectral definition information into a separate |
|---|
| 175 | ! derived type |
|---|
| 176 | call this%spectral_def%read(file); |
|---|
| 177 | this%ng = this%spectral_def%ng |
|---|
| 178 | |
|---|
| 179 | ! Read gases |
|---|
| 180 | call file%get('n_gases', this%ngas) |
|---|
| 181 | allocate(this%single_gas(this%ngas)) |
|---|
| 182 | call file%get_global_attribute('constituent_id', constituent_id) |
|---|
| 183 | nchar = len(trim(constituent_id)) |
|---|
| 184 | istart = 1 |
|---|
| 185 | this%i_gas_mapping = 0 |
|---|
| 186 | do jgas = 1, this%ngas |
|---|
| 187 | if (jgas < this%ngas) then |
|---|
| 188 | inext = istart + scan(constituent_id(istart:nchar), ' ') |
|---|
| 189 | else |
|---|
| 190 | inext = nchar+2 |
|---|
| 191 | end if |
|---|
| 192 | ! Find gas code |
|---|
| 193 | i_gas_code = 0 |
|---|
| 194 | do jjgas = 1, NMaxGases |
|---|
| 195 | if (constituent_id(istart:inext-2) == trim(GasLowerCaseName(jjgas))) then |
|---|
| 196 | i_gas_code = jjgas |
|---|
| 197 | exit |
|---|
| 198 | end if |
|---|
| 199 | end do |
|---|
| 200 | ! if (i_gas_code == 0) then |
|---|
| 201 | ! write(nulerr,'(a,a,a)') '*** Error: Gas "', constituent_id(istart:inext-2), & |
|---|
| 202 | ! & '" not understood' |
|---|
| 203 | ! call radiation_abort('Radiation configuration error') |
|---|
| 204 | ! end if |
|---|
| 205 | this%i_gas_mapping(i_gas_code) = jgas; |
|---|
| 206 | call this%single_gas(jgas)%read(file, constituent_id(istart:inext-2), i_gas_code) |
|---|
| 207 | istart = inext |
|---|
| 208 | end do |
|---|
| 209 | |
|---|
| 210 | if (lhook) call dr_hook('radiation_ecckd:read_ckd_model',1,hook_handle) |
|---|
| 211 | |
|---|
| 212 | end subroutine read_ckd_model |
|---|
| 213 | |
|---|
| 214 | !--------------------------------------------------------------------- |
|---|
| 215 | ! Print a description of the correlated k-distribution model to the |
|---|
| 216 | ! "nulout" unit |
|---|
| 217 | subroutine print_ckd_model(this) |
|---|
| 218 | |
|---|
| 219 | use radiation_io, only : nulout |
|---|
| 220 | use radiation_gas_constants |
|---|
| 221 | |
|---|
| 222 | class(ckd_model_type), intent(in) :: this |
|---|
| 223 | |
|---|
| 224 | integer :: jgas |
|---|
| 225 | |
|---|
| 226 | if (this%is_sw) then |
|---|
| 227 | write(nulout,'(a)',advance='no') 'ecCKD shortwave gas optics model: ' |
|---|
| 228 | else |
|---|
| 229 | write(nulout,'(a)',advance='no') 'ecCKD longwave gas optics model: ' |
|---|
| 230 | end if |
|---|
| 231 | |
|---|
| 232 | write(nulout,'(i0,a,i0,a,i0,a,i0,a)') & |
|---|
| 233 | & nint(this%spectral_def%wavenumber1(1)), '-', & |
|---|
| 234 | & nint(this%spectral_def%wavenumber2(size(this%spectral_def%wavenumber2))), & |
|---|
| 235 | & ' cm-1, ', this%ng, ' g-points in ', this%spectral_def%nband, ' bands' |
|---|
| 236 | write(nulout,'(a,i0,a,i0,a,i0,a)') ' Look-up table sizes: ', this%npress, ' pressures, ', & |
|---|
| 237 | & this%ntemp, ' temperatures, ', this%nplanck, ' planck-function entries' |
|---|
| 238 | write(nulout, '(a)') ' Gases:' |
|---|
| 239 | do jgas = 1,this%ngas |
|---|
| 240 | if (this%single_gas(jgas)%i_gas_code > 0) then |
|---|
| 241 | write(nulout, '(a,a,a)', advance='no') ' ', & |
|---|
| 242 | & trim(GasName(this%single_gas(jgas)%i_gas_code)), ': ' |
|---|
| 243 | else |
|---|
| 244 | write(nulout, '(a)', advance='no') ' Composite of well-mixed background gases: ' |
|---|
| 245 | end if |
|---|
| 246 | select case (this%single_gas(jgas)%i_conc_dependence) |
|---|
| 247 | case (IConcDependenceNone) |
|---|
| 248 | write(nulout, '(a)') 'no concentration dependence' |
|---|
| 249 | case (IConcDependenceLinear) |
|---|
| 250 | write(nulout, '(a)') 'linear concentration dependence' |
|---|
| 251 | case (IConcDependenceRelativeLinear) |
|---|
| 252 | write(nulout, '(a,e14.6)') 'linear concentration dependence relative to a mole fraction of ', & |
|---|
| 253 | & this%single_gas(jgas)%reference_mole_frac |
|---|
| 254 | case (IConcDependenceLUT) |
|---|
| 255 | write(nulout, '(a,i0,a,e14.6,a,e13.6)') 'look-up table with ', this%single_gas(jgas)%n_mole_frac, & |
|---|
| 256 | & ' log-spaced mole fractions in range ', exp(this%single_gas(jgas)%log_mole_frac1), & |
|---|
| 257 | & ' to ', exp(this%single_gas(jgas)%log_mole_frac1 & |
|---|
| 258 | & + this%single_gas(jgas)%n_mole_frac*this%single_gas(jgas)%d_log_mole_frac) |
|---|
| 259 | end select |
|---|
| 260 | end do |
|---|
| 261 | |
|---|
| 262 | end subroutine print_ckd_model |
|---|
| 263 | |
|---|
| 264 | |
|---|
| 265 | !--------------------------------------------------------------------- |
|---|
| 266 | ! Compute layerwise optical depth for each g point for ncol columns |
|---|
| 267 | ! at nlev layers |
|---|
| 268 | subroutine calc_optical_depth_ckd_model(this, ncol, nlev, istartcol, iendcol, nmaxgas, & |
|---|
| 269 | & pressure_hl, temperature_fl, mole_fraction_fl, & |
|---|
| 270 | & optical_depth_fl, rayleigh_od_fl) |
|---|
| 271 | |
|---|
| 272 | use yomhook, only : lhook, dr_hook |
|---|
| 273 | use radiation_constants, only : AccelDueToGravity |
|---|
| 274 | |
|---|
| 275 | ! Input variables |
|---|
| 276 | |
|---|
| 277 | class(ckd_model_type), intent(in), target :: this |
|---|
| 278 | ! Number of columns, levels and input gases |
|---|
| 279 | integer, intent(in) :: ncol, nlev, nmaxgas, istartcol, iendcol |
|---|
| 280 | ! Pressure at half levels (Pa), dimensioned (ncol,nlev+1) |
|---|
| 281 | real(jprb), intent(in) :: pressure_hl(ncol,nlev+1) |
|---|
| 282 | ! Temperature at full levels (K), dimensioned (ncol,nlev) |
|---|
| 283 | real(jprb), intent(in) :: temperature_fl(istartcol:iendcol,nlev) |
|---|
| 284 | ! Gas mole fractions at full levels (mol mol-1), dimensioned (ncol,nlev,nmaxgas) |
|---|
| 285 | real(jprb), intent(in) :: mole_fraction_fl(ncol,nlev,nmaxgas) |
|---|
| 286 | |
|---|
| 287 | ! Output variables |
|---|
| 288 | |
|---|
| 289 | ! Layer absorption optical depth for each g point |
|---|
| 290 | real(jprb), intent(out) :: optical_depth_fl(this%ng,nlev,istartcol:iendcol) |
|---|
| 291 | ! In the shortwave only, the Rayleigh scattering optical depth |
|---|
| 292 | real(jprb), optional, intent(out) :: rayleigh_od_fl(this%ng,nlev,istartcol:iendcol) |
|---|
| 293 | |
|---|
| 294 | ! Local variables |
|---|
| 295 | |
|---|
| 296 | real(jprb), pointer :: molar_abs(:,:,:), molar_abs_conc(:,:,:,:) |
|---|
| 297 | |
|---|
| 298 | ! Natural logarithm of pressure at full levels |
|---|
| 299 | real(jprb) :: log_pressure_fl(nlev) |
|---|
| 300 | |
|---|
| 301 | ! Optical depth of single gas at one point in space versus |
|---|
| 302 | ! spectral interval |
|---|
| 303 | !real(jprb) :: od_single_gas(this%ng) |
|---|
| 304 | |
|---|
| 305 | real(jprb) :: multiplier(nlev), simple_multiplier(nlev), global_multiplier, temperature1 |
|---|
| 306 | |
|---|
| 307 | ! Indices and weights in temperature, pressure and concentration interpolation |
|---|
| 308 | real(jprb) :: pindex1, tindex1, cindex1 |
|---|
| 309 | real(jprb) :: pw1(nlev), pw2(nlev), tw1(nlev), tw2(nlev), cw1(nlev), cw2(nlev) |
|---|
| 310 | integer :: ip1(nlev), it1(nlev), ic1(nlev) |
|---|
| 311 | |
|---|
| 312 | ! Natural logarithm of mole fraction at one point |
|---|
| 313 | real(jprb) :: log_conc |
|---|
| 314 | |
|---|
| 315 | ! Minimum mole fraction in look-up-table |
|---|
| 316 | real(jprb) :: mole_frac1 |
|---|
| 317 | |
|---|
| 318 | integer :: jcol, jlev, jgas, igascode |
|---|
| 319 | |
|---|
| 320 | real(jprb) :: hook_handle |
|---|
| 321 | |
|---|
| 322 | if (lhook) call dr_hook('radiation_ecckd:calc_optical_depth',0,hook_handle) |
|---|
| 323 | |
|---|
| 324 | global_multiplier = 1.0_jprb / (AccelDueToGravity * 0.001_jprb * AirMolarMass) |
|---|
| 325 | |
|---|
| 326 | do jcol = istartcol,iendcol |
|---|
| 327 | |
|---|
| 328 | log_pressure_fl = log(0.5_jprb * (pressure_hl(jcol,1:nlev)+pressure_hl(jcol,2:nlev+1))) |
|---|
| 329 | |
|---|
| 330 | do jlev = 1,nlev |
|---|
| 331 | ! Find interpolation points in pressure |
|---|
| 332 | pindex1 = (log_pressure_fl(jlev)-this%log_pressure1) & |
|---|
| 333 | & / this%d_log_pressure |
|---|
| 334 | pindex1 = 1.0_jprb + max(0.0_jprb, min(pindex1, this%npress-1.0001_jprb)) |
|---|
| 335 | ip1(jlev) = int(pindex1) |
|---|
| 336 | pw2(jlev) = pindex1 - ip1(jlev) |
|---|
| 337 | pw1(jlev) = 1.0_jprb - pw2(jlev) |
|---|
| 338 | |
|---|
| 339 | ! Find interpolation points in temperature |
|---|
| 340 | temperature1 = pw1(jlev)*this%temperature1(ip1(jlev)) & |
|---|
| 341 | & + pw2(jlev)*this%temperature1(ip1(jlev)+1) |
|---|
| 342 | tindex1 = (temperature_fl(jcol,jlev) - temperature1) & |
|---|
| 343 | & / this%d_temperature |
|---|
| 344 | tindex1 = 1.0_jprb + max(0.0_jprb, min(tindex1, this%ntemp-1.0001_jprb)) |
|---|
| 345 | it1(jlev) = int(tindex1) |
|---|
| 346 | tw2(jlev) = tindex1 - it1(jlev) |
|---|
| 347 | tw1(jlev) = 1.0_jprb - tw2(jlev) |
|---|
| 348 | |
|---|
| 349 | ! Concentration multiplier |
|---|
| 350 | simple_multiplier(jlev) = global_multiplier & |
|---|
| 351 | & * (pressure_hl(jcol,jlev+1) - pressure_hl(jcol,jlev)) |
|---|
| 352 | end do |
|---|
| 353 | |
|---|
| 354 | optical_depth_fl(:,:,jcol) = 0.0_jprb |
|---|
| 355 | |
|---|
| 356 | do jgas = 1,this%ngas |
|---|
| 357 | |
|---|
| 358 | associate (single_gas => this%single_gas(jgas)) |
|---|
| 359 | igascode = this%single_gas(jgas)%i_gas_code |
|---|
| 360 | |
|---|
| 361 | select case (single_gas%i_conc_dependence) |
|---|
| 362 | |
|---|
| 363 | case (IConcDependenceLinear) |
|---|
| 364 | molar_abs => this%single_gas(jgas)%molar_abs |
|---|
| 365 | multiplier = simple_multiplier * mole_fraction_fl(jcol,:,igascode) |
|---|
| 366 | |
|---|
| 367 | do jlev = 1,nlev |
|---|
| 368 | optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) & |
|---|
| 369 | & + (multiplier(jlev)*tw1(jlev)) * (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)) & |
|---|
| 370 | & +pw2(jlev) * molar_abs(:,ip1(jlev)+1,it1(jlev))) & |
|---|
| 371 | & + (multiplier(jlev)*tw2(jlev)) * (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)+1) & |
|---|
| 372 | & +pw2(jlev) * molar_abs(:,ip1(jlev)+1,it1(jlev)+1)) |
|---|
| 373 | end do |
|---|
| 374 | |
|---|
| 375 | case (IConcDependenceRelativeLinear) |
|---|
| 376 | molar_abs => this%single_gas(jgas)%molar_abs |
|---|
| 377 | multiplier = simple_multiplier * (mole_fraction_fl(jcol,:,igascode) & |
|---|
| 378 | & - single_gas%reference_mole_frac) |
|---|
| 379 | do jlev = 1,nlev |
|---|
| 380 | optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) & |
|---|
| 381 | & + (multiplier(jlev)*tw1(jlev)) * (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)) & |
|---|
| 382 | & +pw2(jlev) * molar_abs(:,ip1(jlev)+1,it1(jlev))) & |
|---|
| 383 | & + (multiplier(jlev)*tw2(jlev)) * (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)+1) & |
|---|
| 384 | & +pw2(jlev) * molar_abs(:,ip1(jlev)+1,it1(jlev)+1)) |
|---|
| 385 | end do |
|---|
| 386 | |
|---|
| 387 | case (IConcDependenceNone) |
|---|
| 388 | ! Composite gases |
|---|
| 389 | molar_abs => this%single_gas(jgas)%molar_abs |
|---|
| 390 | do jlev = 1,nlev |
|---|
| 391 | optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) & |
|---|
| 392 | & + (simple_multiplier(jlev)*tw1(jlev)) * (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)) & |
|---|
| 393 | & +pw2(jlev) * molar_abs(:,ip1(jlev)+1,it1(jlev))) & |
|---|
| 394 | & + (simple_multiplier(jlev)*tw2(jlev)) * (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)+1) & |
|---|
| 395 | & +pw2(jlev) * molar_abs(:,ip1(jlev)+1,it1(jlev)+1)) |
|---|
| 396 | end do |
|---|
| 397 | |
|---|
| 398 | case (IConcDependenceLUT) |
|---|
| 399 | ! Logarithmic interpolation in concentration space |
|---|
| 400 | molar_abs_conc => this%single_gas(jgas)%molar_abs_conc |
|---|
| 401 | mole_frac1 = exp(single_gas%log_mole_frac1) |
|---|
| 402 | do jlev = 1,nlev |
|---|
| 403 | ! Take care of mole_fraction == 0 |
|---|
| 404 | log_conc = log(max(mole_fraction_fl(jcol,jlev,igascode), mole_frac1)) |
|---|
| 405 | cindex1 = (log_conc - single_gas%log_mole_frac1) / single_gas%d_log_mole_frac |
|---|
| 406 | cindex1 = 1.0_jprb + max(0.0_jprb, min(cindex1, single_gas%n_mole_frac-1.0001_jprb)) |
|---|
| 407 | ic1(jlev) = int(cindex1) |
|---|
| 408 | cw2(jlev) = cindex1 - ic1(jlev) |
|---|
| 409 | cw1(jlev) = 1.0_jprb - cw2(jlev) |
|---|
| 410 | end do |
|---|
| 411 | ! od_single_gas = cw1 * (tw1 * (pw1 * molar_abs_conc(:,ip1,it1,ic1) & |
|---|
| 412 | ! & +pw2 * molar_abs_conc(:,ip1+1,it1,ic1)) & |
|---|
| 413 | ! & +tw2 * (pw1 * molar_abs_conc(:,ip1,it1+1,ic1) & |
|---|
| 414 | ! & +pw2 * molar_abs_conc(:,ip1+1,it1+1,ic1))) & |
|---|
| 415 | ! & +cw2 * (tw1 * (pw1 * molar_abs_conc(:,ip1,it1,ic1+1) & |
|---|
| 416 | ! & +pw2 * molar_abs_conc(:,ip1+1,it1,ic1+1)) & |
|---|
| 417 | ! & +tw2 * (pw1 * molar_abs_conc(:,ip1,it1+1,ic1+1) & |
|---|
| 418 | ! & +pw2 * molar_abs_conc(:,ip1+1,it1+1,ic1+1))) |
|---|
| 419 | do jlev = 1,nlev |
|---|
| 420 | optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) & |
|---|
| 421 | & + (simple_multiplier(jlev) * mole_fraction_fl(jcol,jlev,igascode)) * ( & |
|---|
| 422 | & (cw1(jlev) * tw1(jlev) * pw1(jlev)) * molar_abs_conc(:,ip1(jlev),it1(jlev),ic1(jlev)) & |
|---|
| 423 | & +(cw1(jlev) * tw1(jlev) * pw2(jlev)) * molar_abs_conc(:,ip1(jlev)+1,it1(jlev),ic1(jlev)) & |
|---|
| 424 | & +(cw1(jlev) * tw2(jlev) * pw1(jlev)) * molar_abs_conc(:,ip1(jlev),it1(jlev)+1,ic1(jlev)) & |
|---|
| 425 | & +(cw1(jlev) * tw2(jlev) * pw2(jlev)) * molar_abs_conc(:,ip1(jlev)+1,it1(jlev)+1,ic1(jlev)) & |
|---|
| 426 | & +(cw2(jlev) * tw1(jlev) * pw1(jlev)) * molar_abs_conc(:,ip1(jlev),it1(jlev),ic1(jlev)+1) & |
|---|
| 427 | & +(cw2(jlev) * tw1(jlev) * pw2(jlev)) * molar_abs_conc(:,ip1(jlev)+1,it1(jlev),ic1(jlev)+1) & |
|---|
| 428 | & +(cw2(jlev) * tw2(jlev) * pw1(jlev)) * molar_abs_conc(:,ip1(jlev),it1(jlev)+1,ic1(jlev)+1) & |
|---|
| 429 | & +(cw2(jlev) * tw2(jlev) * pw2(jlev)) * molar_abs_conc(:,ip1(jlev)+1,it1(jlev)+1,ic1(jlev)+1)) |
|---|
| 430 | end do |
|---|
| 431 | end select |
|---|
| 432 | |
|---|
| 433 | end associate |
|---|
| 434 | |
|---|
| 435 | end do |
|---|
| 436 | |
|---|
| 437 | ! Ensure the optical depth is not negative |
|---|
| 438 | optical_depth_fl(:,:,jcol) = max(0.0_jprb, optical_depth_fl(:,:,jcol)) |
|---|
| 439 | |
|---|
| 440 | ! Rayleigh scattering |
|---|
| 441 | if (this%is_sw .and. present(rayleigh_od_fl)) then |
|---|
| 442 | do jlev = 1,nlev |
|---|
| 443 | rayleigh_od_fl(:,jlev,jcol) = global_multiplier & |
|---|
| 444 | & * (pressure_hl(jcol,jlev+1) - pressure_hl(jcol,jlev)) * this%rayleigh_molar_scat |
|---|
| 445 | end do |
|---|
| 446 | end if |
|---|
| 447 | |
|---|
| 448 | end do |
|---|
| 449 | |
|---|
| 450 | if (lhook) call dr_hook('radiation_ecckd:calc_optical_depth',1,hook_handle) |
|---|
| 451 | |
|---|
| 452 | end subroutine calc_optical_depth_ckd_model |
|---|
| 453 | |
|---|
| 454 | !--------------------------------------------------------------------- |
|---|
| 455 | ! Calculate the Planck function integrated across each of the g |
|---|
| 456 | ! points of this correlated k-distribution model, for a given |
|---|
| 457 | ! temperature, where Planck function is defined as the flux emitted |
|---|
| 458 | ! by a black body (rather than radiance) |
|---|
| 459 | subroutine calc_planck_function(this, nt, temperature, planck) |
|---|
| 460 | |
|---|
| 461 | class(ckd_model_type), intent(in) :: this |
|---|
| 462 | integer, intent(in) :: nt |
|---|
| 463 | real(jprb), intent(in) :: temperature(:) ! K |
|---|
| 464 | real(jprb), intent(out) :: planck(this%ng,nt) ! W m-2 |
|---|
| 465 | |
|---|
| 466 | real(jprb) :: tindex1, tw1, tw2 |
|---|
| 467 | integer :: it1, jt |
|---|
| 468 | |
|---|
| 469 | do jt = 1,nt |
|---|
| 470 | tindex1 = (temperature(jt) - this%temperature1_planck) & |
|---|
| 471 | & * (1.0_jprb / this%d_temperature_planck) |
|---|
| 472 | if (tindex1 >= 0) then |
|---|
| 473 | ! Normal interpolation, and extrapolation for high temperatures |
|---|
| 474 | tindex1 = 1.0_jprb + tindex1 |
|---|
| 475 | it1 = min(int(tindex1), this%nplanck-1) |
|---|
| 476 | tw2 = tindex1 - it1 |
|---|
| 477 | tw1 = 1.0_jprb - tw2 |
|---|
| 478 | planck(:,jt) = tw1 * this%planck_function(:,it1) & |
|---|
| 479 | & + tw2 * this%planck_function(:,it1+1) |
|---|
| 480 | else |
|---|
| 481 | ! Interpolate linearly to zero |
|---|
| 482 | planck(:,jt) = this%planck_function(:,1) & |
|---|
| 483 | & * (temperature(jt)/this%temperature1_planck) |
|---|
| 484 | end if |
|---|
| 485 | end do |
|---|
| 486 | |
|---|
| 487 | end subroutine calc_planck_function |
|---|
| 488 | |
|---|
| 489 | |
|---|
| 490 | !--------------------------------------------------------------------- |
|---|
| 491 | ! Return the spectral solar irradiance integrated over each g point |
|---|
| 492 | ! of a solar correlated k-distribution model, given the |
|---|
| 493 | ! total_solar_irradiance |
|---|
| 494 | subroutine calc_incoming_sw(this, total_solar_irradiance, & |
|---|
| 495 | & spectral_solar_irradiance) |
|---|
| 496 | |
|---|
| 497 | class(ckd_model_type), intent(in) :: this |
|---|
| 498 | real(jprb), intent(in) :: total_solar_irradiance ! W m-2 |
|---|
| 499 | real(jprb), intent(inout) :: spectral_solar_irradiance(:,:) ! W m-2 |
|---|
| 500 | |
|---|
| 501 | spectral_solar_irradiance & |
|---|
| 502 | & = spread(total_solar_irradiance * this%norm_solar_irradiance, & |
|---|
| 503 | & 2, size(spectral_solar_irradiance,2)) |
|---|
| 504 | |
|---|
| 505 | end subroutine calc_incoming_sw |
|---|
| 506 | |
|---|
| 507 | end module radiation_ecckd |
|---|
