[3880] | 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 | |
---|
| 102 | contains |
---|
| 103 | |
---|
| 104 | !--------------------------------------------------------------------- |
---|
| 105 | ! Allocate arrays for flux profiles, using config to define which |
---|
| 106 | ! fluxes are needed. The arrays are dimensioned for columns between |
---|
| 107 | ! istartcol, iendcol and levels from 1 to nlev+1 |
---|
| 108 | subroutine allocate_flux_type(this, config, istartcol, iendcol, nlev) |
---|
| 109 | |
---|
| 110 | use yomhook, only : lhook, dr_hook |
---|
| 111 | use radiation_io, only : nulerr, radiation_abort |
---|
| 112 | use radiation_config, only : config_type |
---|
| 113 | |
---|
| 114 | integer, intent(in) :: istartcol, iendcol, nlev |
---|
| 115 | class(flux_type), intent(inout) :: this |
---|
| 116 | type(config_type), intent(in) :: config |
---|
| 117 | |
---|
| 118 | real(jprb) :: hook_handle |
---|
| 119 | |
---|
| 120 | if (lhook) call dr_hook('radiation_flux:allocate',0,hook_handle) |
---|
| 121 | |
---|
| 122 | ! Allocate longwave arrays |
---|
| 123 | if (config%do_lw) then |
---|
| 124 | allocate(this%lw_up(istartcol:iendcol,nlev+1)) |
---|
| 125 | allocate(this%lw_dn(istartcol:iendcol,nlev+1)) |
---|
| 126 | if (config%do_clear) then |
---|
| 127 | allocate(this%lw_up_clear(istartcol:iendcol,nlev+1)) |
---|
| 128 | allocate(this%lw_dn_clear(istartcol:iendcol,nlev+1)) |
---|
| 129 | end if |
---|
| 130 | |
---|
| 131 | if (config%do_save_spectral_flux) then |
---|
| 132 | if (config%n_spec_lw == 0) then |
---|
| 133 | write(nulerr,'(a)') '*** Error: number of LW spectral points to save not yet defined ' & |
---|
| 134 | & // 'so cannot allocated spectral flux arrays' |
---|
| 135 | call radiation_abort() |
---|
| 136 | end if |
---|
| 137 | |
---|
| 138 | allocate(this%lw_up_band(config%n_spec_lw,istartcol:iendcol,nlev+1)) |
---|
| 139 | allocate(this%lw_dn_band(config%n_spec_lw,istartcol:iendcol,nlev+1)) |
---|
| 140 | if (config%do_clear) then |
---|
| 141 | allocate(this%lw_up_clear_band(config%n_spec_lw, & |
---|
| 142 | & istartcol:iendcol,nlev+1)) |
---|
| 143 | allocate(this%lw_dn_clear_band(config%n_spec_lw, & |
---|
| 144 | & istartcol:iendcol,nlev+1)) |
---|
| 145 | end if |
---|
| 146 | end if |
---|
| 147 | |
---|
| 148 | if (config%do_lw_derivatives) then |
---|
| 149 | allocate(this%lw_derivatives(istartcol:iendcol,nlev+1)) |
---|
| 150 | end if |
---|
| 151 | |
---|
| 152 | ! Allocate g-point downwelling fluxes at surface passed from |
---|
| 153 | ! solver to surface_intermediate%partition |
---|
| 154 | allocate(this%lw_dn_surf_g(config%n_g_lw,istartcol:iendcol)) |
---|
| 155 | if (config%do_clear) then |
---|
| 156 | allocate(this%lw_dn_surf_clear_g(config%n_g_lw,istartcol:iendcol)) |
---|
| 157 | end if |
---|
| 158 | |
---|
| 159 | if (config%do_canopy_fluxes_lw) then |
---|
| 160 | ! Downward fluxes at top of canopy at the spectral resolution |
---|
| 161 | ! used in the canopy radiative transfer scheme |
---|
| 162 | allocate(this%lw_dn_surf_canopy(config%n_canopy_bands_lw,istartcol:iendcol)) |
---|
| 163 | end if |
---|
| 164 | end if |
---|
| 165 | |
---|
| 166 | ! Allocate shortwave arrays |
---|
| 167 | if (config%do_sw) then |
---|
| 168 | allocate(this%sw_up(istartcol:iendcol,nlev+1)) |
---|
| 169 | allocate(this%sw_dn(istartcol:iendcol,nlev+1)) |
---|
| 170 | if (config%do_sw_direct) then |
---|
| 171 | allocate(this%sw_dn_direct(istartcol:iendcol,nlev+1)) |
---|
| 172 | end if |
---|
| 173 | if (config%do_clear) then |
---|
| 174 | allocate(this%sw_up_clear(istartcol:iendcol,nlev+1)) |
---|
| 175 | allocate(this%sw_dn_clear(istartcol:iendcol,nlev+1)) |
---|
| 176 | if (config%do_sw_direct) then |
---|
| 177 | allocate(this%sw_dn_direct_clear(istartcol:iendcol,nlev+1)) |
---|
| 178 | end if |
---|
| 179 | end if |
---|
| 180 | |
---|
| 181 | if (config%do_save_spectral_flux) then |
---|
| 182 | if (config%n_spec_sw == 0) then |
---|
| 183 | write(nulerr,'(a)') '*** Error: number of SW spectral points to save not yet defined ' & |
---|
| 184 | & // 'so cannot allocate spectral flux arrays' |
---|
| 185 | call radiation_abort() |
---|
| 186 | end if |
---|
| 187 | |
---|
| 188 | allocate(this%sw_up_band(config%n_spec_sw,istartcol:iendcol,nlev+1)) |
---|
| 189 | allocate(this%sw_dn_band(config%n_spec_sw,istartcol:iendcol,nlev+1)) |
---|
| 190 | |
---|
| 191 | if (config%do_sw_direct) then |
---|
| 192 | allocate(this%sw_dn_direct_band(config%n_spec_sw, & |
---|
| 193 | & istartcol:iendcol,nlev+1)) |
---|
| 194 | end if |
---|
| 195 | if (config%do_clear) then |
---|
| 196 | allocate(this%sw_up_clear_band(config%n_spec_sw, & |
---|
| 197 | & istartcol:iendcol,nlev+1)) |
---|
| 198 | allocate(this%sw_dn_clear_band(config%n_spec_sw, & |
---|
| 199 | & istartcol:iendcol,nlev+1)) |
---|
| 200 | if (config%do_sw_direct) then |
---|
| 201 | allocate(this%sw_dn_direct_clear_band(config%n_spec_sw, & |
---|
| 202 | & istartcol:iendcol, nlev+1)) |
---|
| 203 | end if |
---|
| 204 | end if |
---|
| 205 | end if |
---|
| 206 | |
---|
| 207 | if (config%do_surface_sw_spectral_flux) then |
---|
| 208 | if (config%n_bands_sw == 0) then |
---|
| 209 | write(nulerr,'(a)') '*** Error: number of SW bands not yet defined ' & |
---|
| 210 | & // 'so cannot allocate surface spectral flux arrays' |
---|
| 211 | call radiation_abort() |
---|
| 212 | end if |
---|
| 213 | allocate(this%sw_dn_surf_band(config%n_bands_sw,istartcol:iendcol)) |
---|
| 214 | allocate(this%sw_dn_direct_surf_band(config%n_bands_sw,istartcol:iendcol)) |
---|
| 215 | if (config%do_clear) then |
---|
| 216 | allocate(this%sw_dn_surf_clear_band(config%n_bands_sw, & |
---|
| 217 | & istartcol:iendcol)) |
---|
| 218 | allocate(this%sw_dn_direct_surf_clear_band(config%n_bands_sw, & |
---|
| 219 | & istartcol:iendcol)) |
---|
| 220 | end if |
---|
| 221 | end if |
---|
| 222 | |
---|
| 223 | ! Allocate g-point downwelling fluxes at surface passed from |
---|
| 224 | ! solver to surface_intermediate%partition |
---|
| 225 | allocate(this%sw_dn_diffuse_surf_g(config%n_g_sw,istartcol:iendcol)) |
---|
| 226 | allocate(this%sw_dn_direct_surf_g (config%n_g_sw,istartcol:iendcol)) |
---|
| 227 | if (config%do_clear) then |
---|
| 228 | allocate(this%sw_dn_diffuse_surf_clear_g(config%n_g_sw,istartcol:iendcol)) |
---|
| 229 | allocate(this%sw_dn_direct_surf_clear_g (config%n_g_sw,istartcol:iendcol)) |
---|
| 230 | end if |
---|
| 231 | |
---|
| 232 | if (config%do_canopy_fluxes_sw) then |
---|
| 233 | ! Downward fluxes at top of canopy at the spectral resolution |
---|
| 234 | ! used in the canopy radiative transfer scheme |
---|
| 235 | allocate(this%sw_dn_diffuse_surf_canopy(config%n_canopy_bands_sw,istartcol:iendcol)) |
---|
| 236 | allocate(this%sw_dn_direct_surf_canopy (config%n_canopy_bands_sw,istartcol:iendcol)) |
---|
| 237 | end if |
---|
| 238 | end if |
---|
| 239 | |
---|
| 240 | ! Allocate cloud cover arrays |
---|
| 241 | allocate(this%cloud_cover_lw(istartcol:iendcol)) |
---|
| 242 | allocate(this%cloud_cover_sw(istartcol:iendcol)) |
---|
| 243 | |
---|
| 244 | ! Some solvers may not write to cloud cover, so we initialize to |
---|
| 245 | ! an unphysical value |
---|
| 246 | this%cloud_cover_lw = -1.0_jprb |
---|
| 247 | this%cloud_cover_sw = -1.0_jprb |
---|
| 248 | |
---|
| 249 | if (lhook) call dr_hook('radiation_flux:allocate',1,hook_handle) |
---|
| 250 | |
---|
| 251 | end subroutine allocate_flux_type |
---|
| 252 | |
---|
| 253 | |
---|
| 254 | !--------------------------------------------------------------------- |
---|
| 255 | ! Deallocate flux arrays |
---|
| 256 | subroutine deallocate_flux_type(this) |
---|
| 257 | |
---|
| 258 | use yomhook, only : lhook, dr_hook |
---|
| 259 | |
---|
| 260 | class(flux_type), intent(inout) :: this |
---|
| 261 | real(jprb) :: hook_handle |
---|
| 262 | |
---|
| 263 | if (lhook) call dr_hook('radiation_flux:deallocate',0,hook_handle) |
---|
| 264 | |
---|
| 265 | if (allocated(this%lw_up)) then |
---|
| 266 | deallocate(this%lw_up) |
---|
| 267 | if (allocated(this%lw_dn)) deallocate(this%lw_dn) |
---|
| 268 | if (allocated(this%lw_up_clear)) deallocate(this%lw_up_clear) |
---|
| 269 | if (allocated(this%lw_dn_clear)) deallocate(this%lw_dn_clear) |
---|
| 270 | end if |
---|
| 271 | |
---|
| 272 | if (allocated(this%sw_up)) then |
---|
| 273 | deallocate(this%sw_up) |
---|
| 274 | if (allocated(this%sw_dn)) deallocate(this%sw_dn) |
---|
| 275 | if (allocated(this%sw_up_clear)) deallocate(this%sw_up_clear) |
---|
| 276 | if (allocated(this%sw_dn_clear)) deallocate(this%sw_dn_clear) |
---|
| 277 | if (allocated(this%sw_dn_direct)) deallocate(this%sw_dn_direct) |
---|
| 278 | if (allocated(this%sw_dn_direct_clear)) & |
---|
| 279 | & deallocate(this%sw_dn_direct_clear) |
---|
| 280 | end if |
---|
| 281 | |
---|
| 282 | if (allocated(this%lw_up_band)) then |
---|
| 283 | deallocate(this%lw_up_band) |
---|
| 284 | if (allocated(this%lw_dn_band)) deallocate(this%lw_dn_band) |
---|
| 285 | if (allocated(this%lw_up_clear_band)) deallocate(this%lw_up_clear_band) |
---|
| 286 | if (allocated(this%lw_dn_clear_band)) deallocate(this%lw_dn_clear_band) |
---|
| 287 | end if |
---|
| 288 | |
---|
| 289 | if (allocated(this%sw_up_band)) then |
---|
| 290 | deallocate(this%sw_up_band) |
---|
| 291 | if (allocated(this%sw_dn_band)) deallocate(this%sw_dn_band) |
---|
| 292 | if (allocated(this%sw_up_clear_band)) deallocate(this%sw_up_clear_band) |
---|
| 293 | if (allocated(this%sw_dn_clear_band)) deallocate(this%sw_dn_clear_band) |
---|
| 294 | if (allocated(this%sw_dn_direct_band)) deallocate(this%sw_dn_direct_band) |
---|
| 295 | if (allocated(this%sw_dn_direct_clear_band)) & |
---|
| 296 | & deallocate(this%sw_dn_direct_clear_band) |
---|
| 297 | end if |
---|
| 298 | |
---|
| 299 | if (allocated(this%sw_dn_surf_band)) then |
---|
| 300 | deallocate(this%sw_dn_surf_band) |
---|
| 301 | deallocate(this%sw_dn_direct_surf_band) |
---|
| 302 | end if |
---|
| 303 | if (allocated(this%sw_dn_surf_clear_band)) then |
---|
| 304 | deallocate(this%sw_dn_surf_clear_band) |
---|
| 305 | deallocate(this%sw_dn_direct_surf_clear_band) |
---|
| 306 | end if |
---|
| 307 | |
---|
| 308 | if (allocated(this%lw_dn_surf_canopy)) deallocate(this%lw_dn_surf_canopy) |
---|
| 309 | if (allocated(this%sw_dn_diffuse_surf_canopy)) deallocate(this%sw_dn_diffuse_surf_canopy) |
---|
| 310 | if (allocated(this%sw_dn_direct_surf_canopy)) deallocate(this%sw_dn_direct_surf_canopy) |
---|
| 311 | |
---|
| 312 | if (allocated(this%cloud_cover_sw)) then |
---|
| 313 | deallocate(this%cloud_cover_sw) |
---|
| 314 | end if |
---|
| 315 | if (allocated(this%cloud_cover_lw)) then |
---|
| 316 | deallocate(this%cloud_cover_lw) |
---|
| 317 | end if |
---|
| 318 | |
---|
| 319 | if (allocated(this%lw_derivatives)) then |
---|
| 320 | deallocate(this%lw_derivatives) |
---|
| 321 | end if |
---|
| 322 | |
---|
| 323 | if (lhook) call dr_hook('radiation_flux:deallocate',1,hook_handle) |
---|
| 324 | |
---|
| 325 | end subroutine deallocate_flux_type |
---|
| 326 | |
---|
| 327 | !--------------------------------------------------------------------- |
---|
| 328 | ! Calculate surface downwelling fluxes in each band using the |
---|
| 329 | ! downwelling surface fluxes at each g point |
---|
| 330 | subroutine calc_surface_spectral(this, config, istartcol, iendcol) |
---|
| 331 | |
---|
| 332 | use yomhook, only : lhook, dr_hook |
---|
| 333 | use radiation_config, only : config_type |
---|
| 334 | |
---|
| 335 | class(flux_type), intent(inout) :: this |
---|
| 336 | type(config_type), intent(in) :: config |
---|
| 337 | integer, intent(in) :: istartcol, iendcol |
---|
| 338 | |
---|
| 339 | integer :: jcol, jband, jalbedoband, nalbedoband |
---|
| 340 | |
---|
| 341 | ! Longwave surface downwelling in each band needed to compute |
---|
| 342 | ! canopy fluxes |
---|
| 343 | real(jprb) :: lw_dn_surf_band(config%n_bands_lw,istartcol:iendcol) |
---|
| 344 | |
---|
| 345 | real(jprb) :: hook_handle |
---|
| 346 | |
---|
| 347 | if (lhook) call dr_hook('radiation_flux:calc_surface_spectral',0,hook_handle) |
---|
| 348 | |
---|
| 349 | if (config%do_sw .and. config%do_surface_sw_spectral_flux) then |
---|
| 350 | |
---|
| 351 | do jcol = istartcol,iendcol |
---|
| 352 | call indexed_sum(this%sw_dn_direct_surf_g(:,jcol), & |
---|
| 353 | & config%i_band_from_reordered_g_sw, & |
---|
| 354 | & this%sw_dn_direct_surf_band(:,jcol)) |
---|
| 355 | call indexed_sum(this%sw_dn_diffuse_surf_g(:,jcol), & |
---|
| 356 | & config%i_band_from_reordered_g_sw, & |
---|
| 357 | & this%sw_dn_surf_band(:,jcol)) |
---|
| 358 | this%sw_dn_surf_band(:,jcol) & |
---|
| 359 | & = this%sw_dn_surf_band(:,jcol) & |
---|
| 360 | & + this%sw_dn_direct_surf_band(:,jcol) |
---|
| 361 | end do |
---|
| 362 | |
---|
| 363 | if (config%do_clear) then |
---|
| 364 | do jcol = istartcol,iendcol |
---|
| 365 | call indexed_sum(this%sw_dn_direct_surf_clear_g(:,jcol), & |
---|
| 366 | & config%i_band_from_reordered_g_sw, & |
---|
| 367 | & this%sw_dn_direct_surf_clear_band(:,jcol)) |
---|
| 368 | call indexed_sum(this%sw_dn_diffuse_surf_clear_g(:,jcol), & |
---|
| 369 | & config%i_band_from_reordered_g_sw, & |
---|
| 370 | & this%sw_dn_surf_clear_band(:,jcol)) |
---|
| 371 | this%sw_dn_surf_clear_band(:,jcol) & |
---|
| 372 | & = this%sw_dn_surf_clear_band(:,jcol) & |
---|
| 373 | & + this%sw_dn_direct_surf_clear_band(:,jcol) |
---|
| 374 | end do |
---|
| 375 | end if |
---|
| 376 | |
---|
| 377 | end if ! do_surface_sw_spectral_flux |
---|
| 378 | |
---|
| 379 | ! Fluxes in bands required for canopy radiative transfer |
---|
| 380 | if (config%do_sw .and. config%do_canopy_fluxes_sw) then |
---|
| 381 | if (config%use_canopy_full_spectrum_sw) then |
---|
| 382 | this%sw_dn_diffuse_surf_canopy(:,istartcol:iendcol) = this%sw_dn_diffuse_surf_g(:,istartcol:iendcol) |
---|
| 383 | this%sw_dn_direct_surf_canopy (:,istartcol:iendcol) = this%sw_dn_direct_surf_g (:,istartcol:iendcol) |
---|
| 384 | else if (config%do_nearest_spectral_sw_albedo) then |
---|
| 385 | do jcol = istartcol,iendcol |
---|
| 386 | call indexed_sum(this%sw_dn_direct_surf_g(:,jcol), & |
---|
| 387 | & config%i_albedo_from_band_sw(config%i_band_from_reordered_g_sw), & |
---|
| 388 | & this%sw_dn_direct_surf_canopy(:,jcol)) |
---|
| 389 | call indexed_sum(this%sw_dn_diffuse_surf_g(:,jcol), & |
---|
| 390 | & config%i_albedo_from_band_sw(config%i_band_from_reordered_g_sw), & |
---|
| 391 | & this%sw_dn_diffuse_surf_canopy(:,jcol)) |
---|
| 392 | end do |
---|
| 393 | else |
---|
| 394 | ! More accurate calculations using weights, but requires |
---|
| 395 | ! this%sw_dn_[direct_]surf_band to be defined, i.e. |
---|
| 396 | ! config%do_surface_sw_spectral_flux == .true. |
---|
| 397 | nalbedoband = size(config%sw_albedo_weights,1) |
---|
| 398 | this%sw_dn_diffuse_surf_canopy(:,istartcol:iendcol) = 0.0_jprb |
---|
| 399 | this%sw_dn_direct_surf_canopy (:,istartcol:iendcol) = 0.0_jprb |
---|
| 400 | do jband = 1,config%n_bands_sw |
---|
| 401 | do jalbedoband = 1,nalbedoband |
---|
| 402 | if (config%sw_albedo_weights(jalbedoband,jband) /= 0.0_jprb) then |
---|
| 403 | ! Initially, "diffuse" is actually "total" |
---|
| 404 | this%sw_dn_diffuse_surf_canopy(jalbedoband,istartcol:iendcol) & |
---|
| 405 | & = this%sw_dn_diffuse_surf_canopy(jalbedoband,istartcol:iendcol) & |
---|
| 406 | & + config%sw_albedo_weights(jalbedoband,jband) & |
---|
| 407 | & * this%sw_dn_surf_band(jband,istartcol:iendcol) |
---|
| 408 | this%sw_dn_direct_surf_canopy(jalbedoband,istartcol:iendcol) & |
---|
| 409 | & = this%sw_dn_direct_surf_canopy(jalbedoband,istartcol:iendcol) & |
---|
| 410 | & + config%sw_albedo_weights(jalbedoband,jband) & |
---|
| 411 | & * this%sw_dn_direct_surf_band(jband,istartcol:iendcol) |
---|
| 412 | end if |
---|
| 413 | end do |
---|
| 414 | end do |
---|
| 415 | ! Subtract the direct from total to get diffuse |
---|
| 416 | this%sw_dn_diffuse_surf_canopy(:,istartcol:iendcol) & |
---|
| 417 | & = this%sw_dn_diffuse_surf_canopy(:,istartcol:iendcol) & |
---|
| 418 | & - this%sw_dn_direct_surf_canopy(:,istartcol:iendcol) |
---|
| 419 | end if |
---|
| 420 | |
---|
| 421 | end if ! do_canopy_fluxes_sw |
---|
| 422 | |
---|
| 423 | if (config%do_lw .and. config%do_canopy_fluxes_lw) then |
---|
| 424 | if (config%use_canopy_full_spectrum_lw) then |
---|
| 425 | this%lw_dn_surf_canopy(:,istartcol:iendcol) = this%lw_dn_surf_g(:,istartcol:iendcol) |
---|
| 426 | else if (config%do_nearest_spectral_lw_emiss) then |
---|
| 427 | do jcol = istartcol,iendcol |
---|
| 428 | call indexed_sum(this%lw_dn_surf_g(:,jcol), & |
---|
| 429 | & config%i_emiss_from_band_lw(config%i_band_from_reordered_g_lw), & |
---|
| 430 | & this%lw_dn_surf_canopy(:,jcol)) |
---|
| 431 | end do |
---|
| 432 | else |
---|
| 433 | ! Compute fluxes in each longwave emissivity interval using |
---|
| 434 | ! weights; first sum over g points to get the values in bands |
---|
| 435 | do jcol = istartcol,iendcol |
---|
| 436 | call indexed_sum(this%lw_dn_surf_g(:,jcol), & |
---|
| 437 | & config%i_band_from_reordered_g_lw, & |
---|
| 438 | & lw_dn_surf_band(:,jcol)) |
---|
| 439 | end do |
---|
| 440 | nalbedoband = size(config%lw_emiss_weights,1) |
---|
| 441 | this%lw_dn_surf_canopy(:,istartcol:iendcol) = 0.0_jprb |
---|
| 442 | do jband = 1,config%n_bands_lw |
---|
| 443 | do jalbedoband = 1,nalbedoband |
---|
| 444 | if (config%lw_emiss_weights(jalbedoband,jband) /= 0.0_jprb) then |
---|
| 445 | this%lw_dn_surf_canopy(jalbedoband,istartcol:iendcol) & |
---|
| 446 | & = this%lw_dn_surf_canopy(jalbedoband,istartcol:iendcol) & |
---|
| 447 | & + config%lw_emiss_weights(jalbedoband,jband) & |
---|
| 448 | & * lw_dn_surf_band(jband,istartcol:iendcol) |
---|
| 449 | end if |
---|
| 450 | end do |
---|
| 451 | end do |
---|
| 452 | end if |
---|
| 453 | end if |
---|
| 454 | |
---|
| 455 | if (lhook) call dr_hook('radiation_flux:calc_surface_spectral',1,hook_handle) |
---|
| 456 | |
---|
| 457 | end subroutine calc_surface_spectral |
---|
| 458 | |
---|
| 459 | |
---|
| 460 | !--------------------------------------------------------------------- |
---|
| 461 | ! Return .true. if the most important flux variables are out of a |
---|
| 462 | ! physically sensible range, optionally only considering columns |
---|
| 463 | ! between istartcol and iendcol |
---|
| 464 | function out_of_physical_bounds(this, istartcol, iendcol) result(is_bad) |
---|
| 465 | |
---|
| 466 | use yomhook, only : lhook, dr_hook |
---|
| 467 | use radiation_config, only : out_of_bounds_2d |
---|
| 468 | |
---|
| 469 | class(flux_type), intent(inout) :: this |
---|
| 470 | integer, optional,intent(in) :: istartcol, iendcol |
---|
| 471 | logical :: is_bad |
---|
| 472 | |
---|
| 473 | real(jprb) :: hook_handle |
---|
| 474 | |
---|
| 475 | if (lhook) call dr_hook('radiation_flux:out_of_physical_bounds',0,hook_handle) |
---|
| 476 | |
---|
| 477 | is_bad = out_of_bounds_2d(this%lw_up, 'lw_up', 10.0_jprb, 900.0_jprb, .false., i1=istartcol, i2=iendcol) & |
---|
| 478 | & .or. out_of_bounds_2d(this%lw_dn, 'lw_dn', 0.0_jprb, 800.0_jprb, .false., i1=istartcol, i2=iendcol) & |
---|
| 479 | & .or. out_of_bounds_2d(this%sw_up, 'sw_up', 0.0_jprb, 1500.0_jprb, .false., i1=istartcol, i2=iendcol) & |
---|
| 480 | & .or. out_of_bounds_2d(this%sw_dn, 'sw_dn', 0.0_jprb, 1500.0_jprb, .false., i1=istartcol, i2=iendcol) & |
---|
| 481 | & .or. out_of_bounds_2d(this%sw_dn_direct, 'sw_dn_direct', 0.0_jprb, 1500.0_jprb, .false., i1=istartcol, i2=iendcol) & |
---|
| 482 | & .or. out_of_bounds_2d(this%lw_derivatives, 'lw_derivatives', 0.0_jprb, 1.0_jprb, .false., i1=istartcol, i2=iendcol) & |
---|
| 483 | & .or. out_of_bounds_2d(this%sw_dn_surf_band, 'sw_dn_surf_band', 0.0_jprb, 1500.0_jprb, & |
---|
| 484 | & .false., j1=istartcol, j2=iendcol) & |
---|
| 485 | & .or. out_of_bounds_2d(this%sw_dn_surf_clear_band, 'sw_dn_surf_clear_band', 0.0_jprb, 1500.0_jprb, & |
---|
| 486 | & .false., j1=istartcol, j2=iendcol) |
---|
| 487 | |
---|
| 488 | if (lhook) call dr_hook('radiation_flux:out_of_physical_bounds',1,hook_handle) |
---|
| 489 | |
---|
| 490 | end function out_of_physical_bounds |
---|
| 491 | |
---|
| 492 | !--------------------------------------------------------------------- |
---|
| 493 | ! Return .true. if the heating rates are out of a physically |
---|
| 494 | ! sensible range, optionally only considering columns between |
---|
| 495 | ! istartcol and iendcol. This function allocates and deallocates |
---|
| 496 | ! memory due to the requirements for inputs of out_of_bounds_2d. |
---|
| 497 | function heating_rate_out_of_physical_bounds(this, nlev, istartcol, iendcol, pressure_hl) result(is_bad) |
---|
| 498 | |
---|
| 499 | use radiation_config, only : out_of_bounds_2d |
---|
| 500 | use radiation_constants, only : AccelDueToGravity |
---|
| 501 | |
---|
| 502 | ! "Cp" (J kg-1 K-1) |
---|
| 503 | real(jprb), parameter :: SpecificHeatDryAir = 1004.0 |
---|
| 504 | |
---|
| 505 | class(flux_type), intent(inout) :: this |
---|
| 506 | integer, intent(in) :: istartcol, iendcol, nlev |
---|
| 507 | logical :: is_bad |
---|
| 508 | |
---|
| 509 | real(jprb), intent(in) :: pressure_hl(:,:) |
---|
| 510 | |
---|
| 511 | real(jprb), allocatable :: hr_K_day(:,:) |
---|
| 512 | |
---|
| 513 | real(jprb) :: scaling(istartcol:iendcol,nlev) |
---|
| 514 | |
---|
| 515 | allocate(hr_K_day(istartcol:iendcol,nlev)) |
---|
| 516 | |
---|
| 517 | scaling = -(24.0_jprb * 3600.0_jprb * AccelDueToGravity / SpecificHeatDryAir) & |
---|
| 518 | & / (pressure_hl(istartcol:iendcol,2:nlev+1) - pressure_hl(istartcol:iendcol,1:nlev)) |
---|
| 519 | ! Shortwave |
---|
| 520 | hr_K_day = scaling * (this%sw_dn(istartcol:iendcol,2:nlev+1) - this%sw_up(istartcol:iendcol,2:nlev+1) & |
---|
| 521 | & -this%sw_dn(istartcol:iendcol,1:nlev) + this%sw_up(istartcol:iendcol,1:nlev)) |
---|
| 522 | is_bad = out_of_bounds_2d(hr_K_day, 'sw_heating_rate_K_day', 0.0_jprb, 200.0_jprb, & |
---|
| 523 | & .false., i1=istartcol, i2=iendcol) |
---|
| 524 | |
---|
| 525 | ! Longwave |
---|
| 526 | hr_K_day = scaling * (this%lw_dn(istartcol:iendcol,2:nlev+1) - this%lw_up(istartcol:iendcol,2:nlev+1) & |
---|
| 527 | & -this%lw_dn(istartcol:iendcol,1:nlev) + this%lw_up(istartcol:iendcol,1:nlev)) |
---|
| 528 | is_bad = is_bad .or. out_of_bounds_2d(hr_K_day, 'lw_heating_rate_K_day', -250.0_jprb, 150.0_jprb, & |
---|
| 529 | & .false., i1=istartcol, i2=iendcol) |
---|
| 530 | |
---|
| 531 | deallocate(hr_K_day) |
---|
| 532 | |
---|
| 533 | end function heating_rate_out_of_physical_bounds |
---|
| 534 | |
---|
| 535 | |
---|
| 536 | !--------------------------------------------------------------------- |
---|
| 537 | ! Sum elements of "source" into "dest" according to index "ind". |
---|
| 538 | ! "source" and "ind" should have the same size and bounds, and no |
---|
| 539 | ! element of "ind" should refer outside the bounds of "dest". This |
---|
| 540 | ! version increments existing contents of "dest". |
---|
| 541 | pure subroutine add_indexed_sum(source, ind, dest) |
---|
| 542 | |
---|
| 543 | real(jprb), intent(in) :: source(:) |
---|
| 544 | integer, intent(in) :: ind(:) |
---|
| 545 | real(jprb), intent(inout) :: dest(:) |
---|
| 546 | |
---|
| 547 | integer :: ig, jg, istart, iend |
---|
| 548 | |
---|
| 549 | istart = lbound(source,1) |
---|
| 550 | iend = ubound(source,1) |
---|
| 551 | |
---|
| 552 | do jg = istart, iend |
---|
| 553 | ig = ind(jg) |
---|
| 554 | dest(ig) = dest(ig) + source(jg) |
---|
| 555 | end do |
---|
| 556 | |
---|
| 557 | end subroutine add_indexed_sum |
---|
| 558 | |
---|
| 559 | |
---|
| 560 | !--------------------------------------------------------------------- |
---|
| 561 | ! As "add_indexed_sum" but this version overwrites existing contents |
---|
| 562 | ! of "dest" |
---|
| 563 | pure subroutine indexed_sum(source, ind, dest) |
---|
| 564 | |
---|
| 565 | real(jprb), intent(in) :: source(:) |
---|
| 566 | integer, intent(in) :: ind(:) |
---|
| 567 | real(jprb), intent(out) :: dest(:) |
---|
| 568 | |
---|
| 569 | integer :: ig, jg, istart, iend |
---|
| 570 | |
---|
| 571 | dest = 0.0 |
---|
| 572 | |
---|
| 573 | istart = lbound(source,1) |
---|
| 574 | iend = ubound(source,1) |
---|
| 575 | |
---|
| 576 | do jg = istart, iend |
---|
| 577 | ig = ind(jg) |
---|
| 578 | dest(ig) = dest(ig) + source(jg) |
---|
| 579 | end do |
---|
| 580 | |
---|
| 581 | end subroutine indexed_sum |
---|
| 582 | |
---|
| 583 | |
---|
| 584 | !--------------------------------------------------------------------- |
---|
| 585 | ! As "add_indexed_sum" but a whole vertical profiles |
---|
| 586 | pure subroutine add_indexed_sum_profile(source, ind, dest) |
---|
| 587 | |
---|
| 588 | real(jprb), intent(in) :: source(:,:) |
---|
| 589 | integer, intent(in) :: ind(:) |
---|
| 590 | real(jprb), intent(out) :: dest(:,:) |
---|
| 591 | |
---|
| 592 | integer :: ig, jg, istart, iend, jlev, nlev |
---|
| 593 | |
---|
| 594 | istart = lbound(source,1) |
---|
| 595 | iend = ubound(source,1) |
---|
| 596 | nlev = size(source,2) |
---|
| 597 | |
---|
| 598 | do jlev = 1,nlev |
---|
| 599 | do jg = istart, iend |
---|
| 600 | ig = ind(jg) |
---|
| 601 | dest(ig,jlev) = dest(ig,jlev) + source(jg,jlev) |
---|
| 602 | end do |
---|
| 603 | end do |
---|
| 604 | |
---|
| 605 | end subroutine add_indexed_sum_profile |
---|
| 606 | |
---|
| 607 | |
---|
| 608 | !--------------------------------------------------------------------- |
---|
| 609 | ! As "indexed_sum" but a whole vertical profiles |
---|
| 610 | pure subroutine indexed_sum_profile(source, ind, dest) |
---|
| 611 | |
---|
| 612 | real(jprb), intent(in) :: source(:,:) |
---|
| 613 | integer, intent(in) :: ind(:) |
---|
| 614 | real(jprb), intent(out) :: dest(:,:) |
---|
| 615 | |
---|
| 616 | integer :: ig, jg, istart, iend, jlev, nlev |
---|
| 617 | |
---|
| 618 | dest = 0.0 |
---|
| 619 | |
---|
| 620 | istart = lbound(source,1) |
---|
| 621 | iend = ubound(source,1) |
---|
| 622 | nlev = size(source,2) |
---|
| 623 | |
---|
| 624 | do jlev = 1,nlev |
---|
| 625 | do jg = istart, iend |
---|
| 626 | ig = ind(jg) |
---|
| 627 | dest(ig,jlev) = dest(ig,jlev) + source(jg,jlev) |
---|
| 628 | end do |
---|
| 629 | end do |
---|
| 630 | |
---|
| 631 | end subroutine indexed_sum_profile |
---|
| 632 | |
---|
| 633 | end module radiation_flux |
---|