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