[3908] | 1 | ! radiation_flux.F90 - Derived type to store the output fluxes |
---|
| 2 | ! |
---|
| 3 | ! (C) Copyright 2014- ECMWF. |
---|
| 4 | ! |
---|
| 5 | ! This software is licensed under the terms of the Apache Licence Version 2.0 |
---|
| 6 | ! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. |
---|
| 7 | ! |
---|
| 8 | ! In applying this licence, ECMWF does not waive the privileges and immunities |
---|
| 9 | ! granted to it by virtue of its status as an intergovernmental organisation |
---|
| 10 | ! nor does it submit to any jurisdiction. |
---|
| 11 | ! |
---|
| 12 | ! Author: Robin Hogan |
---|
| 13 | ! Email: r.j.hogan@ecmwf.int |
---|
| 14 | ! |
---|
| 15 | ! Modifications |
---|
| 16 | ! 2017-09-08 R. Hogan Store g-point fluxes |
---|
| 17 | ! 2017-10-23 R. Hogan Renamed single-character variables |
---|
| 18 | ! 2019-01-08 R. Hogan Added "indexed_sum_profile" |
---|
| 19 | ! 2019-01-14 R. Hogan out_of_physical_bounds calls routine in radiation_config |
---|
| 20 | ! 2021-01-20 R. Hogan Added heating_rate_out_of_physical_bounds function |
---|
| 21 | |
---|
| 22 | module radiation_flux |
---|
| 23 | |
---|
| 24 | use parkind1, only : jprb |
---|
| 25 | |
---|
| 26 | implicit none |
---|
| 27 | public |
---|
| 28 | |
---|
| 29 | !--------------------------------------------------------------------- |
---|
| 30 | ! This derived type contains the output from the radiation |
---|
| 31 | ! calculation. Currently this is solely flux profiles, but in |
---|
| 32 | ! future surface fluxes in each band may be stored in order that the |
---|
| 33 | ! calling program can compute surface-radiation such as |
---|
| 34 | ! photosynthetically active radiation and UV index. |
---|
| 35 | type flux_type |
---|
| 36 | ! All the following are broad-band fluxes in W m-2 with |
---|
| 37 | ! dimensions (ncol,nlev+1). Note that only those fluxes that are |
---|
| 38 | ! requested will be used, so clear-sky and direct-beam arrays may |
---|
| 39 | ! not be allocated |
---|
| 40 | real(jprb), allocatable, dimension(:,:) :: & |
---|
| 41 | & lw_up, lw_dn, & ! Upwelling and downwelling longwave |
---|
| 42 | & sw_up, sw_dn, & ! Upwelling and downwelling shortwave |
---|
| 43 | & sw_dn_direct, & ! Direct-beam shortwave into a horizontal plane |
---|
| 44 | & lw_up_clear, lw_dn_clear, & ! Clear-sky quantities... |
---|
| 45 | & sw_up_clear, sw_dn_clear, & |
---|
| 46 | & sw_dn_direct_clear |
---|
| 47 | ! As above but fluxes in each spectral band in W m-2 with |
---|
| 48 | ! dimensions (nband,ncol,nlev+1). These are only allocated if |
---|
| 49 | ! config%do_save_spectral_flux==.true. |
---|
| 50 | real(jprb), allocatable, dimension(:,:,:) :: & |
---|
| 51 | & lw_up_band, lw_dn_band, & ! Upwelling and downwelling longwave |
---|
| 52 | & sw_up_band, sw_dn_band, & ! Upwelling and downwelling shortwave |
---|
| 53 | & sw_dn_direct_band, & ! Direct-beam shortwave |
---|
| 54 | & lw_up_clear_band, lw_dn_clear_band, & ! Clear-sky quantities... |
---|
| 55 | & sw_up_clear_band, sw_dn_clear_band, & |
---|
| 56 | & sw_dn_direct_clear_band |
---|
| 57 | ! Surface downwelling quantaties at each g point, dimensioned |
---|
| 58 | ! (ng,ncol), that are always saved by the solver, except for the |
---|
| 59 | ! clear-sky ones that are only produced if |
---|
| 60 | ! config%do_clear==.true. |
---|
| 61 | real(jprb), allocatable, dimension(:,:) :: & |
---|
| 62 | & lw_dn_surf_g, lw_dn_surf_clear_g, & |
---|
| 63 | & sw_dn_diffuse_surf_g, sw_dn_direct_surf_g, & |
---|
| 64 | & sw_dn_diffuse_surf_clear_g, sw_dn_direct_surf_clear_g |
---|
| 65 | ! Shortwave downwelling spectral fluxes in W m-2 at the surface, |
---|
| 66 | ! from which quantities such as photosynthetically active and UV |
---|
| 67 | ! radiation can be computed. Only allocated in |
---|
| 68 | ! config%do_surface_sw_spectral_flux==.true. Note that the |
---|
| 69 | ! clear-sky quantities are only computed if |
---|
| 70 | ! config%do_clear==.true., but direct fluxes are computed whether |
---|
| 71 | ! or not do_direct==.true.. The dimensions are (nband,ncol). |
---|
| 72 | real(jprb), allocatable, dimension(:,:) :: & |
---|
| 73 | & sw_dn_surf_band, sw_dn_direct_surf_band, & |
---|
| 74 | & sw_dn_surf_clear_band, sw_dn_direct_surf_clear_band |
---|
| 75 | ! Surface downwelling fluxes in W m-2 at the spectral resolution |
---|
| 76 | ! needed by any subsequent canopy radiative transfer. If |
---|
| 77 | ! config%use_canopy_full_spectrum_[sw|lw] then these will be at |
---|
| 78 | ! g-point resolution; otherwise they will be at |
---|
| 79 | ! config%n_albedo_bands and config%n_emiss_bands resolution. |
---|
| 80 | real(jprb), allocatable, dimension(:,:) :: & |
---|
| 81 | & lw_dn_surf_canopy, & |
---|
| 82 | & sw_dn_diffuse_surf_canopy, sw_dn_direct_surf_canopy |
---|
| 83 | |
---|
| 84 | ! Diagnosed cloud cover from the short- and long-wave solvers |
---|
| 85 | real(jprb), allocatable, dimension(:) :: & |
---|
| 86 | & cloud_cover_lw, cloud_cover_sw |
---|
| 87 | ! Longwave derivatives needed by Hogan and Bozzo (2015) method |
---|
| 88 | ! for approximate longwave updates in between the full radiation |
---|
| 89 | ! calls: rate of change of upwelling broad-band flux with respect |
---|
| 90 | ! to surface value, dimensioned (ncol,nlev+1) |
---|
| 91 | real(jprb), allocatable, dimension(:,:) :: & |
---|
| 92 | & lw_derivatives |
---|
| 93 | |
---|
| 94 | contains |
---|
| 95 | procedure :: allocate => allocate_flux_type |
---|
| 96 | procedure :: deallocate => deallocate_flux_type |
---|
| 97 | procedure :: calc_surface_spectral |
---|
| 98 | procedure :: out_of_physical_bounds |
---|
| 99 | procedure :: heating_rate_out_of_physical_bounds |
---|
| 100 | end type flux_type |
---|
| 101 | |
---|
[4584] | 102 | ! Added for DWD (2020) |
---|
| 103 | #ifdef __SX__ |
---|
| 104 | logical, parameter :: use_indexed_sum_vec = .true. |
---|
| 105 | #else |
---|
| 106 | logical, parameter :: use_indexed_sum_vec = .false. |
---|
| 107 | #endif |
---|
| 108 | |
---|
[3908] | 109 | contains |
---|
| 110 | |
---|
| 111 | !--------------------------------------------------------------------- |
---|
| 112 | ! Allocate arrays for flux profiles, using config to define which |
---|
| 113 | ! fluxes are needed. The arrays are dimensioned for columns between |
---|
| 114 | ! istartcol, iendcol and levels from 1 to nlev+1 |
---|
| 115 | subroutine allocate_flux_type(this, config, istartcol, iendcol, nlev) |
---|
| 116 | |
---|
| 117 | use yomhook, only : lhook, dr_hook |
---|
| 118 | use radiation_io, only : nulerr, radiation_abort |
---|
| 119 | use radiation_config, only : config_type |
---|
| 120 | |
---|
| 121 | integer, intent(in) :: istartcol, iendcol, nlev |
---|
| 122 | class(flux_type), intent(inout) :: this |
---|
| 123 | type(config_type), intent(in) :: config |
---|
| 124 | |
---|
| 125 | real(jprb) :: hook_handle |
---|
| 126 | |
---|
| 127 | if (lhook) call dr_hook('radiation_flux:allocate',0,hook_handle) |
---|
| 128 | |
---|
| 129 | ! Allocate longwave arrays |
---|
| 130 | if (config%do_lw) then |
---|
| 131 | allocate(this%lw_up(istartcol:iendcol,nlev+1)) |
---|
| 132 | allocate(this%lw_dn(istartcol:iendcol,nlev+1)) |
---|
| 133 | if (config%do_clear) then |
---|
| 134 | allocate(this%lw_up_clear(istartcol:iendcol,nlev+1)) |
---|
| 135 | allocate(this%lw_dn_clear(istartcol:iendcol,nlev+1)) |
---|
| 136 | end if |
---|
| 137 | |
---|
| 138 | if (config%do_save_spectral_flux) then |
---|
| 139 | if (config%n_spec_lw == 0) then |
---|
| 140 | write(nulerr,'(a)') '*** Error: number of LW spectral points to save not yet defined ' & |
---|
[4584] | 141 | & // 'so cannot allocate spectral flux arrays' |
---|
[3908] | 142 | call radiation_abort() |
---|
| 143 | end if |
---|
| 144 | |
---|
| 145 | allocate(this%lw_up_band(config%n_spec_lw,istartcol:iendcol,nlev+1)) |
---|
| 146 | allocate(this%lw_dn_band(config%n_spec_lw,istartcol:iendcol,nlev+1)) |
---|
| 147 | if (config%do_clear) then |
---|
| 148 | allocate(this%lw_up_clear_band(config%n_spec_lw, & |
---|
| 149 | & istartcol:iendcol,nlev+1)) |
---|
| 150 | allocate(this%lw_dn_clear_band(config%n_spec_lw, & |
---|
| 151 | & istartcol:iendcol,nlev+1)) |
---|
| 152 | end if |
---|
| 153 | end if |
---|
| 154 | |
---|
| 155 | if (config%do_lw_derivatives) then |
---|
| 156 | allocate(this%lw_derivatives(istartcol:iendcol,nlev+1)) |
---|
| 157 | end if |
---|
| 158 | |
---|
| 159 | ! Allocate g-point downwelling fluxes at surface passed from |
---|
| 160 | ! solver to surface_intermediate%partition |
---|
| 161 | allocate(this%lw_dn_surf_g(config%n_g_lw,istartcol:iendcol)) |
---|
| 162 | if (config%do_clear) then |
---|
| 163 | allocate(this%lw_dn_surf_clear_g(config%n_g_lw,istartcol:iendcol)) |
---|
| 164 | end if |
---|
| 165 | |
---|
| 166 | if (config%do_canopy_fluxes_lw) then |
---|
| 167 | ! Downward fluxes at top of canopy at the spectral resolution |
---|
| 168 | ! used in the canopy radiative transfer scheme |
---|
| 169 | allocate(this%lw_dn_surf_canopy(config%n_canopy_bands_lw,istartcol:iendcol)) |
---|
| 170 | end if |
---|
| 171 | end if |
---|
| 172 | |
---|
| 173 | ! Allocate shortwave arrays |
---|
| 174 | if (config%do_sw) then |
---|
| 175 | allocate(this%sw_up(istartcol:iendcol,nlev+1)) |
---|
| 176 | allocate(this%sw_dn(istartcol:iendcol,nlev+1)) |
---|
| 177 | if (config%do_sw_direct) then |
---|
| 178 | allocate(this%sw_dn_direct(istartcol:iendcol,nlev+1)) |
---|
| 179 | end if |
---|
| 180 | if (config%do_clear) then |
---|
| 181 | allocate(this%sw_up_clear(istartcol:iendcol,nlev+1)) |
---|
| 182 | allocate(this%sw_dn_clear(istartcol:iendcol,nlev+1)) |
---|
| 183 | if (config%do_sw_direct) then |
---|
| 184 | allocate(this%sw_dn_direct_clear(istartcol:iendcol,nlev+1)) |
---|
| 185 | end if |
---|
| 186 | end if |
---|
| 187 | |
---|
| 188 | if (config%do_save_spectral_flux) then |
---|
| 189 | if (config%n_spec_sw == 0) then |
---|
| 190 | write(nulerr,'(a)') '*** Error: number of SW spectral points to save not yet defined ' & |
---|
| 191 | & // 'so cannot allocate spectral flux arrays' |
---|
| 192 | call radiation_abort() |
---|
| 193 | end if |
---|
| 194 | |
---|
| 195 | allocate(this%sw_up_band(config%n_spec_sw,istartcol:iendcol,nlev+1)) |
---|
| 196 | allocate(this%sw_dn_band(config%n_spec_sw,istartcol:iendcol,nlev+1)) |
---|
| 197 | |
---|
| 198 | if (config%do_sw_direct) then |
---|
| 199 | allocate(this%sw_dn_direct_band(config%n_spec_sw, & |
---|
| 200 | & istartcol:iendcol,nlev+1)) |
---|
| 201 | end if |
---|
| 202 | if (config%do_clear) then |
---|
| 203 | allocate(this%sw_up_clear_band(config%n_spec_sw, & |
---|
| 204 | & istartcol:iendcol,nlev+1)) |
---|
| 205 | allocate(this%sw_dn_clear_band(config%n_spec_sw, & |
---|
| 206 | & istartcol:iendcol,nlev+1)) |
---|
| 207 | if (config%do_sw_direct) then |
---|
| 208 | allocate(this%sw_dn_direct_clear_band(config%n_spec_sw, & |
---|
| 209 | & istartcol:iendcol, nlev+1)) |
---|
| 210 | end if |
---|
| 211 | end if |
---|
| 212 | end if |
---|
| 213 | |
---|
| 214 | if (config%do_surface_sw_spectral_flux) then |
---|
| 215 | if (config%n_bands_sw == 0) then |
---|
| 216 | write(nulerr,'(a)') '*** Error: number of SW bands not yet defined ' & |
---|
| 217 | & // 'so cannot allocate surface spectral flux arrays' |
---|
| 218 | call radiation_abort() |
---|
| 219 | end if |
---|
| 220 | allocate(this%sw_dn_surf_band(config%n_bands_sw,istartcol:iendcol)) |
---|
| 221 | allocate(this%sw_dn_direct_surf_band(config%n_bands_sw,istartcol:iendcol)) |
---|
| 222 | if (config%do_clear) then |
---|
| 223 | allocate(this%sw_dn_surf_clear_band(config%n_bands_sw, & |
---|
| 224 | & istartcol:iendcol)) |
---|
| 225 | allocate(this%sw_dn_direct_surf_clear_band(config%n_bands_sw, & |
---|
| 226 | & istartcol:iendcol)) |
---|
| 227 | end if |
---|
| 228 | end if |
---|
| 229 | |
---|
| 230 | ! Allocate g-point downwelling fluxes at surface passed from |
---|
| 231 | ! solver to surface_intermediate%partition |
---|
| 232 | allocate(this%sw_dn_diffuse_surf_g(config%n_g_sw,istartcol:iendcol)) |
---|
| 233 | allocate(this%sw_dn_direct_surf_g (config%n_g_sw,istartcol:iendcol)) |
---|
| 234 | if (config%do_clear) then |
---|
| 235 | allocate(this%sw_dn_diffuse_surf_clear_g(config%n_g_sw,istartcol:iendcol)) |
---|
| 236 | allocate(this%sw_dn_direct_surf_clear_g (config%n_g_sw,istartcol:iendcol)) |
---|
| 237 | end if |
---|
| 238 | |
---|
| 239 | if (config%do_canopy_fluxes_sw) then |
---|
| 240 | ! Downward fluxes at top of canopy at the spectral resolution |
---|
| 241 | ! used in the canopy radiative transfer scheme |
---|
| 242 | allocate(this%sw_dn_diffuse_surf_canopy(config%n_canopy_bands_sw,istartcol:iendcol)) |
---|
| 243 | allocate(this%sw_dn_direct_surf_canopy (config%n_canopy_bands_sw,istartcol:iendcol)) |
---|
| 244 | end if |
---|
| 245 | end if |
---|
| 246 | |
---|
| 247 | ! Allocate cloud cover arrays |
---|
| 248 | allocate(this%cloud_cover_lw(istartcol:iendcol)) |
---|
| 249 | allocate(this%cloud_cover_sw(istartcol:iendcol)) |
---|
| 250 | |
---|
| 251 | ! Some solvers may not write to cloud cover, so we initialize to |
---|
| 252 | ! an unphysical value |
---|
| 253 | this%cloud_cover_lw = -1.0_jprb |
---|
| 254 | this%cloud_cover_sw = -1.0_jprb |
---|
| 255 | |
---|
| 256 | if (lhook) call dr_hook('radiation_flux:allocate',1,hook_handle) |
---|
| 257 | |
---|
| 258 | end subroutine allocate_flux_type |
---|
| 259 | |
---|
| 260 | |
---|
| 261 | !--------------------------------------------------------------------- |
---|
| 262 | ! Deallocate flux arrays |
---|
| 263 | subroutine deallocate_flux_type(this) |
---|
| 264 | |
---|
| 265 | use yomhook, only : lhook, dr_hook |
---|
| 266 | |
---|
| 267 | class(flux_type), intent(inout) :: this |
---|
| 268 | real(jprb) :: hook_handle |
---|
| 269 | |
---|
| 270 | if (lhook) call dr_hook('radiation_flux:deallocate',0,hook_handle) |
---|
| 271 | |
---|
| 272 | if (allocated(this%lw_up)) then |
---|
| 273 | deallocate(this%lw_up) |
---|
| 274 | if (allocated(this%lw_dn)) deallocate(this%lw_dn) |
---|
| 275 | if (allocated(this%lw_up_clear)) deallocate(this%lw_up_clear) |
---|
| 276 | if (allocated(this%lw_dn_clear)) deallocate(this%lw_dn_clear) |
---|
| 277 | end if |
---|
| 278 | |
---|
| 279 | if (allocated(this%sw_up)) then |
---|
| 280 | deallocate(this%sw_up) |
---|
| 281 | if (allocated(this%sw_dn)) deallocate(this%sw_dn) |
---|
| 282 | if (allocated(this%sw_up_clear)) deallocate(this%sw_up_clear) |
---|
| 283 | if (allocated(this%sw_dn_clear)) deallocate(this%sw_dn_clear) |
---|
| 284 | if (allocated(this%sw_dn_direct)) deallocate(this%sw_dn_direct) |
---|
| 285 | if (allocated(this%sw_dn_direct_clear)) & |
---|
| 286 | & deallocate(this%sw_dn_direct_clear) |
---|
| 287 | end if |
---|
| 288 | |
---|
| 289 | if (allocated(this%lw_up_band)) then |
---|
| 290 | deallocate(this%lw_up_band) |
---|
| 291 | if (allocated(this%lw_dn_band)) deallocate(this%lw_dn_band) |
---|
| 292 | if (allocated(this%lw_up_clear_band)) deallocate(this%lw_up_clear_band) |
---|
| 293 | if (allocated(this%lw_dn_clear_band)) deallocate(this%lw_dn_clear_band) |
---|
| 294 | end if |
---|
| 295 | |
---|
| 296 | if (allocated(this%sw_up_band)) then |
---|
| 297 | deallocate(this%sw_up_band) |
---|
| 298 | if (allocated(this%sw_dn_band)) deallocate(this%sw_dn_band) |
---|
| 299 | if (allocated(this%sw_up_clear_band)) deallocate(this%sw_up_clear_band) |
---|
| 300 | if (allocated(this%sw_dn_clear_band)) deallocate(this%sw_dn_clear_band) |
---|
| 301 | if (allocated(this%sw_dn_direct_band)) deallocate(this%sw_dn_direct_band) |
---|
| 302 | if (allocated(this%sw_dn_direct_clear_band)) & |
---|
| 303 | & deallocate(this%sw_dn_direct_clear_band) |
---|
| 304 | end if |
---|
| 305 | |
---|
| 306 | if (allocated(this%sw_dn_surf_band)) then |
---|
| 307 | deallocate(this%sw_dn_surf_band) |
---|
| 308 | deallocate(this%sw_dn_direct_surf_band) |
---|
| 309 | end if |
---|
| 310 | if (allocated(this%sw_dn_surf_clear_band)) then |
---|
| 311 | deallocate(this%sw_dn_surf_clear_band) |
---|
| 312 | deallocate(this%sw_dn_direct_surf_clear_band) |
---|
| 313 | end if |
---|
| 314 | |
---|
| 315 | if (allocated(this%lw_dn_surf_canopy)) deallocate(this%lw_dn_surf_canopy) |
---|
| 316 | if (allocated(this%sw_dn_diffuse_surf_canopy)) deallocate(this%sw_dn_diffuse_surf_canopy) |
---|
| 317 | if (allocated(this%sw_dn_direct_surf_canopy)) deallocate(this%sw_dn_direct_surf_canopy) |
---|
| 318 | |
---|
| 319 | if (allocated(this%cloud_cover_sw)) then |
---|
| 320 | deallocate(this%cloud_cover_sw) |
---|
| 321 | end if |
---|
| 322 | if (allocated(this%cloud_cover_lw)) then |
---|
| 323 | deallocate(this%cloud_cover_lw) |
---|
| 324 | end if |
---|
| 325 | |
---|
| 326 | if (allocated(this%lw_derivatives)) then |
---|
| 327 | deallocate(this%lw_derivatives) |
---|
| 328 | end if |
---|
| 329 | |
---|
[4584] | 330 | if (allocated(this%lw_dn_surf_g)) deallocate(this%lw_dn_surf_g) |
---|
| 331 | if (allocated(this%lw_dn_surf_clear_g)) deallocate(this%lw_dn_surf_clear_g) |
---|
| 332 | if (allocated(this%sw_dn_diffuse_surf_g)) deallocate(this%sw_dn_diffuse_surf_g) |
---|
| 333 | if (allocated(this%sw_dn_direct_surf_g)) deallocate(this%sw_dn_direct_surf_g) |
---|
| 334 | if (allocated(this%sw_dn_diffuse_surf_clear_g)) deallocate(this%sw_dn_diffuse_surf_clear_g) |
---|
| 335 | if (allocated(this%sw_dn_direct_surf_clear_g)) deallocate(this%sw_dn_direct_surf_clear_g) |
---|
| 336 | |
---|
[3908] | 337 | if (lhook) call dr_hook('radiation_flux:deallocate',1,hook_handle) |
---|
| 338 | |
---|
| 339 | end subroutine deallocate_flux_type |
---|
| 340 | |
---|
| 341 | !--------------------------------------------------------------------- |
---|
| 342 | ! Calculate surface downwelling fluxes in each band using the |
---|
| 343 | ! downwelling surface fluxes at each g point |
---|
| 344 | subroutine calc_surface_spectral(this, config, istartcol, iendcol) |
---|
| 345 | |
---|
| 346 | use yomhook, only : lhook, dr_hook |
---|
| 347 | use radiation_config, only : config_type |
---|
| 348 | |
---|
| 349 | class(flux_type), intent(inout) :: this |
---|
| 350 | type(config_type), intent(in) :: config |
---|
| 351 | integer, intent(in) :: istartcol, iendcol |
---|
| 352 | |
---|
| 353 | integer :: jcol, jband, jalbedoband, nalbedoband |
---|
| 354 | |
---|
| 355 | ! Longwave surface downwelling in each band needed to compute |
---|
| 356 | ! canopy fluxes |
---|
| 357 | real(jprb) :: lw_dn_surf_band(config%n_bands_lw,istartcol:iendcol) |
---|
| 358 | |
---|
| 359 | real(jprb) :: hook_handle |
---|
| 360 | |
---|
| 361 | if (lhook) call dr_hook('radiation_flux:calc_surface_spectral',0,hook_handle) |
---|
| 362 | |
---|
| 363 | if (config%do_sw .and. config%do_surface_sw_spectral_flux) then |
---|
| 364 | |
---|
[4584] | 365 | if (use_indexed_sum_vec) then |
---|
| 366 | call indexed_sum_vec(this%sw_dn_direct_surf_g, & |
---|
| 367 | & config%i_band_from_reordered_g_sw, & |
---|
| 368 | & this%sw_dn_direct_surf_band, istartcol, iendcol) |
---|
| 369 | call indexed_sum_vec(this%sw_dn_diffuse_surf_g, & |
---|
| 370 | & config%i_band_from_reordered_g_sw, & |
---|
| 371 | & this%sw_dn_surf_band, istartcol, iendcol) |
---|
[3908] | 372 | do jcol = istartcol,iendcol |
---|
[4584] | 373 | this%sw_dn_surf_band(:,jcol) & |
---|
| 374 | & = this%sw_dn_surf_band(:,jcol) & |
---|
| 375 | & + this%sw_dn_direct_surf_band(:,jcol) |
---|
| 376 | end do |
---|
| 377 | else |
---|
| 378 | do jcol = istartcol,iendcol |
---|
| 379 | call indexed_sum(this%sw_dn_direct_surf_g(:,jcol), & |
---|
[3908] | 380 | & config%i_band_from_reordered_g_sw, & |
---|
[4584] | 381 | & this%sw_dn_direct_surf_band(:,jcol)) |
---|
| 382 | call indexed_sum(this%sw_dn_diffuse_surf_g(:,jcol), & |
---|
[3908] | 383 | & config%i_band_from_reordered_g_sw, & |
---|
[4584] | 384 | & this%sw_dn_surf_band(:,jcol)) |
---|
| 385 | this%sw_dn_surf_band(:,jcol) & |
---|
| 386 | & = this%sw_dn_surf_band(:,jcol) & |
---|
| 387 | & + this%sw_dn_direct_surf_band(:,jcol) |
---|
[3908] | 388 | end do |
---|
| 389 | end if |
---|
| 390 | |
---|
[4584] | 391 | if (config%do_clear) then |
---|
| 392 | if (use_indexed_sum_vec) then |
---|
| 393 | call indexed_sum_vec(this%sw_dn_direct_surf_clear_g, & |
---|
| 394 | & config%i_band_from_reordered_g_sw, & |
---|
| 395 | & this%sw_dn_direct_surf_clear_band, istartcol, iendcol) |
---|
| 396 | call indexed_sum_vec(this%sw_dn_diffuse_surf_clear_g, & |
---|
| 397 | & config%i_band_from_reordered_g_sw, & |
---|
| 398 | & this%sw_dn_surf_clear_band, istartcol, iendcol) |
---|
| 399 | do jcol = istartcol,iendcol |
---|
| 400 | this%sw_dn_surf_clear_band(:,jcol) & |
---|
| 401 | & = this%sw_dn_surf_clear_band(:,jcol) & |
---|
| 402 | & + this%sw_dn_direct_surf_clear_band(:,jcol) |
---|
| 403 | end do |
---|
| 404 | else |
---|
| 405 | do jcol = istartcol,iendcol |
---|
| 406 | call indexed_sum(this%sw_dn_direct_surf_clear_g(:,jcol), & |
---|
| 407 | & config%i_band_from_reordered_g_sw, & |
---|
| 408 | & this%sw_dn_direct_surf_clear_band(:,jcol)) |
---|
| 409 | call indexed_sum(this%sw_dn_diffuse_surf_clear_g(:,jcol), & |
---|
| 410 | & config%i_band_from_reordered_g_sw, & |
---|
| 411 | & this%sw_dn_surf_clear_band(:,jcol)) |
---|
| 412 | this%sw_dn_surf_clear_band(:,jcol) & |
---|
| 413 | & = this%sw_dn_surf_clear_band(:,jcol) & |
---|
| 414 | & + this%sw_dn_direct_surf_clear_band(:,jcol) |
---|
| 415 | end do |
---|
| 416 | end if |
---|
| 417 | end if |
---|
| 418 | |
---|
[3908] | 419 | end if ! do_surface_sw_spectral_flux |
---|
| 420 | |
---|
| 421 | ! Fluxes in bands required for canopy radiative transfer |
---|
| 422 | if (config%do_sw .and. config%do_canopy_fluxes_sw) then |
---|
| 423 | if (config%use_canopy_full_spectrum_sw) then |
---|
| 424 | this%sw_dn_diffuse_surf_canopy(:,istartcol:iendcol) = this%sw_dn_diffuse_surf_g(:,istartcol:iendcol) |
---|
| 425 | this%sw_dn_direct_surf_canopy (:,istartcol:iendcol) = this%sw_dn_direct_surf_g (:,istartcol:iendcol) |
---|
| 426 | else if (config%do_nearest_spectral_sw_albedo) then |
---|
[4584] | 427 | if (use_indexed_sum_vec) then |
---|
| 428 | call indexed_sum_vec(this%sw_dn_direct_surf_g, & |
---|
| 429 | & config%i_albedo_from_band_sw(config%i_band_from_reordered_g_sw), & |
---|
| 430 | & this%sw_dn_direct_surf_canopy, istartcol, iendcol) |
---|
| 431 | call indexed_sum_vec(this%sw_dn_diffuse_surf_g, & |
---|
| 432 | & config%i_albedo_from_band_sw(config%i_band_from_reordered_g_sw), & |
---|
| 433 | & this%sw_dn_diffuse_surf_canopy, istartcol, iendcol) |
---|
| 434 | else |
---|
| 435 | do jcol = istartcol,iendcol |
---|
| 436 | call indexed_sum(this%sw_dn_direct_surf_g(:,jcol), & |
---|
| 437 | & config%i_albedo_from_band_sw(config%i_band_from_reordered_g_sw), & |
---|
| 438 | & this%sw_dn_direct_surf_canopy(:,jcol)) |
---|
| 439 | call indexed_sum(this%sw_dn_diffuse_surf_g(:,jcol), & |
---|
| 440 | & config%i_albedo_from_band_sw(config%i_band_from_reordered_g_sw), & |
---|
| 441 | & this%sw_dn_diffuse_surf_canopy(:,jcol)) |
---|
| 442 | end do |
---|
| 443 | end if |
---|
[3908] | 444 | else |
---|
| 445 | ! More accurate calculations using weights, but requires |
---|
| 446 | ! this%sw_dn_[direct_]surf_band to be defined, i.e. |
---|
| 447 | ! config%do_surface_sw_spectral_flux == .true. |
---|
| 448 | nalbedoband = size(config%sw_albedo_weights,1) |
---|
| 449 | this%sw_dn_diffuse_surf_canopy(:,istartcol:iendcol) = 0.0_jprb |
---|
| 450 | this%sw_dn_direct_surf_canopy (:,istartcol:iendcol) = 0.0_jprb |
---|
| 451 | do jband = 1,config%n_bands_sw |
---|
| 452 | do jalbedoband = 1,nalbedoband |
---|
| 453 | if (config%sw_albedo_weights(jalbedoband,jband) /= 0.0_jprb) then |
---|
| 454 | ! Initially, "diffuse" is actually "total" |
---|
| 455 | this%sw_dn_diffuse_surf_canopy(jalbedoband,istartcol:iendcol) & |
---|
| 456 | & = this%sw_dn_diffuse_surf_canopy(jalbedoband,istartcol:iendcol) & |
---|
| 457 | & + config%sw_albedo_weights(jalbedoband,jband) & |
---|
| 458 | & * this%sw_dn_surf_band(jband,istartcol:iendcol) |
---|
| 459 | this%sw_dn_direct_surf_canopy(jalbedoband,istartcol:iendcol) & |
---|
| 460 | & = this%sw_dn_direct_surf_canopy(jalbedoband,istartcol:iendcol) & |
---|
| 461 | & + config%sw_albedo_weights(jalbedoband,jband) & |
---|
| 462 | & * this%sw_dn_direct_surf_band(jband,istartcol:iendcol) |
---|
| 463 | end if |
---|
| 464 | end do |
---|
| 465 | end do |
---|
| 466 | ! Subtract the direct from total to get diffuse |
---|
| 467 | this%sw_dn_diffuse_surf_canopy(:,istartcol:iendcol) & |
---|
| 468 | & = this%sw_dn_diffuse_surf_canopy(:,istartcol:iendcol) & |
---|
| 469 | & - this%sw_dn_direct_surf_canopy(:,istartcol:iendcol) |
---|
| 470 | end if |
---|
| 471 | |
---|
| 472 | end if ! do_canopy_fluxes_sw |
---|
| 473 | |
---|
| 474 | if (config%do_lw .and. config%do_canopy_fluxes_lw) then |
---|
| 475 | if (config%use_canopy_full_spectrum_lw) then |
---|
| 476 | this%lw_dn_surf_canopy(:,istartcol:iendcol) = this%lw_dn_surf_g(:,istartcol:iendcol) |
---|
| 477 | else if (config%do_nearest_spectral_lw_emiss) then |
---|
[4584] | 478 | if (use_indexed_sum_vec) then |
---|
| 479 | call indexed_sum_vec(this%lw_dn_surf_g, & |
---|
| 480 | & config%i_emiss_from_band_lw(config%i_band_from_reordered_g_lw), & |
---|
| 481 | & this%lw_dn_surf_canopy, istartcol, iendcol) |
---|
| 482 | else |
---|
| 483 | do jcol = istartcol,iendcol |
---|
| 484 | call indexed_sum(this%lw_dn_surf_g(:,jcol), & |
---|
| 485 | & config%i_emiss_from_band_lw(config%i_band_from_reordered_g_lw), & |
---|
| 486 | & this%lw_dn_surf_canopy(:,jcol)) |
---|
| 487 | end do |
---|
| 488 | end if |
---|
[3908] | 489 | else |
---|
| 490 | ! Compute fluxes in each longwave emissivity interval using |
---|
| 491 | ! weights; first sum over g points to get the values in bands |
---|
[4584] | 492 | if (use_indexed_sum_vec) then |
---|
| 493 | call indexed_sum_vec(this%lw_dn_surf_g, & |
---|
| 494 | & config%i_band_from_reordered_g_lw, & |
---|
| 495 | & lw_dn_surf_band, istartcol, iendcol) |
---|
| 496 | else |
---|
| 497 | do jcol = istartcol,iendcol |
---|
| 498 | call indexed_sum(this%lw_dn_surf_g(:,jcol), & |
---|
| 499 | & config%i_band_from_reordered_g_lw, & |
---|
| 500 | & lw_dn_surf_band(:,jcol)) |
---|
| 501 | end do |
---|
| 502 | end if |
---|
[3908] | 503 | nalbedoband = size(config%lw_emiss_weights,1) |
---|
| 504 | this%lw_dn_surf_canopy(:,istartcol:iendcol) = 0.0_jprb |
---|
| 505 | do jband = 1,config%n_bands_lw |
---|
| 506 | do jalbedoband = 1,nalbedoband |
---|
| 507 | if (config%lw_emiss_weights(jalbedoband,jband) /= 0.0_jprb) then |
---|
| 508 | this%lw_dn_surf_canopy(jalbedoband,istartcol:iendcol) & |
---|
| 509 | & = this%lw_dn_surf_canopy(jalbedoband,istartcol:iendcol) & |
---|
| 510 | & + config%lw_emiss_weights(jalbedoband,jband) & |
---|
| 511 | & * lw_dn_surf_band(jband,istartcol:iendcol) |
---|
| 512 | end if |
---|
| 513 | end do |
---|
| 514 | end do |
---|
| 515 | end if |
---|
| 516 | end if |
---|
| 517 | |
---|
| 518 | if (lhook) call dr_hook('radiation_flux:calc_surface_spectral',1,hook_handle) |
---|
| 519 | |
---|
| 520 | end subroutine calc_surface_spectral |
---|
| 521 | |
---|
| 522 | |
---|
| 523 | !--------------------------------------------------------------------- |
---|
| 524 | ! Return .true. if the most important flux variables are out of a |
---|
| 525 | ! physically sensible range, optionally only considering columns |
---|
| 526 | ! between istartcol and iendcol |
---|
| 527 | function out_of_physical_bounds(this, istartcol, iendcol) result(is_bad) |
---|
| 528 | |
---|
| 529 | use yomhook, only : lhook, dr_hook |
---|
[4584] | 530 | use radiation_check, only : out_of_bounds_2d |
---|
[3908] | 531 | |
---|
| 532 | class(flux_type), intent(inout) :: this |
---|
| 533 | integer, optional,intent(in) :: istartcol, iendcol |
---|
| 534 | logical :: is_bad |
---|
| 535 | |
---|
| 536 | real(jprb) :: hook_handle |
---|
| 537 | |
---|
| 538 | if (lhook) call dr_hook('radiation_flux:out_of_physical_bounds',0,hook_handle) |
---|
| 539 | |
---|
| 540 | is_bad = out_of_bounds_2d(this%lw_up, 'lw_up', 10.0_jprb, 900.0_jprb, .false., i1=istartcol, i2=iendcol) & |
---|
| 541 | & .or. out_of_bounds_2d(this%lw_dn, 'lw_dn', 0.0_jprb, 800.0_jprb, .false., i1=istartcol, i2=iendcol) & |
---|
| 542 | & .or. out_of_bounds_2d(this%sw_up, 'sw_up', 0.0_jprb, 1500.0_jprb, .false., i1=istartcol, i2=iendcol) & |
---|
| 543 | & .or. out_of_bounds_2d(this%sw_dn, 'sw_dn', 0.0_jprb, 1500.0_jprb, .false., i1=istartcol, i2=iendcol) & |
---|
| 544 | & .or. out_of_bounds_2d(this%sw_dn_direct, 'sw_dn_direct', 0.0_jprb, 1500.0_jprb, .false., i1=istartcol, i2=iendcol) & |
---|
| 545 | & .or. out_of_bounds_2d(this%lw_derivatives, 'lw_derivatives', 0.0_jprb, 1.0_jprb, .false., i1=istartcol, i2=iendcol) & |
---|
| 546 | & .or. out_of_bounds_2d(this%sw_dn_surf_band, 'sw_dn_surf_band', 0.0_jprb, 1500.0_jprb, & |
---|
| 547 | & .false., j1=istartcol, j2=iendcol) & |
---|
| 548 | & .or. out_of_bounds_2d(this%sw_dn_surf_clear_band, 'sw_dn_surf_clear_band', 0.0_jprb, 1500.0_jprb, & |
---|
| 549 | & .false., j1=istartcol, j2=iendcol) |
---|
| 550 | |
---|
| 551 | if (lhook) call dr_hook('radiation_flux:out_of_physical_bounds',1,hook_handle) |
---|
| 552 | |
---|
| 553 | end function out_of_physical_bounds |
---|
| 554 | |
---|
| 555 | !--------------------------------------------------------------------- |
---|
| 556 | ! Return .true. if the heating rates are out of a physically |
---|
| 557 | ! sensible range, optionally only considering columns between |
---|
| 558 | ! istartcol and iendcol. This function allocates and deallocates |
---|
| 559 | ! memory due to the requirements for inputs of out_of_bounds_2d. |
---|
| 560 | function heating_rate_out_of_physical_bounds(this, nlev, istartcol, iendcol, pressure_hl) result(is_bad) |
---|
| 561 | |
---|
[4584] | 562 | use radiation_check, only : out_of_bounds_2d |
---|
[3908] | 563 | use radiation_constants, only : AccelDueToGravity |
---|
| 564 | |
---|
| 565 | ! "Cp" (J kg-1 K-1) |
---|
| 566 | real(jprb), parameter :: SpecificHeatDryAir = 1004.0 |
---|
| 567 | |
---|
| 568 | class(flux_type), intent(inout) :: this |
---|
| 569 | integer, intent(in) :: istartcol, iendcol, nlev |
---|
| 570 | logical :: is_bad |
---|
| 571 | |
---|
| 572 | real(jprb), intent(in) :: pressure_hl(:,:) |
---|
| 573 | |
---|
| 574 | real(jprb), allocatable :: hr_K_day(:,:) |
---|
| 575 | |
---|
| 576 | real(jprb) :: scaling(istartcol:iendcol,nlev) |
---|
| 577 | |
---|
| 578 | allocate(hr_K_day(istartcol:iendcol,nlev)) |
---|
| 579 | |
---|
| 580 | scaling = -(24.0_jprb * 3600.0_jprb * AccelDueToGravity / SpecificHeatDryAir) & |
---|
| 581 | & / (pressure_hl(istartcol:iendcol,2:nlev+1) - pressure_hl(istartcol:iendcol,1:nlev)) |
---|
| 582 | ! Shortwave |
---|
| 583 | hr_K_day = scaling * (this%sw_dn(istartcol:iendcol,2:nlev+1) - this%sw_up(istartcol:iendcol,2:nlev+1) & |
---|
| 584 | & -this%sw_dn(istartcol:iendcol,1:nlev) + this%sw_up(istartcol:iendcol,1:nlev)) |
---|
| 585 | is_bad = out_of_bounds_2d(hr_K_day, 'sw_heating_rate_K_day', 0.0_jprb, 200.0_jprb, & |
---|
| 586 | & .false., i1=istartcol, i2=iendcol) |
---|
| 587 | |
---|
| 588 | ! Longwave |
---|
| 589 | hr_K_day = scaling * (this%lw_dn(istartcol:iendcol,2:nlev+1) - this%lw_up(istartcol:iendcol,2:nlev+1) & |
---|
| 590 | & -this%lw_dn(istartcol:iendcol,1:nlev) + this%lw_up(istartcol:iendcol,1:nlev)) |
---|
| 591 | is_bad = is_bad .or. out_of_bounds_2d(hr_K_day, 'lw_heating_rate_K_day', -250.0_jprb, 150.0_jprb, & |
---|
| 592 | & .false., i1=istartcol, i2=iendcol) |
---|
| 593 | |
---|
| 594 | deallocate(hr_K_day) |
---|
| 595 | |
---|
| 596 | end function heating_rate_out_of_physical_bounds |
---|
| 597 | |
---|
| 598 | |
---|
| 599 | !--------------------------------------------------------------------- |
---|
| 600 | ! Sum elements of "source" into "dest" according to index "ind". |
---|
| 601 | ! "source" and "ind" should have the same size and bounds, and no |
---|
| 602 | ! element of "ind" should refer outside the bounds of "dest". This |
---|
| 603 | ! version increments existing contents of "dest". |
---|
| 604 | pure subroutine add_indexed_sum(source, ind, dest) |
---|
| 605 | |
---|
| 606 | real(jprb), intent(in) :: source(:) |
---|
| 607 | integer, intent(in) :: ind(:) |
---|
| 608 | real(jprb), intent(inout) :: dest(:) |
---|
| 609 | |
---|
| 610 | integer :: ig, jg, istart, iend |
---|
| 611 | |
---|
| 612 | istart = lbound(source,1) |
---|
| 613 | iend = ubound(source,1) |
---|
| 614 | |
---|
| 615 | do jg = istart, iend |
---|
| 616 | ig = ind(jg) |
---|
| 617 | dest(ig) = dest(ig) + source(jg) |
---|
| 618 | end do |
---|
| 619 | |
---|
| 620 | end subroutine add_indexed_sum |
---|
| 621 | |
---|
| 622 | |
---|
| 623 | !--------------------------------------------------------------------- |
---|
| 624 | ! As "add_indexed_sum" but this version overwrites existing contents |
---|
| 625 | ! of "dest" |
---|
| 626 | pure subroutine indexed_sum(source, ind, dest) |
---|
| 627 | |
---|
| 628 | real(jprb), intent(in) :: source(:) |
---|
| 629 | integer, intent(in) :: ind(:) |
---|
| 630 | real(jprb), intent(out) :: dest(:) |
---|
| 631 | |
---|
| 632 | integer :: ig, jg, istart, iend |
---|
| 633 | |
---|
| 634 | dest = 0.0 |
---|
| 635 | |
---|
| 636 | istart = lbound(source,1) |
---|
| 637 | iend = ubound(source,1) |
---|
| 638 | |
---|
| 639 | do jg = istart, iend |
---|
| 640 | ig = ind(jg) |
---|
| 641 | dest(ig) = dest(ig) + source(jg) |
---|
| 642 | end do |
---|
| 643 | |
---|
| 644 | end subroutine indexed_sum |
---|
| 645 | |
---|
[4584] | 646 | !--------------------------------------------------------------------- |
---|
| 647 | ! Vectorized version of "add_indexed_sum" |
---|
| 648 | subroutine indexed_sum_vec(source, ind, dest, ist, iend) |
---|
[3908] | 649 | |
---|
[4584] | 650 | real(jprb), intent(in) :: source(:,:) |
---|
| 651 | integer, intent(in) :: ind(:) |
---|
| 652 | real(jprb), intent(out) :: dest(:,:) |
---|
| 653 | integer, intent(in) :: ist, iend |
---|
| 654 | |
---|
| 655 | integer :: ig, jg, jc |
---|
| 656 | |
---|
| 657 | dest = 0.0 |
---|
| 658 | |
---|
| 659 | do jg = lbound(source,1), ubound(source,1) |
---|
| 660 | ig = ind(jg) |
---|
| 661 | do jc = ist, iend |
---|
| 662 | dest(ig,jc) = dest(ig,jc) + source(jg,jc) |
---|
| 663 | end do |
---|
| 664 | end do |
---|
| 665 | |
---|
| 666 | end subroutine indexed_sum_vec |
---|
| 667 | |
---|
[3908] | 668 | !--------------------------------------------------------------------- |
---|
| 669 | ! As "add_indexed_sum" but a whole vertical profiles |
---|
| 670 | pure subroutine add_indexed_sum_profile(source, ind, dest) |
---|
| 671 | |
---|
| 672 | real(jprb), intent(in) :: source(:,:) |
---|
| 673 | integer, intent(in) :: ind(:) |
---|
| 674 | real(jprb), intent(out) :: dest(:,:) |
---|
| 675 | |
---|
| 676 | integer :: ig, jg, istart, iend, jlev, nlev |
---|
| 677 | |
---|
| 678 | istart = lbound(source,1) |
---|
| 679 | iend = ubound(source,1) |
---|
| 680 | nlev = size(source,2) |
---|
| 681 | |
---|
| 682 | do jlev = 1,nlev |
---|
| 683 | do jg = istart, iend |
---|
| 684 | ig = ind(jg) |
---|
| 685 | dest(ig,jlev) = dest(ig,jlev) + source(jg,jlev) |
---|
| 686 | end do |
---|
| 687 | end do |
---|
| 688 | |
---|
| 689 | end subroutine add_indexed_sum_profile |
---|
| 690 | |
---|
| 691 | |
---|
| 692 | !--------------------------------------------------------------------- |
---|
| 693 | ! As "indexed_sum" but a whole vertical profiles |
---|
| 694 | pure subroutine indexed_sum_profile(source, ind, dest) |
---|
| 695 | |
---|
| 696 | real(jprb), intent(in) :: source(:,:) |
---|
| 697 | integer, intent(in) :: ind(:) |
---|
| 698 | real(jprb), intent(out) :: dest(:,:) |
---|
| 699 | |
---|
| 700 | integer :: ig, jg, istart, iend, jlev, nlev |
---|
| 701 | |
---|
| 702 | dest = 0.0 |
---|
| 703 | |
---|
| 704 | istart = lbound(source,1) |
---|
| 705 | iend = ubound(source,1) |
---|
| 706 | nlev = size(source,2) |
---|
| 707 | |
---|
| 708 | do jlev = 1,nlev |
---|
| 709 | do jg = istart, iend |
---|
| 710 | ig = ind(jg) |
---|
| 711 | dest(ig,jlev) = dest(ig,jlev) + source(jg,jlev) |
---|
| 712 | end do |
---|
| 713 | end do |
---|
| 714 | |
---|
| 715 | end subroutine indexed_sum_profile |
---|
| 716 | |
---|
| 717 | end module radiation_flux |
---|