| 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 | ! Normalized amplitude of variations in the solar irradiance |
|---|
| 72 | ! through the solar cycle in each g point, dimensioned (ng). |
|---|
| 73 | ! Since the user always provides the solar irradiance SI |
|---|
| 74 | ! integrated across the spectrum, this variable must sum to zero: |
|---|
| 75 | ! this ensures that the solar irradiance in each g-point is |
|---|
| 76 | ! SSI=SI*(norm_solar_irradiance + |
|---|
| 77 | ! A*norm_amplitude_solar_irradiance) for any A, where A denotes |
|---|
| 78 | ! the amplitude of deviations from the mean solar spectrum, |
|---|
| 79 | ! typically between -1.0 and 1.0 and provided by |
|---|
| 80 | ! single_level%solar_spectral_multiplier. |
|---|
| 81 | real(jprb), allocatable :: norm_amplitude_solar_irradiance(:) |
|---|
| 82 | |
|---|
| 83 | ! Rayleigh molar scattering coefficient in m2 mol-1 in each g |
|---|
| 84 | ! point |
|---|
| 85 | real(jprb), allocatable :: rayleigh_molar_scat(:) |
|---|
| 86 | |
|---|
| 87 | ! ! Spectral mapping of g points |
|---|
| 88 | |
|---|
| 89 | ! ! Number of wavenumber intervals |
|---|
| 90 | ! integer :: nwav = 0 |
|---|
| 91 | |
|---|
| 92 | ! Number of k terms / g points |
|---|
| 93 | integer :: ng = 0 |
|---|
| 94 | |
|---|
| 95 | ! Spectral definition describing bands and g points |
|---|
| 96 | type(spectral_definition_type) :: spectral_def |
|---|
| 97 | |
|---|
| 98 | ! Shortwave: true, longwave: false |
|---|
| 99 | logical :: is_sw |
|---|
| 100 | |
|---|
| 101 | contains |
|---|
| 102 | |
|---|
| 103 | procedure :: read => read_ckd_model |
|---|
| 104 | procedure :: read_spectral_solar_cycle |
|---|
| 105 | ! Vectorized version of the optical depth look-up performs better on |
|---|
| 106 | ! NEC, but slower on x86 |
|---|
| 107 | #ifdef __SX__ |
|---|
| 108 | procedure :: calc_optical_depth => calc_optical_depth_ckd_model_vec |
|---|
| 109 | #else |
|---|
| 110 | procedure :: calc_optical_depth => calc_optical_depth_ckd_model |
|---|
| 111 | #endif |
|---|
| 112 | procedure :: print => print_ckd_model |
|---|
| 113 | procedure :: calc_planck_function |
|---|
| 114 | procedure :: calc_incoming_sw |
|---|
| 115 | ! procedure :: deallocate => deallocate_ckd_model |
|---|
| 116 | |
|---|
| 117 | end type ckd_model_type |
|---|
| 118 | |
|---|
| 119 | |
|---|
| 120 | contains |
|---|
| 121 | |
|---|
| 122 | !--------------------------------------------------------------------- |
|---|
| 123 | ! Read a complete ecCKD gas optics model from a NetCDF file |
|---|
| 124 | ! "filename" |
|---|
| 125 | subroutine read_ckd_model(this, filename, iverbose) |
|---|
| 126 | |
|---|
| 127 | use easy_netcdf, only : netcdf_file |
|---|
| 128 | !use radiation_io, only : nulerr, radiation_abort |
|---|
| 129 | use yomhook, only : lhook, dr_hook, jphook |
|---|
| 130 | |
|---|
| 131 | class(ckd_model_type), intent(inout) :: this |
|---|
| 132 | character(len=*), intent(in) :: filename |
|---|
| 133 | integer, optional, intent(in) :: iverbose |
|---|
| 134 | |
|---|
| 135 | type(netcdf_file) :: file |
|---|
| 136 | |
|---|
| 137 | real(jprb), allocatable :: pressure_lut(:) |
|---|
| 138 | real(jprb), allocatable :: temperature_full(:,:) |
|---|
| 139 | real(jprb), allocatable :: temperature_planck(:) |
|---|
| 140 | |
|---|
| 141 | character(len=512) :: constituent_id |
|---|
| 142 | |
|---|
| 143 | integer :: iverbose_local |
|---|
| 144 | |
|---|
| 145 | ! Loop counters |
|---|
| 146 | integer :: jgas, jjgas |
|---|
| 147 | |
|---|
| 148 | integer :: istart, inext, nchar, i_gas_code |
|---|
| 149 | |
|---|
| 150 | real(jphook) :: hook_handle |
|---|
| 151 | |
|---|
| 152 | if (lhook) call dr_hook('radiation_ecckd:read_ckd_model',0,hook_handle) |
|---|
| 153 | |
|---|
| 154 | if (present(iverbose)) then |
|---|
| 155 | iverbose_local = iverbose |
|---|
| 156 | else |
|---|
| 157 | iverbose_local = 3 |
|---|
| 158 | end if |
|---|
| 159 | |
|---|
| 160 | call file%open(trim(filename), iverbose=iverbose_local) |
|---|
| 161 | |
|---|
| 162 | ! Read temperature and pressure coordinate variables |
|---|
| 163 | call file%get('pressure', pressure_lut) |
|---|
| 164 | this%log_pressure1 = log(pressure_lut(1)) |
|---|
| 165 | this%npress = size(pressure_lut) |
|---|
| 166 | this%d_log_pressure = log(pressure_lut(2)) - this%log_pressure1 |
|---|
| 167 | call file%get('temperature', temperature_full) |
|---|
| 168 | ! AI oct 2023 ajout pour le double appel de ecrad |
|---|
| 169 | if (allocated(this%temperature1)) deallocate(this%temperature1) |
|---|
| 170 | allocate(this%temperature1(this%npress)); |
|---|
| 171 | this%temperature1 = temperature_full(:,1) |
|---|
| 172 | this%d_temperature = temperature_full(1,2)-temperature_full(1,1) |
|---|
| 173 | this%ntemp = size(temperature_full,2) |
|---|
| 174 | deallocate(temperature_full) |
|---|
| 175 | |
|---|
| 176 | ! Read Planck function, or solar irradiance and Rayleigh |
|---|
| 177 | ! scattering coefficient |
|---|
| 178 | if (file%exists('solar_irradiance')) then |
|---|
| 179 | this%is_sw = .true. |
|---|
| 180 | call file%get('solar_irradiance', this%norm_solar_irradiance) |
|---|
| 181 | this%norm_solar_irradiance = this%norm_solar_irradiance & |
|---|
| 182 | & / sum(this%norm_solar_irradiance) |
|---|
| 183 | call file%get('rayleigh_molar_scattering_coeff', & |
|---|
| 184 | & this%rayleigh_molar_scat) |
|---|
| 185 | else |
|---|
| 186 | this%is_sw = .false. |
|---|
| 187 | call file%get('temperature_planck', temperature_planck) |
|---|
| 188 | this%nplanck = size(temperature_planck) |
|---|
| 189 | this%temperature1_planck = temperature_planck(1) |
|---|
| 190 | this%d_temperature_planck = temperature_planck(2) - temperature_planck(1) |
|---|
| 191 | deallocate(temperature_planck) |
|---|
| 192 | call file%get('planck_function', this%planck_function) |
|---|
| 193 | end if |
|---|
| 194 | |
|---|
| 195 | ! Read the spectral definition information into a separate |
|---|
| 196 | ! derived type |
|---|
| 197 | call this%spectral_def%read(file); |
|---|
| 198 | this%ng = this%spectral_def%ng |
|---|
| 199 | |
|---|
| 200 | ! Read gases |
|---|
| 201 | call file%get('n_gases', this%ngas) |
|---|
| 202 | if (allocated(this%single_gas)) deallocate(this%single_gas) |
|---|
| 203 | allocate(this%single_gas(this%ngas)) |
|---|
| 204 | call file%get_global_attribute('constituent_id', constituent_id) |
|---|
| 205 | nchar = len(trim(constituent_id)) |
|---|
| 206 | istart = 1 |
|---|
| 207 | this%i_gas_mapping = 0 |
|---|
| 208 | do jgas = 1, this%ngas |
|---|
| 209 | if (jgas < this%ngas) then |
|---|
| 210 | inext = istart + scan(constituent_id(istart:nchar), ' ') |
|---|
| 211 | else |
|---|
| 212 | inext = nchar+2 |
|---|
| 213 | end if |
|---|
| 214 | ! Find gas code |
|---|
| 215 | i_gas_code = 0 |
|---|
| 216 | do jjgas = 1, NMaxGases |
|---|
| 217 | if (constituent_id(istart:inext-2) == trim(GasLowerCaseName(jjgas))) then |
|---|
| 218 | i_gas_code = jjgas |
|---|
| 219 | exit |
|---|
| 220 | end if |
|---|
| 221 | end do |
|---|
| 222 | ! if (i_gas_code == 0) then |
|---|
| 223 | ! write(nulerr,'(a,a,a)') '*** Error: Gas "', constituent_id(istart:inext-2), & |
|---|
| 224 | ! & '" not understood' |
|---|
| 225 | ! call radiation_abort('Radiation configuration error') |
|---|
| 226 | ! end if |
|---|
| 227 | this%i_gas_mapping(i_gas_code) = jgas; |
|---|
| 228 | call this%single_gas(jgas)%read(file, constituent_id(istart:inext-2), i_gas_code) |
|---|
| 229 | istart = inext |
|---|
| 230 | end do |
|---|
| 231 | |
|---|
| 232 | call file%close() |
|---|
| 233 | |
|---|
| 234 | if (lhook) call dr_hook('radiation_ecckd:read_ckd_model',1,hook_handle) |
|---|
| 235 | |
|---|
| 236 | end subroutine read_ckd_model |
|---|
| 237 | |
|---|
| 238 | !--------------------------------------------------------------------- |
|---|
| 239 | ! Print a description of the correlated k-distribution model to the |
|---|
| 240 | ! "nulout" unit |
|---|
| 241 | subroutine print_ckd_model(this) |
|---|
| 242 | |
|---|
| 243 | use radiation_io, only : nulout |
|---|
| 244 | use radiation_gas_constants |
|---|
| 245 | |
|---|
| 246 | class(ckd_model_type), intent(in) :: this |
|---|
| 247 | |
|---|
| 248 | integer :: jgas |
|---|
| 249 | |
|---|
| 250 | if (this%is_sw) then |
|---|
| 251 | write(nulout,'(a)',advance='no') 'ecCKD shortwave gas optics model: ' |
|---|
| 252 | else |
|---|
| 253 | write(nulout,'(a)',advance='no') 'ecCKD longwave gas optics model: ' |
|---|
| 254 | end if |
|---|
| 255 | |
|---|
| 256 | write(nulout,'(i0,a,i0,a,i0,a,i0,a)') & |
|---|
| 257 | & nint(this%spectral_def%wavenumber1(1)), '-', & |
|---|
| 258 | & nint(this%spectral_def%wavenumber2(size(this%spectral_def%wavenumber2))), & |
|---|
| 259 | & ' cm-1, ', this%ng, ' g-points in ', this%spectral_def%nband, ' bands' |
|---|
| 260 | write(nulout,'(a,i0,a,i0,a,i0,a)') ' Look-up table sizes: ', this%npress, ' pressures, ', & |
|---|
| 261 | & this%ntemp, ' temperatures, ', this%nplanck, ' planck-function entries' |
|---|
| 262 | write(nulout, '(a)') ' Gases:' |
|---|
| 263 | do jgas = 1,this%ngas |
|---|
| 264 | if (this%single_gas(jgas)%i_gas_code > 0) then |
|---|
| 265 | write(nulout, '(a,a,a)', advance='no') ' ', & |
|---|
| 266 | & trim(GasName(this%single_gas(jgas)%i_gas_code)), ': ' |
|---|
| 267 | else |
|---|
| 268 | write(nulout, '(a)', advance='no') ' Composite of well-mixed background gases: ' |
|---|
| 269 | end if |
|---|
| 270 | select case (this%single_gas(jgas)%i_conc_dependence) |
|---|
| 271 | case (IConcDependenceNone) |
|---|
| 272 | write(nulout, '(a)') 'no concentration dependence' |
|---|
| 273 | case (IConcDependenceLinear) |
|---|
| 274 | write(nulout, '(a)') 'linear concentration dependence' |
|---|
| 275 | case (IConcDependenceRelativeLinear) |
|---|
| 276 | write(nulout, '(a,e14.6)') 'linear concentration dependence relative to a mole fraction of ', & |
|---|
| 277 | & this%single_gas(jgas)%reference_mole_frac |
|---|
| 278 | case (IConcDependenceLUT) |
|---|
| 279 | write(nulout, '(a,i0,a,e14.6,a,e13.6)') 'look-up table with ', this%single_gas(jgas)%n_mole_frac, & |
|---|
| 280 | & ' log-spaced mole fractions in range ', exp(this%single_gas(jgas)%log_mole_frac1), & |
|---|
| 281 | & ' to ', exp(this%single_gas(jgas)%log_mole_frac1 & |
|---|
| 282 | & + this%single_gas(jgas)%n_mole_frac*this%single_gas(jgas)%d_log_mole_frac) |
|---|
| 283 | end select |
|---|
| 284 | end do |
|---|
| 285 | |
|---|
| 286 | end subroutine print_ckd_model |
|---|
| 287 | |
|---|
| 288 | |
|---|
| 289 | !--------------------------------------------------------------------- |
|---|
| 290 | ! Read the amplitude of the spectral variations associated with the |
|---|
| 291 | ! solar cycle and map to g-points |
|---|
| 292 | subroutine read_spectral_solar_cycle(this, filename, iverbose, use_updated_solar_spectrum) |
|---|
| 293 | |
|---|
| 294 | use easy_netcdf, only : netcdf_file |
|---|
| 295 | use radiation_io, only : nulout, nulerr, radiation_abort |
|---|
| 296 | use yomhook, only : lhook, dr_hook, jphook |
|---|
| 297 | |
|---|
| 298 | ! Reference total solar irradiance (W m-2) |
|---|
| 299 | real(jprb), parameter :: ReferenceTSI = 1361.0_jprb |
|---|
| 300 | |
|---|
| 301 | class(ckd_model_type), intent(inout) :: this |
|---|
| 302 | character(len=*), intent(in) :: filename |
|---|
| 303 | integer, optional, intent(in) :: iverbose |
|---|
| 304 | ! Do we update the mean solar spectral irradiance for each g-point |
|---|
| 305 | ! based on the contents of the file? |
|---|
| 306 | logical, optional, intent(in) :: use_updated_solar_spectrum |
|---|
| 307 | |
|---|
| 308 | type(netcdf_file) :: file |
|---|
| 309 | |
|---|
| 310 | ! Solar spectral irradiance, its amplitude and wavenumber |
|---|
| 311 | ! coordinate variable, read from NetCDF file |
|---|
| 312 | real(jprb), allocatable :: wavenumber(:) ! cm-1 |
|---|
| 313 | real(jprb), allocatable :: ssi(:) ! W m-2 cm |
|---|
| 314 | real(jprb), allocatable :: ssi_amplitude(:) ! W m-2 cm |
|---|
| 315 | |
|---|
| 316 | ! As above but on the wavenumber grid delimited by |
|---|
| 317 | ! this%wavenumber1 and this%wavenumber2 |
|---|
| 318 | real(jprb), allocatable :: ssi_grid(:) |
|---|
| 319 | real(jprb), allocatable :: ssi_amplitude_grid(:) |
|---|
| 320 | real(jprb), allocatable :: wavenumber_grid(:) |
|---|
| 321 | |
|---|
| 322 | ! Old normalized solar irradiance in case it gets changed and we |
|---|
| 323 | ! need to report the amplitude of the change |
|---|
| 324 | real(jprb), allocatable :: old_norm_solar_irradiance(:) |
|---|
| 325 | |
|---|
| 326 | real(jprb) :: dwav_grid |
|---|
| 327 | |
|---|
| 328 | ! Number of input wavenumbers, and number on ecCKD model's grid |
|---|
| 329 | integer :: nwav, nwav_grid |
|---|
| 330 | ! Corresponding loop indices |
|---|
| 331 | integer :: jwav, jwav_grid, jg |
|---|
| 332 | |
|---|
| 333 | integer :: iband |
|---|
| 334 | |
|---|
| 335 | integer :: iverbose_local |
|---|
| 336 | |
|---|
| 337 | real(jphook) :: hook_handle |
|---|
| 338 | |
|---|
| 339 | if (lhook) call dr_hook('radiation_ecckd:read_spectral_solar_cycle',0,hook_handle) |
|---|
| 340 | |
|---|
| 341 | if (present(iverbose)) then |
|---|
| 342 | iverbose_local = iverbose |
|---|
| 343 | else |
|---|
| 344 | iverbose_local = 3 |
|---|
| 345 | end if |
|---|
| 346 | |
|---|
| 347 | call file%open(trim(filename), iverbose=iverbose_local) |
|---|
| 348 | |
|---|
| 349 | call file%get('wavenumber', wavenumber) |
|---|
| 350 | call file%get('mean_solar_spectral_irradiance', ssi) |
|---|
| 351 | call file%get('ssi_solar_cycle_amplitude', ssi_amplitude) |
|---|
| 352 | |
|---|
| 353 | call file%close() |
|---|
| 354 | |
|---|
| 355 | nwav = size(wavenumber) |
|---|
| 356 | |
|---|
| 357 | nwav_grid = size(this%spectral_def%wavenumber1) |
|---|
| 358 | allocate(ssi_grid(nwav_grid)) |
|---|
| 359 | allocate(ssi_amplitude_grid(nwav_grid)) |
|---|
| 360 | allocate(wavenumber_grid(nwav_grid)) |
|---|
| 361 | wavenumber_grid = 0.5_jprb * (this%spectral_def%wavenumber1+this%spectral_def%wavenumber2) |
|---|
| 362 | dwav_grid = this%spectral_def%wavenumber2(1)-this%spectral_def%wavenumber1(1) |
|---|
| 363 | |
|---|
| 364 | ssi_grid = 0.0_jprb |
|---|
| 365 | ssi_amplitude_grid = 0.0_jprb |
|---|
| 366 | |
|---|
| 367 | ! Interpolate input SSI to regular wavenumber grid |
|---|
| 368 | do jwav_grid = 1,nwav_grid |
|---|
| 369 | do jwav = 1,nwav-1 |
|---|
| 370 | if (wavenumber(jwav) < wavenumber_grid(jwav_grid) & |
|---|
| 371 | & .and. wavenumber(jwav+1) >= wavenumber_grid(jwav_grid)) then |
|---|
| 372 | ! Linear interpolation - this is not perfect |
|---|
| 373 | ssi_grid(jwav_grid) = (ssi(jwav)*(wavenumber(jwav+1)-wavenumber_grid(jwav_grid)) & |
|---|
| 374 | & +ssi(jwav+1)*(wavenumber_grid(jwav_grid)-wavenumber(jwav))) & |
|---|
| 375 | & * dwav_grid / (wavenumber(jwav+1)-wavenumber(jwav)) |
|---|
| 376 | ssi_amplitude_grid(jwav_grid) = (ssi_amplitude(jwav)*(wavenumber(jwav+1)-wavenumber_grid(jwav_grid)) & |
|---|
| 377 | & +ssi_amplitude(jwav+1)*(wavenumber_grid(jwav_grid)-wavenumber(jwav))) & |
|---|
| 378 | & * dwav_grid / (wavenumber(jwav+1)-wavenumber(jwav)) |
|---|
| 379 | exit |
|---|
| 380 | end if |
|---|
| 381 | end do |
|---|
| 382 | end do |
|---|
| 383 | |
|---|
| 384 | ! Optionally update the solar irradiances in each g-point, and the |
|---|
| 385 | ! spectral solar irradiance on the wavenumber grid corresponding |
|---|
| 386 | ! to gpoint_fraction |
|---|
| 387 | allocate(old_norm_solar_irradiance(nwav_grid)) |
|---|
| 388 | old_norm_solar_irradiance = this%norm_solar_irradiance |
|---|
| 389 | if (present(use_updated_solar_spectrum)) then |
|---|
| 390 | if (use_updated_solar_spectrum) then |
|---|
| 391 | if (.not. allocated(this%spectral_def%solar_spectral_irradiance)) then |
|---|
| 392 | write(nulerr,'(a)') 'Cannot use_updated_solar_spectrum unless gas optics model is from ecCKD >= 1.4' |
|---|
| 393 | call radiation_abort() |
|---|
| 394 | end if |
|---|
| 395 | this%norm_solar_irradiance = old_norm_solar_irradiance & |
|---|
| 396 | & * matmul(ssi_grid,this%spectral_def%gpoint_fraction) & |
|---|
| 397 | & / matmul(this%spectral_def%solar_spectral_irradiance,this%spectral_def%gpoint_fraction) |
|---|
| 398 | this%norm_solar_irradiance = this%norm_solar_irradiance / sum(this%norm_solar_irradiance) |
|---|
| 399 | this%spectral_def%solar_spectral_irradiance = ssi_grid |
|---|
| 400 | end if |
|---|
| 401 | end if |
|---|
| 402 | |
|---|
| 403 | ! Map on to g-points |
|---|
| 404 | this%norm_amplitude_solar_irradiance & |
|---|
| 405 | & = this%norm_solar_irradiance & |
|---|
| 406 | & * matmul(ssi_amplitude_grid, this%spectral_def%gpoint_fraction) & |
|---|
| 407 | & / matmul(ssi_grid,this%spectral_def%gpoint_fraction) |
|---|
| 408 | |
|---|
| 409 | ! Remove the mean from the solar-cycle fluctuations, since the |
|---|
| 410 | ! user will scale with total solar irradiance |
|---|
| 411 | this%norm_amplitude_solar_irradiance & |
|---|
| 412 | & = (this%norm_solar_irradiance+this%norm_amplitude_solar_irradiance) & |
|---|
| 413 | & / sum(this%norm_solar_irradiance+this%norm_amplitude_solar_irradiance) & |
|---|
| 414 | & - this%norm_solar_irradiance |
|---|
| 415 | |
|---|
| 416 | ! Print the spectral solar irradiance per g point, and solar cycle amplitude |
|---|
| 417 | if (iverbose_local >= 2) then |
|---|
| 418 | write(nulout,'(a,f6.1,a)') 'G-point, solar irradiance for nominal TSI = ', & |
|---|
| 419 | & ReferenceTSI, ' W m-2, solar cycle amplitude (at solar maximum), update to original solar irradiance' |
|---|
| 420 | iband = 0 |
|---|
| 421 | do jg = 1,this%ng |
|---|
| 422 | if (this%spectral_def%i_band_number(jg) > iband) then |
|---|
| 423 | iband = this%spectral_def%i_band_number(jg) |
|---|
| 424 | write(nulout, '(i2,f10.4,f7.3,a,f8.4,a,i2,a,f7.1,a,f7.1,a)') & |
|---|
| 425 | & jg, ReferenceTSI*this%norm_solar_irradiance(jg), & |
|---|
| 426 | & 100.0_jprb * this%norm_amplitude_solar_irradiance(jg) & |
|---|
| 427 | & / this%norm_solar_irradiance(jg), '% ', & |
|---|
| 428 | & 100.0_jprb * (this%norm_solar_irradiance(jg) & |
|---|
| 429 | & / old_norm_solar_irradiance(jg) - 1.0_jprb), '% Band ', iband, ': ', & |
|---|
| 430 | & this%spectral_def%wavenumber1_band(iband), '-', & |
|---|
| 431 | & this%spectral_def%wavenumber2_band(iband), ' cm-1' |
|---|
| 432 | else |
|---|
| 433 | write(nulout, '(i2,f10.4,f7.3,a,f8.4,a)') jg, ReferenceTSI*this%norm_solar_irradiance(jg), & |
|---|
| 434 | & 100.0_jprb * this%norm_amplitude_solar_irradiance(jg) & |
|---|
| 435 | & / this%norm_solar_irradiance(jg), '% ', & |
|---|
| 436 | & 100.0_jprb * (this%norm_solar_irradiance(jg) & |
|---|
| 437 | & / old_norm_solar_irradiance(jg) - 1.0_jprb), '%' |
|---|
| 438 | end if |
|---|
| 439 | end do |
|---|
| 440 | end if |
|---|
| 441 | |
|---|
| 442 | if (lhook) call dr_hook('radiation_ecckd:read_spectral_solar_cycle',1,hook_handle) |
|---|
| 443 | |
|---|
| 444 | end subroutine read_spectral_solar_cycle |
|---|
| 445 | |
|---|
| 446 | |
|---|
| 447 | !--------------------------------------------------------------------- |
|---|
| 448 | ! Compute layerwise optical depth for each g point for ncol columns |
|---|
| 449 | ! at nlev layers |
|---|
| 450 | subroutine calc_optical_depth_ckd_model(this, ncol, nlev, istartcol, iendcol, nmaxgas, & |
|---|
| 451 | & pressure_hl, temperature_fl, mole_fraction_fl, & |
|---|
| 452 | & optical_depth_fl, rayleigh_od_fl) |
|---|
| 453 | |
|---|
| 454 | use yomhook, only : lhook, dr_hook, jphook |
|---|
| 455 | use radiation_constants, only : AccelDueToGravity |
|---|
| 456 | |
|---|
| 457 | ! Input variables |
|---|
| 458 | |
|---|
| 459 | class(ckd_model_type), intent(in), target :: this |
|---|
| 460 | ! Number of columns, levels and input gases |
|---|
| 461 | integer, intent(in) :: ncol, nlev, nmaxgas, istartcol, iendcol |
|---|
| 462 | ! Pressure at half levels (Pa), dimensioned (ncol,nlev+1) |
|---|
| 463 | real(jprb), intent(in) :: pressure_hl(ncol,nlev+1) |
|---|
| 464 | ! Temperature at full levels (K), dimensioned (ncol,nlev) |
|---|
| 465 | real(jprb), intent(in) :: temperature_fl(istartcol:iendcol,nlev) |
|---|
| 466 | ! Gas mole fractions at full levels (mol mol-1), dimensioned (ncol,nlev,nmaxgas) |
|---|
| 467 | real(jprb), intent(in) :: mole_fraction_fl(ncol,nlev,nmaxgas) |
|---|
| 468 | |
|---|
| 469 | ! Output variables |
|---|
| 470 | |
|---|
| 471 | ! Layer absorption optical depth for each g point |
|---|
| 472 | real(jprb), intent(out) :: optical_depth_fl(this%ng,nlev,istartcol:iendcol) |
|---|
| 473 | ! In the shortwave only, the Rayleigh scattering optical depth |
|---|
| 474 | real(jprb), optional, intent(out) :: rayleigh_od_fl(this%ng,nlev,istartcol:iendcol) |
|---|
| 475 | |
|---|
| 476 | ! Local variables |
|---|
| 477 | |
|---|
| 478 | real(jprb), pointer :: molar_abs(:,:,:), molar_abs_conc(:,:,:,:) |
|---|
| 479 | |
|---|
| 480 | ! Natural logarithm of pressure at full levels |
|---|
| 481 | real(jprb) :: log_pressure_fl(nlev) |
|---|
| 482 | |
|---|
| 483 | ! Optical depth of single gas at one point in space versus |
|---|
| 484 | ! spectral interval |
|---|
| 485 | !real(jprb) :: od_single_gas(this%ng) |
|---|
| 486 | |
|---|
| 487 | real(jprb) :: multiplier(nlev), simple_multiplier(nlev), global_multiplier, temperature1 |
|---|
| 488 | |
|---|
| 489 | ! Indices and weights in temperature, pressure and concentration interpolation |
|---|
| 490 | real(jprb) :: pindex1, tindex1, cindex1 |
|---|
| 491 | real(jprb) :: pw1(nlev), pw2(nlev), tw1(nlev), tw2(nlev), cw1(nlev), cw2(nlev) |
|---|
| 492 | integer :: ip1(nlev), it1(nlev), ic1(nlev) |
|---|
| 493 | |
|---|
| 494 | ! Natural logarithm of mole fraction at one point |
|---|
| 495 | real(jprb) :: log_conc |
|---|
| 496 | |
|---|
| 497 | ! Minimum mole fraction in look-up-table |
|---|
| 498 | real(jprb) :: mole_frac1 |
|---|
| 499 | |
|---|
| 500 | integer :: jcol, jlev, jgas, igascode |
|---|
| 501 | |
|---|
| 502 | real(jphook) :: hook_handle |
|---|
| 503 | |
|---|
| 504 | if (lhook) call dr_hook('radiation_ecckd:calc_optical_depth',0,hook_handle) |
|---|
| 505 | |
|---|
| 506 | global_multiplier = 1.0_jprb / (AccelDueToGravity * 0.001_jprb * AirMolarMass) |
|---|
| 507 | |
|---|
| 508 | do jcol = istartcol,iendcol |
|---|
| 509 | |
|---|
| 510 | log_pressure_fl = log(0.5_jprb * (pressure_hl(jcol,1:nlev)+pressure_hl(jcol,2:nlev+1))) |
|---|
| 511 | |
|---|
| 512 | do jlev = 1,nlev |
|---|
| 513 | ! Find interpolation points in pressure |
|---|
| 514 | pindex1 = (log_pressure_fl(jlev)-this%log_pressure1) & |
|---|
| 515 | & / this%d_log_pressure |
|---|
| 516 | pindex1 = 1.0_jprb + max(0.0_jprb, min(pindex1, this%npress-1.0001_jprb)) |
|---|
| 517 | ip1(jlev) = int(pindex1) |
|---|
| 518 | pw2(jlev) = pindex1 - ip1(jlev) |
|---|
| 519 | pw1(jlev) = 1.0_jprb - pw2(jlev) |
|---|
| 520 | |
|---|
| 521 | ! Find interpolation points in temperature |
|---|
| 522 | temperature1 = pw1(jlev)*this%temperature1(ip1(jlev)) & |
|---|
| 523 | & + pw2(jlev)*this%temperature1(ip1(jlev)+1) |
|---|
| 524 | tindex1 = (temperature_fl(jcol,jlev) - temperature1) & |
|---|
| 525 | & / this%d_temperature |
|---|
| 526 | tindex1 = 1.0_jprb + max(0.0_jprb, min(tindex1, this%ntemp-1.0001_jprb)) |
|---|
| 527 | it1(jlev) = int(tindex1) |
|---|
| 528 | tw2(jlev) = tindex1 - it1(jlev) |
|---|
| 529 | tw1(jlev) = 1.0_jprb - tw2(jlev) |
|---|
| 530 | |
|---|
| 531 | ! Concentration multiplier |
|---|
| 532 | simple_multiplier(jlev) = global_multiplier & |
|---|
| 533 | & * (pressure_hl(jcol,jlev+1) - pressure_hl(jcol,jlev)) |
|---|
| 534 | end do |
|---|
| 535 | |
|---|
| 536 | optical_depth_fl(:,:,jcol) = 0.0_jprb |
|---|
| 537 | |
|---|
| 538 | do jgas = 1,this%ngas |
|---|
| 539 | |
|---|
| 540 | associate (single_gas => this%single_gas(jgas)) |
|---|
| 541 | igascode = this%single_gas(jgas)%i_gas_code |
|---|
| 542 | |
|---|
| 543 | select case (single_gas%i_conc_dependence) |
|---|
| 544 | |
|---|
| 545 | case (IConcDependenceLinear) |
|---|
| 546 | molar_abs => this%single_gas(jgas)%molar_abs |
|---|
| 547 | multiplier = simple_multiplier * mole_fraction_fl(jcol,:,igascode) |
|---|
| 548 | |
|---|
| 549 | do jlev = 1,nlev |
|---|
| 550 | optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) & |
|---|
| 551 | & + (multiplier(jlev)*tw1(jlev)) * (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)) & |
|---|
| 552 | & +pw2(jlev) * molar_abs(:,ip1(jlev)+1,it1(jlev))) & |
|---|
| 553 | & + (multiplier(jlev)*tw2(jlev)) * (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)+1) & |
|---|
| 554 | & +pw2(jlev) * molar_abs(:,ip1(jlev)+1,it1(jlev)+1)) |
|---|
| 555 | end do |
|---|
| 556 | |
|---|
| 557 | case (IConcDependenceRelativeLinear) |
|---|
| 558 | molar_abs => this%single_gas(jgas)%molar_abs |
|---|
| 559 | multiplier = simple_multiplier * (mole_fraction_fl(jcol,:,igascode) & |
|---|
| 560 | & - single_gas%reference_mole_frac) |
|---|
| 561 | do jlev = 1,nlev |
|---|
| 562 | optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) & |
|---|
| 563 | & + (multiplier(jlev)*tw1(jlev)) * (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)) & |
|---|
| 564 | & +pw2(jlev) * molar_abs(:,ip1(jlev)+1,it1(jlev))) & |
|---|
| 565 | & + (multiplier(jlev)*tw2(jlev)) * (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)+1) & |
|---|
| 566 | & +pw2(jlev) * molar_abs(:,ip1(jlev)+1,it1(jlev)+1)) |
|---|
| 567 | end do |
|---|
| 568 | |
|---|
| 569 | case (IConcDependenceNone) |
|---|
| 570 | ! Composite gases |
|---|
| 571 | molar_abs => this%single_gas(jgas)%molar_abs |
|---|
| 572 | do jlev = 1,nlev |
|---|
| 573 | optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) & |
|---|
| 574 | & + (simple_multiplier(jlev)*tw1(jlev)) * (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)) & |
|---|
| 575 | & +pw2(jlev) * molar_abs(:,ip1(jlev)+1,it1(jlev))) & |
|---|
| 576 | & + (simple_multiplier(jlev)*tw2(jlev)) * (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)+1) & |
|---|
| 577 | & +pw2(jlev) * molar_abs(:,ip1(jlev)+1,it1(jlev)+1)) |
|---|
| 578 | end do |
|---|
| 579 | |
|---|
| 580 | case (IConcDependenceLUT) |
|---|
| 581 | ! Logarithmic interpolation in concentration space |
|---|
| 582 | molar_abs_conc => this%single_gas(jgas)%molar_abs_conc |
|---|
| 583 | mole_frac1 = exp(single_gas%log_mole_frac1) |
|---|
| 584 | do jlev = 1,nlev |
|---|
| 585 | ! Take care of mole_fraction == 0 |
|---|
| 586 | log_conc = log(max(mole_fraction_fl(jcol,jlev,igascode), mole_frac1)) |
|---|
| 587 | cindex1 = (log_conc - single_gas%log_mole_frac1) / single_gas%d_log_mole_frac |
|---|
| 588 | cindex1 = 1.0_jprb + max(0.0_jprb, min(cindex1, single_gas%n_mole_frac-1.0001_jprb)) |
|---|
| 589 | ic1(jlev) = int(cindex1) |
|---|
| 590 | cw2(jlev) = cindex1 - ic1(jlev) |
|---|
| 591 | cw1(jlev) = 1.0_jprb - cw2(jlev) |
|---|
| 592 | end do |
|---|
| 593 | ! od_single_gas = cw1 * (tw1 * (pw1 * molar_abs_conc(:,ip1,it1,ic1) & |
|---|
| 594 | ! & +pw2 * molar_abs_conc(:,ip1+1,it1,ic1)) & |
|---|
| 595 | ! & +tw2 * (pw1 * molar_abs_conc(:,ip1,it1+1,ic1) & |
|---|
| 596 | ! & +pw2 * molar_abs_conc(:,ip1+1,it1+1,ic1))) & |
|---|
| 597 | ! & +cw2 * (tw1 * (pw1 * molar_abs_conc(:,ip1,it1,ic1+1) & |
|---|
| 598 | ! & +pw2 * molar_abs_conc(:,ip1+1,it1,ic1+1)) & |
|---|
| 599 | ! & +tw2 * (pw1 * molar_abs_conc(:,ip1,it1+1,ic1+1) & |
|---|
| 600 | ! & +pw2 * molar_abs_conc(:,ip1+1,it1+1,ic1+1))) |
|---|
| 601 | do jlev = 1,nlev |
|---|
| 602 | optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) & |
|---|
| 603 | & + (simple_multiplier(jlev) * mole_fraction_fl(jcol,jlev,igascode)) * ( & |
|---|
| 604 | & (cw1(jlev) * tw1(jlev) * pw1(jlev)) * molar_abs_conc(:,ip1(jlev),it1(jlev),ic1(jlev)) & |
|---|
| 605 | & +(cw1(jlev) * tw1(jlev) * pw2(jlev)) * molar_abs_conc(:,ip1(jlev)+1,it1(jlev),ic1(jlev)) & |
|---|
| 606 | & +(cw1(jlev) * tw2(jlev) * pw1(jlev)) * molar_abs_conc(:,ip1(jlev),it1(jlev)+1,ic1(jlev)) & |
|---|
| 607 | & +(cw1(jlev) * tw2(jlev) * pw2(jlev)) * molar_abs_conc(:,ip1(jlev)+1,it1(jlev)+1,ic1(jlev)) & |
|---|
| 608 | & +(cw2(jlev) * tw1(jlev) * pw1(jlev)) * molar_abs_conc(:,ip1(jlev),it1(jlev),ic1(jlev)+1) & |
|---|
| 609 | & +(cw2(jlev) * tw1(jlev) * pw2(jlev)) * molar_abs_conc(:,ip1(jlev)+1,it1(jlev),ic1(jlev)+1) & |
|---|
| 610 | & +(cw2(jlev) * tw2(jlev) * pw1(jlev)) * molar_abs_conc(:,ip1(jlev),it1(jlev)+1,ic1(jlev)+1) & |
|---|
| 611 | & +(cw2(jlev) * tw2(jlev) * pw2(jlev)) * molar_abs_conc(:,ip1(jlev)+1,it1(jlev)+1,ic1(jlev)+1)) |
|---|
| 612 | end do |
|---|
| 613 | end select |
|---|
| 614 | |
|---|
| 615 | end associate |
|---|
| 616 | |
|---|
| 617 | end do |
|---|
| 618 | |
|---|
| 619 | ! Ensure the optical depth is not negative |
|---|
| 620 | optical_depth_fl(:,:,jcol) = max(0.0_jprb, optical_depth_fl(:,:,jcol)) |
|---|
| 621 | |
|---|
| 622 | ! Rayleigh scattering |
|---|
| 623 | if (this%is_sw .and. present(rayleigh_od_fl)) then |
|---|
| 624 | do jlev = 1,nlev |
|---|
| 625 | rayleigh_od_fl(:,jlev,jcol) = global_multiplier & |
|---|
| 626 | & * (pressure_hl(jcol,jlev+1) - pressure_hl(jcol,jlev)) * this%rayleigh_molar_scat |
|---|
| 627 | end do |
|---|
| 628 | end if |
|---|
| 629 | |
|---|
| 630 | end do |
|---|
| 631 | |
|---|
| 632 | if (lhook) call dr_hook('radiation_ecckd:calc_optical_depth',1,hook_handle) |
|---|
| 633 | |
|---|
| 634 | end subroutine calc_optical_depth_ckd_model |
|---|
| 635 | |
|---|
| 636 | |
|---|
| 637 | !--------------------------------------------------------------------- |
|---|
| 638 | ! Vectorized variant of above routine |
|---|
| 639 | subroutine calc_optical_depth_ckd_model_vec(this, ncol, nlev, istartcol, iendcol, nmaxgas, & |
|---|
| 640 | & pressure_hl, temperature_fl, mole_fraction_fl, & |
|---|
| 641 | & optical_depth_fl, rayleigh_od_fl) |
|---|
| 642 | |
|---|
| 643 | use yomhook, only : lhook, dr_hook, jphook |
|---|
| 644 | use radiation_constants, only : AccelDueToGravity |
|---|
| 645 | |
|---|
| 646 | ! Input variables |
|---|
| 647 | |
|---|
| 648 | class(ckd_model_type), intent(in), target :: this |
|---|
| 649 | ! Number of columns, levels and input gases |
|---|
| 650 | integer, intent(in) :: ncol, nlev, nmaxgas, istartcol, iendcol |
|---|
| 651 | ! Pressure at half levels (Pa), dimensioned (ncol,nlev+1) |
|---|
| 652 | real(jprb), intent(in) :: pressure_hl(ncol,nlev+1) |
|---|
| 653 | ! Temperature at full levels (K), dimensioned (ncol,nlev) |
|---|
| 654 | real(jprb), intent(in) :: temperature_fl(istartcol:iendcol,nlev) |
|---|
| 655 | ! Gas mole fractions at full levels (mol mol-1), dimensioned (ncol,nlev,nmaxgas) |
|---|
| 656 | real(jprb), intent(in) :: mole_fraction_fl(ncol,nlev,nmaxgas) |
|---|
| 657 | |
|---|
| 658 | ! Output variables |
|---|
| 659 | |
|---|
| 660 | ! Layer absorption optical depth for each g point |
|---|
| 661 | real(jprb), intent(out) :: optical_depth_fl(this%ng,nlev,istartcol:iendcol) |
|---|
| 662 | ! In the shortwave only, the Rayleigh scattering optical depth |
|---|
| 663 | real(jprb), optional, intent(out) :: rayleigh_od_fl(this%ng,nlev,istartcol:iendcol) |
|---|
| 664 | |
|---|
| 665 | ! Local variables |
|---|
| 666 | |
|---|
| 667 | real(jprb), pointer :: molar_abs(:,:,:), molar_abs_conc(:,:,:,:) |
|---|
| 668 | |
|---|
| 669 | ! Natural logarithm of pressure at full levels |
|---|
| 670 | real(jprb) :: log_pressure_fl |
|---|
| 671 | |
|---|
| 672 | ! Optical depth of single gas at one point in space versus |
|---|
| 673 | ! spectral interval |
|---|
| 674 | !real(jprb) :: od_single_gas(this%ng) |
|---|
| 675 | |
|---|
| 676 | real(jprb) :: multiplier, simple_multiplier(ncol,nlev), global_multiplier, temperature1 |
|---|
| 677 | |
|---|
| 678 | ! Indices and weights in temperature, pressure and concentration interpolation |
|---|
| 679 | real(jprb) :: pindex1, tindex1, cindex1 |
|---|
| 680 | real(jprb) :: pw1(ncol,nlev), pw2(ncol,nlev), tw1(ncol,nlev), tw2(ncol,nlev), cw1(ncol,nlev), cw2(ncol,nlev) |
|---|
| 681 | integer :: ip1(ncol,nlev), it1(ncol,nlev), ic1(ncol,nlev) |
|---|
| 682 | |
|---|
| 683 | ! Natural logarithm of mole fraction at one point |
|---|
| 684 | real(jprb) :: log_conc |
|---|
| 685 | |
|---|
| 686 | ! Minimum mole fraction in look-up-table |
|---|
| 687 | real(jprb) :: mole_frac1 |
|---|
| 688 | |
|---|
| 689 | ! Layer absorption optical depth for each g point (memory layout adjusted to vectorization) |
|---|
| 690 | real(jprb) :: od_fl(ncol,this%ng,nlev) |
|---|
| 691 | |
|---|
| 692 | integer :: jcol, jlev, jgas, igascode, jg |
|---|
| 693 | |
|---|
| 694 | real(jphook) :: hook_handle |
|---|
| 695 | |
|---|
| 696 | if (lhook) call dr_hook('radiation_ecckd:calc_optical_depth_vec',0,hook_handle) |
|---|
| 697 | |
|---|
| 698 | global_multiplier = 1.0_jprb / (AccelDueToGravity * 0.001_jprb * AirMolarMass) |
|---|
| 699 | |
|---|
| 700 | od_fl(:,:,:) = 0.0_jprb |
|---|
| 701 | |
|---|
| 702 | do jlev = 1,nlev |
|---|
| 703 | do jcol = istartcol,iendcol |
|---|
| 704 | |
|---|
| 705 | log_pressure_fl = log(0.5_jprb * (pressure_hl(jcol,jlev)+pressure_hl(jcol,jlev+1))) |
|---|
| 706 | |
|---|
| 707 | ! Find interpolation points in pressure |
|---|
| 708 | pindex1 = (log_pressure_fl-this%log_pressure1) & |
|---|
| 709 | & / this%d_log_pressure |
|---|
| 710 | pindex1 = 1.0_jprb + max(0.0_jprb, min(pindex1, this%npress-1.0001_jprb)) |
|---|
| 711 | ip1(jcol,jlev) = int(pindex1) |
|---|
| 712 | pw2(jcol,jlev) = pindex1 - ip1(jcol,jlev) |
|---|
| 713 | pw1(jcol,jlev) = 1.0_jprb - pw2(jcol,jlev) |
|---|
| 714 | |
|---|
| 715 | ! Find interpolation points in temperature |
|---|
| 716 | temperature1 = pw1(jcol,jlev)*this%temperature1(ip1(jcol,jlev)) & |
|---|
| 717 | & + pw2(jcol,jlev)*this%temperature1(ip1(jcol,jlev)+1) |
|---|
| 718 | tindex1 = (temperature_fl(jcol,jlev) - temperature1) & |
|---|
| 719 | & / this%d_temperature |
|---|
| 720 | tindex1 = 1.0_jprb + max(0.0_jprb, min(tindex1, this%ntemp-1.0001_jprb)) |
|---|
| 721 | it1(jcol,jlev) = int(tindex1) |
|---|
| 722 | tw2(jcol,jlev) = tindex1 - it1(jcol,jlev) |
|---|
| 723 | tw1(jcol,jlev) = 1.0_jprb - tw2(jcol,jlev) |
|---|
| 724 | |
|---|
| 725 | ! Concentration multiplier |
|---|
| 726 | simple_multiplier(jcol,jlev) = global_multiplier & |
|---|
| 727 | & * (pressure_hl(jcol,jlev+1) - pressure_hl(jcol,jlev)) |
|---|
| 728 | end do |
|---|
| 729 | end do |
|---|
| 730 | |
|---|
| 731 | do jgas = 1,this%ngas |
|---|
| 732 | |
|---|
| 733 | associate (single_gas => this%single_gas(jgas)) |
|---|
| 734 | igascode = this%single_gas(jgas)%i_gas_code |
|---|
| 735 | |
|---|
| 736 | select case (single_gas%i_conc_dependence) |
|---|
| 737 | |
|---|
| 738 | case (IConcDependenceLinear) |
|---|
| 739 | molar_abs => this%single_gas(jgas)%molar_abs |
|---|
| 740 | |
|---|
| 741 | do jlev = 1,nlev |
|---|
| 742 | do jg = 1, this%ng |
|---|
| 743 | do jcol = istartcol,iendcol |
|---|
| 744 | multiplier = simple_multiplier(jcol,jlev) * mole_fraction_fl(jcol,jlev,igascode) |
|---|
| 745 | |
|---|
| 746 | od_fl(jcol,jg,jlev) = od_fl(jcol,jg,jlev) & |
|---|
| 747 | & + (multiplier*tw1(jcol,jlev)) * (pw1(jcol,jlev) * molar_abs(jg,ip1(jcol,jlev),it1(jcol,jlev)) & |
|---|
| 748 | & +pw2(jcol,jlev) * molar_abs(jg,ip1(jcol,jlev)+1,it1(jcol,jlev))) & |
|---|
| 749 | & + (multiplier*tw2(jcol,jlev)) * (pw1(jcol,jlev) * molar_abs(jg,ip1(jcol,jlev),it1(jcol,jlev)+1) & |
|---|
| 750 | & +pw2(jcol,jlev) * molar_abs(jg,ip1(jcol,jlev)+1,it1(jcol,jlev)+1)) |
|---|
| 751 | end do |
|---|
| 752 | end do |
|---|
| 753 | end do |
|---|
| 754 | |
|---|
| 755 | case (IConcDependenceRelativeLinear) |
|---|
| 756 | molar_abs => this%single_gas(jgas)%molar_abs |
|---|
| 757 | |
|---|
| 758 | do jlev = 1,nlev |
|---|
| 759 | do jg = 1, this%ng |
|---|
| 760 | do jcol = istartcol,iendcol |
|---|
| 761 | multiplier = simple_multiplier(jcol,jlev) * (mole_fraction_fl(jcol,jlev,igascode) & |
|---|
| 762 | & - single_gas%reference_mole_frac) |
|---|
| 763 | |
|---|
| 764 | od_fl(jcol,jg,jlev) = od_fl(jcol,jg,jlev) & |
|---|
| 765 | & + (multiplier*tw1(jcol,jlev)) * (pw1(jcol,jlev) * molar_abs(jg,ip1(jcol,jlev),it1(jcol,jlev)) & |
|---|
| 766 | & +pw2(jcol,jlev) * molar_abs(jg,ip1(jcol,jlev)+1,it1(jcol,jlev))) & |
|---|
| 767 | & + (multiplier*tw2(jcol,jlev)) * (pw1(jcol,jlev) * molar_abs(jg,ip1(jcol,jlev),it1(jcol,jlev)+1) & |
|---|
| 768 | & +pw2(jcol,jlev) * molar_abs(jg,ip1(jcol,jlev)+1,it1(jcol,jlev)+1)) |
|---|
| 769 | end do |
|---|
| 770 | end do |
|---|
| 771 | end do |
|---|
| 772 | |
|---|
| 773 | case (IConcDependenceNone) |
|---|
| 774 | ! Composite gases |
|---|
| 775 | molar_abs => this%single_gas(jgas)%molar_abs |
|---|
| 776 | |
|---|
| 777 | do jlev = 1,nlev |
|---|
| 778 | do jg = 1, this%ng |
|---|
| 779 | do jcol = istartcol,iendcol |
|---|
| 780 | od_fl(jcol,jg,jlev) = od_fl(jcol,jg,jlev) & |
|---|
| 781 | & + (simple_multiplier(jcol,jlev)*tw1(jcol,jlev)) * & |
|---|
| 782 | & (pw1(jcol,jlev) * molar_abs(jg,ip1(jcol,jlev),it1(jcol,jlev)) & |
|---|
| 783 | & +pw2(jcol,jlev) * molar_abs(jg,ip1(jcol,jlev)+1,it1(jcol,jlev))) & |
|---|
| 784 | & + (simple_multiplier(jcol,jlev)*tw2(jcol,jlev)) * & |
|---|
| 785 | & (pw1(jcol,jlev) * molar_abs(jg,ip1(jcol,jlev),it1(jcol,jlev)+1) & |
|---|
| 786 | & +pw2(jcol,jlev) * molar_abs(jg,ip1(jcol,jlev)+1,it1(jcol,jlev)+1)) |
|---|
| 787 | end do |
|---|
| 788 | end do |
|---|
| 789 | end do |
|---|
| 790 | |
|---|
| 791 | case (IConcDependenceLUT) |
|---|
| 792 | ! Logarithmic interpolation in concentration space |
|---|
| 793 | molar_abs_conc => this%single_gas(jgas)%molar_abs_conc |
|---|
| 794 | mole_frac1 = exp(single_gas%log_mole_frac1) |
|---|
| 795 | |
|---|
| 796 | do jlev = 1,nlev |
|---|
| 797 | do jcol = istartcol,iendcol |
|---|
| 798 | ! Take care of mole_fraction == 0 |
|---|
| 799 | log_conc = log(max(mole_fraction_fl(jcol,jlev,igascode), mole_frac1)) |
|---|
| 800 | cindex1 = (log_conc - single_gas%log_mole_frac1) / single_gas%d_log_mole_frac |
|---|
| 801 | cindex1 = 1.0_jprb + max(0.0_jprb, min(cindex1, single_gas%n_mole_frac-1.0001_jprb)) |
|---|
| 802 | ic1(jcol,jlev) = int(cindex1) |
|---|
| 803 | cw2(jcol,jlev) = cindex1 - ic1(jcol,jlev) |
|---|
| 804 | cw1(jcol,jlev) = 1.0_jprb - cw2(jcol,jlev) |
|---|
| 805 | end do |
|---|
| 806 | end do |
|---|
| 807 | |
|---|
| 808 | do jlev = 1,nlev |
|---|
| 809 | do jg = 1, this%ng |
|---|
| 810 | !NEC$ select_vector |
|---|
| 811 | do jcol = istartcol,iendcol |
|---|
| 812 | |
|---|
| 813 | od_fl(jcol,jg,jlev) = od_fl(jcol,jg,jlev) & |
|---|
| 814 | & + (simple_multiplier(jcol,jlev) * mole_fraction_fl(jcol,jlev,igascode)) * ( & |
|---|
| 815 | & (cw1(jcol,jlev) * tw1(jcol,jlev) * pw1(jcol,jlev)) * & |
|---|
| 816 | & molar_abs_conc(jg,ip1(jcol,jlev),it1(jcol,jlev),ic1(jcol,jlev)) & |
|---|
| 817 | & +(cw1(jcol,jlev) * tw1(jcol,jlev) * pw2(jcol,jlev)) * & |
|---|
| 818 | & molar_abs_conc(jg,ip1(jcol,jlev)+1,it1(jcol,jlev),ic1(jcol,jlev)) & |
|---|
| 819 | & +(cw1(jcol,jlev) * tw2(jcol,jlev) * pw1(jcol,jlev)) * & |
|---|
| 820 | & molar_abs_conc(jg,ip1(jcol,jlev),it1(jcol,jlev)+1,ic1(jcol,jlev)) & |
|---|
| 821 | & +(cw1(jcol,jlev) * tw2(jcol,jlev) * pw2(jcol,jlev)) * & |
|---|
| 822 | & molar_abs_conc(jg,ip1(jcol,jlev)+1,it1(jcol,jlev)+1,ic1(jcol,jlev)) & |
|---|
| 823 | & +(cw2(jcol,jlev) * tw1(jcol,jlev) * pw1(jcol,jlev)) * & |
|---|
| 824 | & molar_abs_conc(jg,ip1(jcol,jlev),it1(jcol,jlev),ic1(jcol,jlev)+1) & |
|---|
| 825 | & +(cw2(jcol,jlev) * tw1(jcol,jlev) * pw2(jcol,jlev)) * & |
|---|
| 826 | & molar_abs_conc(jg,ip1(jcol,jlev)+1,it1(jcol,jlev),ic1(jcol,jlev)+1) & |
|---|
| 827 | & +(cw2(jcol,jlev) * tw2(jcol,jlev) * pw1(jcol,jlev)) * & |
|---|
| 828 | & molar_abs_conc(jg,ip1(jcol,jlev),it1(jcol,jlev)+1,ic1(jcol,jlev)+1) & |
|---|
| 829 | & +(cw2(jcol,jlev) * tw2(jcol,jlev) * pw2(jcol,jlev)) * & |
|---|
| 830 | & molar_abs_conc(jg,ip1(jcol,jlev)+1,it1(jcol,jlev)+1,ic1(jcol,jlev)+1)) |
|---|
| 831 | end do |
|---|
| 832 | end do |
|---|
| 833 | end do |
|---|
| 834 | end select |
|---|
| 835 | |
|---|
| 836 | end associate |
|---|
| 837 | |
|---|
| 838 | ! Ensure the optical depth is not negative |
|---|
| 839 | do jcol = istartcol,iendcol |
|---|
| 840 | do jlev = 1,nlev |
|---|
| 841 | do jg = 1, this%ng |
|---|
| 842 | optical_depth_fl(jg,jlev,jcol) = max(0.0_jprb, od_fl(jcol,jg,jlev)) |
|---|
| 843 | end do |
|---|
| 844 | end do |
|---|
| 845 | end do |
|---|
| 846 | |
|---|
| 847 | ! Rayleigh scattering |
|---|
| 848 | if (this%is_sw .and. present(rayleigh_od_fl)) then |
|---|
| 849 | do jcol = istartcol,iendcol |
|---|
| 850 | do jlev = 1,nlev |
|---|
| 851 | do jg = 1, this%ng |
|---|
| 852 | rayleigh_od_fl(jg,jlev,jcol) = global_multiplier & |
|---|
| 853 | & * (pressure_hl(jcol,jlev+1) - pressure_hl(jcol,jlev)) * this%rayleigh_molar_scat(jg) |
|---|
| 854 | end do |
|---|
| 855 | end do |
|---|
| 856 | end do |
|---|
| 857 | end if |
|---|
| 858 | |
|---|
| 859 | end do |
|---|
| 860 | |
|---|
| 861 | if (lhook) call dr_hook('radiation_ecckd:calc_optical_depth_vec',1,hook_handle) |
|---|
| 862 | |
|---|
| 863 | end subroutine calc_optical_depth_ckd_model_vec |
|---|
| 864 | |
|---|
| 865 | |
|---|
| 866 | !--------------------------------------------------------------------- |
|---|
| 867 | ! Calculate the Planck function integrated across each of the g |
|---|
| 868 | ! points of this correlated k-distribution model, for a given |
|---|
| 869 | ! temperature, where Planck function is defined as the flux emitted |
|---|
| 870 | ! by a black body (rather than radiance) |
|---|
| 871 | subroutine calc_planck_function(this, nt, temperature, planck) |
|---|
| 872 | |
|---|
| 873 | class(ckd_model_type), intent(in) :: this |
|---|
| 874 | integer, intent(in) :: nt |
|---|
| 875 | real(jprb), intent(in) :: temperature(:) ! K |
|---|
| 876 | real(jprb), intent(out) :: planck(this%ng,nt) ! W m-2 |
|---|
| 877 | |
|---|
| 878 | real(jprb) :: tindex1, tw1, tw2 |
|---|
| 879 | integer :: it1, jt |
|---|
| 880 | |
|---|
| 881 | do jt = 1,nt |
|---|
| 882 | tindex1 = (temperature(jt) - this%temperature1_planck) & |
|---|
| 883 | & * (1.0_jprb / this%d_temperature_planck) |
|---|
| 884 | if (tindex1 >= 0) then |
|---|
| 885 | ! Normal interpolation, and extrapolation for high temperatures |
|---|
| 886 | tindex1 = 1.0_jprb + tindex1 |
|---|
| 887 | it1 = min(int(tindex1), this%nplanck-1) |
|---|
| 888 | tw2 = tindex1 - it1 |
|---|
| 889 | tw1 = 1.0_jprb - tw2 |
|---|
| 890 | planck(:,jt) = tw1 * this%planck_function(:,it1) & |
|---|
| 891 | & + tw2 * this%planck_function(:,it1+1) |
|---|
| 892 | else |
|---|
| 893 | ! Interpolate linearly to zero |
|---|
| 894 | planck(:,jt) = this%planck_function(:,1) & |
|---|
| 895 | & * (temperature(jt)/this%temperature1_planck) |
|---|
| 896 | end if |
|---|
| 897 | end do |
|---|
| 898 | |
|---|
| 899 | end subroutine calc_planck_function |
|---|
| 900 | |
|---|
| 901 | |
|---|
| 902 | !--------------------------------------------------------------------- |
|---|
| 903 | ! Return the spectral solar irradiance integrated over each g point |
|---|
| 904 | ! of a solar correlated k-distribution model, given the |
|---|
| 905 | ! total_solar_irradiance |
|---|
| 906 | subroutine calc_incoming_sw(this, total_solar_irradiance, & |
|---|
| 907 | & spectral_solar_irradiance, & |
|---|
| 908 | & solar_spectral_multiplier) |
|---|
| 909 | |
|---|
| 910 | use radiation_io, only : nulerr, radiation_abort |
|---|
| 911 | |
|---|
| 912 | class(ckd_model_type), intent(in) :: this |
|---|
| 913 | real(jprb), intent(in) :: total_solar_irradiance ! W m-2 |
|---|
| 914 | real(jprb), intent(inout) :: spectral_solar_irradiance(:,:) ! W m-2 |
|---|
| 915 | real(jprb), optional, intent(in) :: solar_spectral_multiplier |
|---|
| 916 | |
|---|
| 917 | if (.not. present(solar_spectral_multiplier)) then |
|---|
| 918 | spectral_solar_irradiance & |
|---|
| 919 | & = spread(total_solar_irradiance * this%norm_solar_irradiance, & |
|---|
| 920 | & 2, size(spectral_solar_irradiance,2)) |
|---|
| 921 | else if (allocated(this%norm_amplitude_solar_irradiance)) then |
|---|
| 922 | spectral_solar_irradiance & |
|---|
| 923 | & = spread(total_solar_irradiance * (this%norm_solar_irradiance & |
|---|
| 924 | & + solar_spectral_multiplier*this%norm_amplitude_solar_irradiance), & |
|---|
| 925 | & 2, size(spectral_solar_irradiance,2)) |
|---|
| 926 | else if (solar_spectral_multiplier == 0.0_jprb) then |
|---|
| 927 | spectral_solar_irradiance & |
|---|
| 928 | & = spread(total_solar_irradiance * this%norm_solar_irradiance, & |
|---|
| 929 | & 2, size(spectral_solar_irradiance,2)) |
|---|
| 930 | else |
|---|
| 931 | write(nulerr, '(a)') '*** Error in calc_incoming_sw: no information present on solar cycle' |
|---|
| 932 | call radiation_abort() |
|---|
| 933 | end if |
|---|
| 934 | |
|---|
| 935 | end subroutine calc_incoming_sw |
|---|
| 936 | |
|---|
| 937 | end module radiation_ecckd |
|---|
