[4773] | 1 | ! radiation_lw_derivatives.F90 - Compute longwave derivatives for Hogan and Bozzo (2015) method |
---|
| 2 | ! |
---|
| 3 | ! (C) Copyright 2016- 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 | ! This module provides routines to compute the rate of change of |
---|
| 16 | ! broadband upwelling longwave flux at each half level with respect to |
---|
| 17 | ! the surface broadband upwelling flux. This is done from the surface |
---|
| 18 | ! spectral fluxes and the spectral transmittance of each atmospheric |
---|
| 19 | ! layer, assuming no longwave scattering. The result may be used to |
---|
| 20 | ! perform approximate updates to the longwave flux profile in between |
---|
| 21 | ! calls to the full radiation scheme, accounting for the change in |
---|
| 22 | ! skin temperature, following the method of Hogan and Bozzo (JAMES |
---|
| 23 | ! 2015). Separate routines are provided for each solver. |
---|
| 24 | ! |
---|
| 25 | ! Note that currently a more approximate calculation is performed from |
---|
| 26 | ! the exact one in Hogan and Bozzo (2015); here we assume that a |
---|
| 27 | ! change in temperature increases the spectral fluxes in proportion, |
---|
| 28 | ! when in reality there is a change in shape of the Planck function in |
---|
| 29 | ! addition to an overall increase in the total emission. |
---|
| 30 | ! |
---|
| 31 | ! Modifications |
---|
| 32 | ! 2017-10-23 R. Hogan Renamed single-character variables |
---|
| 33 | ! 2022-11-22 P. Ukkonen / R. Hogan Optimized calc_lw_derivatives_region |
---|
| 34 | |
---|
| 35 | module radiation_lw_derivatives |
---|
| 36 | |
---|
| 37 | public |
---|
| 38 | |
---|
| 39 | contains |
---|
| 40 | |
---|
| 41 | !--------------------------------------------------------------------- |
---|
| 42 | ! Calculation for the Independent Column Approximation |
---|
| 43 | subroutine calc_lw_derivatives_ica(ng, nlev, icol, transmittance, flux_up_surf, lw_derivatives) |
---|
| 44 | |
---|
| 45 | use parkind1, only : jprb |
---|
| 46 | use yomhook, only : lhook, dr_hook, jphook |
---|
| 47 | |
---|
| 48 | implicit none |
---|
| 49 | |
---|
| 50 | ! Inputs |
---|
| 51 | integer, intent(in) :: ng ! number of spectral intervals |
---|
| 52 | integer, intent(in) :: nlev ! number of levels |
---|
| 53 | integer, intent(in) :: icol ! Index of column for output |
---|
| 54 | real(jprb), intent(in) :: transmittance(ng,nlev) |
---|
| 55 | real(jprb), intent(in) :: flux_up_surf(ng) ! Upwelling surface spectral flux (W m-2) |
---|
| 56 | |
---|
| 57 | ! Output |
---|
| 58 | real(jprb), intent(out) :: lw_derivatives(:,:) ! dimensioned (ncol,nlev+1) |
---|
| 59 | |
---|
| 60 | ! Rate of change of spectral flux at a given height with respect |
---|
| 61 | ! to the surface value |
---|
| 62 | real(jprb) :: lw_derivatives_g(ng) |
---|
| 63 | |
---|
| 64 | integer :: jlev |
---|
| 65 | |
---|
| 66 | real(jphook) :: hook_handle |
---|
| 67 | |
---|
| 68 | if (lhook) call dr_hook('radiation_lw_derivatives:calc_lw_derivatives_ica',0,hook_handle) |
---|
| 69 | |
---|
| 70 | ! Initialize the derivatives at the surface |
---|
| 71 | lw_derivatives_g = flux_up_surf / sum(flux_up_surf) |
---|
| 72 | lw_derivatives(icol, nlev+1) = 1.0_jprb |
---|
| 73 | |
---|
| 74 | ! Move up through the atmosphere computing the derivatives at each |
---|
| 75 | ! half-level |
---|
| 76 | do jlev = nlev,1,-1 |
---|
| 77 | lw_derivatives_g = lw_derivatives_g * transmittance(:,jlev) |
---|
| 78 | lw_derivatives(icol,jlev) = sum(lw_derivatives_g) |
---|
| 79 | end do |
---|
| 80 | |
---|
| 81 | if (lhook) call dr_hook('radiation_lw_derivatives:calc_lw_derivatives_ica',1,hook_handle) |
---|
| 82 | |
---|
| 83 | end subroutine calc_lw_derivatives_ica |
---|
| 84 | |
---|
| 85 | |
---|
| 86 | !--------------------------------------------------------------------- |
---|
| 87 | ! Calculation for the Independent Column Approximation |
---|
| 88 | subroutine modify_lw_derivatives_ica(ng, nlev, icol, transmittance, & |
---|
| 89 | & flux_up_surf, weight, lw_derivatives) |
---|
| 90 | |
---|
| 91 | use parkind1, only : jprb |
---|
| 92 | use yomhook, only : lhook, dr_hook, jphook |
---|
| 93 | |
---|
| 94 | implicit none |
---|
| 95 | |
---|
| 96 | ! Inputs |
---|
| 97 | integer, intent(in) :: ng ! number of spectral intervals |
---|
| 98 | integer, intent(in) :: nlev ! number of levels |
---|
| 99 | integer, intent(in) :: icol ! Index of column for output |
---|
| 100 | real(jprb), intent(in) :: transmittance(ng,nlev) |
---|
| 101 | real(jprb), intent(in) :: flux_up_surf(ng) ! Upwelling surface spectral flux (W m-2) |
---|
| 102 | real(jprb), intent(in) :: weight ! Weight new values against existing |
---|
| 103 | |
---|
| 104 | ! Output |
---|
| 105 | real(jprb), intent(inout) :: lw_derivatives(:,:) ! dimensioned (ncol,nlev+1) |
---|
| 106 | |
---|
| 107 | ! Rate of change of spectral flux at a given height with respect |
---|
| 108 | ! to the surface value |
---|
| 109 | real(jprb) :: lw_derivatives_g(ng) |
---|
| 110 | |
---|
| 111 | integer :: jlev |
---|
| 112 | |
---|
| 113 | real(jphook) :: hook_handle |
---|
| 114 | |
---|
| 115 | if (lhook) call dr_hook('radiation_lw_derivatives:modify_lw_derivatives_ica',0,hook_handle) |
---|
| 116 | |
---|
| 117 | ! Initialize the derivatives at the surface |
---|
| 118 | lw_derivatives_g = flux_up_surf / sum(flux_up_surf) |
---|
| 119 | ! This value must be 1 so no weighting applied |
---|
| 120 | lw_derivatives(icol, nlev+1) = 1.0_jprb |
---|
| 121 | |
---|
| 122 | ! Move up through the atmosphere computing the derivatives at each |
---|
| 123 | ! half-level |
---|
| 124 | do jlev = nlev,1,-1 |
---|
| 125 | lw_derivatives_g = lw_derivatives_g * transmittance(:,jlev) |
---|
| 126 | lw_derivatives(icol,jlev) = (1.0_jprb - weight) * lw_derivatives(icol,jlev) & |
---|
| 127 | & + weight * sum(lw_derivatives_g) |
---|
| 128 | end do |
---|
| 129 | |
---|
| 130 | if (lhook) call dr_hook('radiation_lw_derivatives:modify_lw_derivatives_ica',1,hook_handle) |
---|
| 131 | |
---|
| 132 | end subroutine modify_lw_derivatives_ica |
---|
| 133 | |
---|
| 134 | |
---|
| 135 | |
---|
| 136 | !--------------------------------------------------------------------- |
---|
| 137 | ! Calculation for solvers involving multiple regions and matrices |
---|
| 138 | subroutine calc_lw_derivatives_matrix(ng, nlev, nreg, icol, transmittance, & |
---|
| 139 | & u_matrix, flux_up_surf, lw_derivatives) |
---|
| 140 | |
---|
| 141 | use parkind1, only : jprb |
---|
| 142 | use yomhook, only : lhook, dr_hook, jphook |
---|
| 143 | |
---|
| 144 | use radiation_matrix |
---|
| 145 | |
---|
| 146 | implicit none |
---|
| 147 | |
---|
| 148 | ! Inputs |
---|
| 149 | integer, intent(in) :: ng ! number of spectral intervals |
---|
| 150 | integer, intent(in) :: nlev ! number of levels |
---|
| 151 | integer, intent(in) :: nreg ! number of regions |
---|
| 152 | integer, intent(in) :: icol ! Index of column for output |
---|
| 153 | real(jprb), intent(in) :: transmittance(ng,nreg,nreg,nlev) |
---|
| 154 | real(jprb), intent(in) :: u_matrix(nreg,nreg,nlev+1) ! Upward overlap matrix |
---|
| 155 | real(jprb), intent(in) :: flux_up_surf(ng) ! Upwelling surface spectral flux (W m-2) |
---|
| 156 | |
---|
| 157 | ! Output |
---|
| 158 | real(jprb), intent(out) :: lw_derivatives(:,:) ! dimensioned (ncol,nlev+1) |
---|
| 159 | |
---|
| 160 | ! Rate of change of spectral flux at a given height with respect |
---|
| 161 | ! to the surface value |
---|
| 162 | real(jprb) :: lw_derivatives_g_reg(ng,nreg) |
---|
| 163 | |
---|
| 164 | integer :: jlev |
---|
| 165 | |
---|
| 166 | real(jphook) :: hook_handle |
---|
| 167 | |
---|
| 168 | if (lhook) call dr_hook('radiation_lw_derivatives:calc_lw_derivatives_matrix',0,hook_handle) |
---|
| 169 | |
---|
| 170 | ! Initialize the derivatives at the surface; the surface is |
---|
| 171 | ! treated as a single clear-sky layer so we only need to put |
---|
| 172 | ! values in region 1. |
---|
| 173 | lw_derivatives_g_reg = 0.0_jprb |
---|
| 174 | lw_derivatives_g_reg(:,1) = flux_up_surf / sum(flux_up_surf) |
---|
| 175 | lw_derivatives(icol, nlev+1) = 1.0_jprb |
---|
| 176 | |
---|
| 177 | ! Move up through the atmosphere computing the derivatives at each |
---|
| 178 | ! half-level |
---|
| 179 | do jlev = nlev,1,-1 |
---|
| 180 | ! Compute effect of overlap at half-level jlev+1, yielding |
---|
| 181 | ! derivatives just above that half-level |
---|
| 182 | lw_derivatives_g_reg = singlemat_x_vec(ng,ng,nreg,u_matrix(:,:,jlev+1),lw_derivatives_g_reg) |
---|
| 183 | |
---|
| 184 | ! Compute effect of transmittance of layer jlev, yielding |
---|
| 185 | ! derivatives just below the half-level above (jlev) |
---|
| 186 | lw_derivatives_g_reg = mat_x_vec(ng,ng,nreg,transmittance(:,:,:,jlev),lw_derivatives_g_reg) |
---|
| 187 | |
---|
| 188 | lw_derivatives(icol, jlev) = sum(lw_derivatives_g_reg) |
---|
| 189 | end do |
---|
| 190 | |
---|
| 191 | if (lhook) call dr_hook('radiation_lw_derivatives:calc_lw_derivatives_matrix',1,hook_handle) |
---|
| 192 | |
---|
| 193 | end subroutine calc_lw_derivatives_matrix |
---|
| 194 | |
---|
| 195 | |
---|
| 196 | !--------------------------------------------------------------------- |
---|
| 197 | ! Calculation for solvers involving multiple regions but no 3D |
---|
| 198 | ! effects: the difference from calc_lw_derivatives_matrix is that transmittance |
---|
| 199 | ! has one fewer dimensions |
---|
| 200 | subroutine calc_lw_derivatives_region(ng, nlev, nreg, icol, transmittance, & |
---|
| 201 | & u_matrix, flux_up_surf, lw_derivatives) |
---|
| 202 | |
---|
| 203 | use parkind1, only : jprb |
---|
| 204 | use yomhook, only : lhook, dr_hook, jphook |
---|
| 205 | |
---|
| 206 | use radiation_matrix |
---|
| 207 | |
---|
| 208 | implicit none |
---|
| 209 | |
---|
| 210 | ! Inputs |
---|
| 211 | integer, intent(in) :: ng ! number of spectral intervals |
---|
| 212 | integer, intent(in) :: nlev ! number of levels |
---|
| 213 | integer, intent(in) :: nreg ! number of regions |
---|
| 214 | integer, intent(in) :: icol ! Index of column for output |
---|
| 215 | real(jprb), intent(in) :: transmittance(ng,nreg,nlev) |
---|
| 216 | real(jprb), intent(in) :: u_matrix(nreg,nreg,nlev+1) ! Upward overlap matrix |
---|
| 217 | real(jprb), intent(in) :: flux_up_surf(ng) ! Upwelling surface spectral flux (W m-2) |
---|
| 218 | |
---|
| 219 | ! Output |
---|
| 220 | real(jprb), intent(out) :: lw_derivatives(:,:) ! dimensioned (ncol,nlev+1) |
---|
| 221 | |
---|
| 222 | ! Rate of change of spectral flux at a given height with respect |
---|
| 223 | ! to the surface value |
---|
| 224 | real(jprb) :: lw_deriv(ng,nreg), lw_deriv_below(ng,nreg) |
---|
| 225 | real(jprb) :: partial_sum(ng) |
---|
| 226 | |
---|
| 227 | integer :: jlev, jg |
---|
| 228 | |
---|
| 229 | real(jphook) :: hook_handle |
---|
| 230 | |
---|
| 231 | if (lhook) call dr_hook('radiation_lw_derivatives:calc_lw_derivatives_region',0,hook_handle) |
---|
| 232 | |
---|
| 233 | ! Initialize the derivatives at the surface; the surface is |
---|
| 234 | ! treated as a single clear-sky layer so we only need to put |
---|
| 235 | ! values in region 1. |
---|
| 236 | lw_deriv = 0.0_jprb |
---|
| 237 | lw_deriv(:,1) = flux_up_surf / sum(flux_up_surf) |
---|
| 238 | lw_derivatives(icol, nlev+1) = 1.0_jprb |
---|
| 239 | |
---|
| 240 | if (nreg == 3) then |
---|
| 241 | ! Optimize the most common case of 3 regions by removing the |
---|
| 242 | ! nested call to singlemat_x_vec and unrolling the matrix |
---|
| 243 | ! multiplication inline |
---|
| 244 | |
---|
| 245 | do jlev = nlev,1,-1 |
---|
| 246 | ! Compute effect of overlap at half-level jlev+1, yielding |
---|
| 247 | ! derivatives just above that half-level |
---|
| 248 | lw_deriv_below = lw_deriv |
---|
| 249 | |
---|
| 250 | associate(A=>u_matrix(:,:,jlev+1), b=>lw_deriv_below) |
---|
| 251 | do jg = 1,ng |
---|
| 252 | ! Both inner and outer loop of the matrix loops j1 and j2 unrolled |
---|
| 253 | ! inner loop: j2=1 j2=2 j2=3 |
---|
| 254 | lw_deriv(jg,1) = A(1,1)*b(jg,1) + A(1,2)*b(jg,2) + A(1,3)*b(jg,3) |
---|
| 255 | lw_deriv(jg,2) = A(2,1)*b(jg,1) + A(2,2)*b(jg,2) + A(2,3)*b(jg,3) |
---|
| 256 | lw_deriv(jg,3) = A(3,1)*b(jg,1) + A(3,2)*b(jg,2) + A(3,3)*b(jg,3) |
---|
| 257 | |
---|
| 258 | ! Compute effect of transmittance of layer jlev, yielding |
---|
| 259 | ! derivatives just below the half-level above (jlev) |
---|
| 260 | lw_deriv(jg,1) = lw_deriv(jg,1) * transmittance(jg,1,jlev) |
---|
| 261 | lw_deriv(jg,2) = lw_deriv(jg,2) * transmittance(jg,2,jlev) |
---|
| 262 | lw_deriv(jg,3) = lw_deriv(jg,3) * transmittance(jg,3,jlev) |
---|
| 263 | |
---|
| 264 | partial_sum(jg) = lw_deriv(jg,1) + lw_deriv(jg,2) + lw_deriv(jg,3) |
---|
| 265 | end do |
---|
| 266 | end associate |
---|
| 267 | |
---|
| 268 | lw_derivatives(icol, jlev) = sum(partial_sum) |
---|
| 269 | end do |
---|
| 270 | else |
---|
| 271 | ! General case when number of regions is not 3 |
---|
| 272 | |
---|
| 273 | ! Move up through the atmosphere computing the derivatives at each |
---|
| 274 | ! half-level |
---|
| 275 | do jlev = nlev,1,-1 |
---|
| 276 | ! Compute effect of overlap at half-level jlev+1, yielding |
---|
| 277 | ! derivatives just above that half-level |
---|
| 278 | lw_deriv = singlemat_x_vec(ng,ng,nreg,u_matrix(:,:,jlev+1),lw_deriv) |
---|
| 279 | |
---|
| 280 | ! Compute effect of transmittance of layer jlev, yielding |
---|
| 281 | ! derivatives just below the half-level above (jlev) |
---|
| 282 | lw_deriv = transmittance(:,:,jlev) * lw_deriv |
---|
| 283 | |
---|
| 284 | lw_derivatives(icol, jlev) = sum(lw_deriv) |
---|
| 285 | end do |
---|
| 286 | end if |
---|
| 287 | |
---|
| 288 | if (lhook) call dr_hook('radiation_lw_derivatives:calc_lw_derivatives_region',1,hook_handle) |
---|
| 289 | |
---|
| 290 | end subroutine calc_lw_derivatives_region |
---|
| 291 | |
---|
| 292 | |
---|
| 293 | end module radiation_lw_derivatives |
---|