| 1 | ! radiation_adding_ica_lw.F90 - Longwave adding method in independent column approximation |
|---|
| 2 | ! |
|---|
| 3 | ! (C) Copyright 2015- 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-04-11 R. Hogan Receive emission/albedo rather than planck/emissivity |
|---|
| 17 | ! 2017-07-12 R. Hogan Fast adding method for if only clouds scatter |
|---|
| 18 | ! 2017-10-23 R. Hogan Renamed single-character variables |
|---|
| 19 | |
|---|
| 20 | module radiation_adding_ica_lw |
|---|
| 21 | |
|---|
| 22 | public |
|---|
| 23 | |
|---|
| 24 | contains |
|---|
| 25 | |
|---|
| 26 | !--------------------------------------------------------------------- |
|---|
| 27 | ! Use the scalar "adding" method to compute longwave flux profiles, |
|---|
| 28 | ! including scattering, by successively adding the contribution of |
|---|
| 29 | ! layers starting from the surface to compute the total albedo and |
|---|
| 30 | ! total upward emission of the increasingly larger block of |
|---|
| 31 | ! atmospheric layers. |
|---|
| 32 | subroutine adding_ica_lw(ncol, nlev, & |
|---|
| 33 | & reflectance, transmittance, source_up, source_dn, emission_surf, albedo_surf, & |
|---|
| 34 | & flux_up, flux_dn) |
|---|
| 35 | |
|---|
| 36 | use parkind1, only : jprb |
|---|
| 37 | use yomhook, only : lhook, dr_hook |
|---|
| 38 | |
|---|
| 39 | implicit none |
|---|
| 40 | |
|---|
| 41 | ! Inputs |
|---|
| 42 | integer, intent(in) :: ncol ! number of columns (may be spectral intervals) |
|---|
| 43 | integer, intent(in) :: nlev ! number of levels |
|---|
| 44 | |
|---|
| 45 | ! Surface emission (W m-2) and albedo |
|---|
| 46 | real(jprb), intent(in), dimension(ncol) :: emission_surf, albedo_surf |
|---|
| 47 | |
|---|
| 48 | ! Diffuse reflectance and transmittance of each layer |
|---|
| 49 | real(jprb), intent(in), dimension(ncol, nlev) :: reflectance, transmittance |
|---|
| 50 | |
|---|
| 51 | ! Emission from each layer in an upward and downward direction |
|---|
| 52 | real(jprb), intent(in), dimension(ncol, nlev) :: source_up, source_dn |
|---|
| 53 | |
|---|
| 54 | ! Resulting fluxes (W m-2) at half-levels: diffuse upwelling and |
|---|
| 55 | ! downwelling |
|---|
| 56 | real(jprb), intent(out), dimension(ncol, nlev+1) :: flux_up, flux_dn |
|---|
| 57 | |
|---|
| 58 | ! Albedo of the entire earth/atmosphere system below each half |
|---|
| 59 | ! level |
|---|
| 60 | real(jprb), dimension(ncol, nlev+1) :: albedo |
|---|
| 61 | |
|---|
| 62 | ! Upwelling radiation at each half-level due to emission below |
|---|
| 63 | ! that half-level (W m-2) |
|---|
| 64 | real(jprb), dimension(ncol, nlev+1) :: source |
|---|
| 65 | |
|---|
| 66 | ! Equal to 1/(1-albedo*reflectance) |
|---|
| 67 | real(jprb), dimension(ncol, nlev) :: inv_denominator |
|---|
| 68 | |
|---|
| 69 | ! Loop index for model level and column |
|---|
| 70 | integer :: jlev, jcol |
|---|
| 71 | |
|---|
| 72 | real(jprb) :: hook_handle |
|---|
| 73 | |
|---|
| 74 | if (lhook) call dr_hook('radiation_adding_ica_lw:adding_ica_lw',0,hook_handle) |
|---|
| 75 | |
|---|
| 76 | albedo(:,nlev+1) = albedo_surf |
|---|
| 77 | |
|---|
| 78 | ! At the surface, the source is thermal emission |
|---|
| 79 | source(:,nlev+1) = emission_surf |
|---|
| 80 | |
|---|
| 81 | ! Work back up through the atmosphere and compute the albedo of |
|---|
| 82 | ! the entire earth/atmosphere system below that half-level, and |
|---|
| 83 | ! also the "source", which is the upwelling flux due to emission |
|---|
| 84 | ! below that level |
|---|
| 85 | do jlev = nlev,1,-1 |
|---|
| 86 | ! Next loop over columns. We could do this by indexing the |
|---|
| 87 | ! entire inner dimension as follows, e.g. for the first line: |
|---|
| 88 | ! inv_denominator(:,jlev) = 1.0_jprb / (1.0_jprb-albedo(:,jlev+1)*reflectance(:,jlev)) |
|---|
| 89 | ! and similarly for subsequent lines, but this slows down the |
|---|
| 90 | ! routine by a factor of 2! Rather, we do it with an explicit |
|---|
| 91 | ! loop. |
|---|
| 92 | do jcol = 1,ncol |
|---|
| 93 | ! Lacis and Hansen (1974) Eq 33, Shonk & Hogan (2008) Eq 10: |
|---|
| 94 | inv_denominator(jcol,jlev) = 1.0_jprb & |
|---|
| 95 | & / (1.0_jprb-albedo(jcol,jlev+1)*reflectance(jcol,jlev)) |
|---|
| 96 | ! Shonk & Hogan (2008) Eq 9, Petty (2006) Eq 13.81: |
|---|
| 97 | albedo(jcol,jlev) = reflectance(jcol,jlev) + transmittance(jcol,jlev)*transmittance(jcol,jlev) & |
|---|
| 98 | & * albedo(jcol,jlev+1) * inv_denominator(jcol,jlev) |
|---|
| 99 | ! Shonk & Hogan (2008) Eq 11: |
|---|
| 100 | source(jcol,jlev) = source_up(jcol,jlev) & |
|---|
| 101 | & + transmittance(jcol,jlev) * (source(jcol,jlev+1) & |
|---|
| 102 | & + albedo(jcol,jlev+1)*source_dn(jcol,jlev)) & |
|---|
| 103 | & * inv_denominator(jcol,jlev) |
|---|
| 104 | end do |
|---|
| 105 | end do |
|---|
| 106 | |
|---|
| 107 | ! At top-of-atmosphere there is no diffuse downwelling radiation |
|---|
| 108 | flux_dn(:,1) = 0.0_jprb |
|---|
| 109 | |
|---|
| 110 | ! At top-of-atmosphere, all upwelling radiation is due to emission |
|---|
| 111 | ! below that level |
|---|
| 112 | flux_up(:,1) = source(:,1) |
|---|
| 113 | |
|---|
| 114 | ! Work back down through the atmosphere computing the fluxes at |
|---|
| 115 | ! each half-level |
|---|
| 116 | do jlev = 1,nlev |
|---|
| 117 | do jcol = 1,ncol |
|---|
| 118 | ! Shonk & Hogan (2008) Eq 14 (after simplification): |
|---|
| 119 | flux_dn(jcol,jlev+1) & |
|---|
| 120 | & = (transmittance(jcol,jlev)*flux_dn(jcol,jlev) & |
|---|
| 121 | & + reflectance(jcol,jlev)*source(jcol,jlev+1) & |
|---|
| 122 | & + source_dn(jcol,jlev)) * inv_denominator(jcol,jlev) |
|---|
| 123 | ! Shonk & Hogan (2008) Eq 12: |
|---|
| 124 | flux_up(jcol,jlev+1) = albedo(jcol,jlev+1)*flux_dn(jcol,jlev+1) & |
|---|
| 125 | & + source(jcol,jlev+1) |
|---|
| 126 | end do |
|---|
| 127 | end do |
|---|
| 128 | |
|---|
| 129 | if (lhook) call dr_hook('radiation_adding_ica_lw:adding_ica_lw',1,hook_handle) |
|---|
| 130 | |
|---|
| 131 | end subroutine adding_ica_lw |
|---|
| 132 | |
|---|
| 133 | |
|---|
| 134 | !--------------------------------------------------------------------- |
|---|
| 135 | ! Use the scalar "adding" method to compute longwave flux profiles, |
|---|
| 136 | ! including scattering in cloudy layers only. |
|---|
| 137 | subroutine fast_adding_ica_lw(ncol, nlev, & |
|---|
| 138 | & reflectance, transmittance, source_up, source_dn, emission_surf, albedo_surf, & |
|---|
| 139 | & is_clear_sky_layer, i_cloud_top, flux_dn_clear, & |
|---|
| 140 | & flux_up, flux_dn) |
|---|
| 141 | |
|---|
| 142 | use parkind1, only : jprb |
|---|
| 143 | use yomhook, only : lhook, dr_hook |
|---|
| 144 | |
|---|
| 145 | implicit none |
|---|
| 146 | |
|---|
| 147 | ! Inputs |
|---|
| 148 | integer, intent(in) :: ncol ! number of columns (may be spectral intervals) |
|---|
| 149 | integer, intent(in) :: nlev ! number of levels |
|---|
| 150 | |
|---|
| 151 | ! Surface emission (W m-2) and albedo |
|---|
| 152 | real(jprb), intent(in), dimension(ncol) :: emission_surf, albedo_surf |
|---|
| 153 | |
|---|
| 154 | ! Diffuse reflectance and transmittance of each layer |
|---|
| 155 | real(jprb), intent(in), dimension(ncol, nlev) :: reflectance, transmittance |
|---|
| 156 | |
|---|
| 157 | ! Emission from each layer in an upward and downward direction |
|---|
| 158 | real(jprb), intent(in), dimension(ncol, nlev) :: source_up, source_dn |
|---|
| 159 | |
|---|
| 160 | ! Determine which layers are cloud-free |
|---|
| 161 | logical, intent(in) :: is_clear_sky_layer(nlev) |
|---|
| 162 | |
|---|
| 163 | ! Index to highest cloudy layer |
|---|
| 164 | integer, intent(in) :: i_cloud_top |
|---|
| 165 | |
|---|
| 166 | ! Pre-computed clear-sky downwelling fluxes (W m-2) at half-levels |
|---|
| 167 | real(jprb), intent(in), dimension(ncol, nlev+1) :: flux_dn_clear |
|---|
| 168 | |
|---|
| 169 | ! Resulting fluxes (W m-2) at half-levels: diffuse upwelling and |
|---|
| 170 | ! downwelling |
|---|
| 171 | real(jprb), intent(out), dimension(ncol, nlev+1) :: flux_up, flux_dn |
|---|
| 172 | |
|---|
| 173 | ! Albedo of the entire earth/atmosphere system below each half |
|---|
| 174 | ! level |
|---|
| 175 | real(jprb), dimension(ncol, nlev+1) :: albedo |
|---|
| 176 | |
|---|
| 177 | ! Upwelling radiation at each half-level due to emission below |
|---|
| 178 | ! that half-level (W m-2) |
|---|
| 179 | real(jprb), dimension(ncol, nlev+1) :: source |
|---|
| 180 | |
|---|
| 181 | ! Equal to 1/(1-albedo*reflectance) |
|---|
| 182 | real(jprb), dimension(ncol, nlev) :: inv_denominator |
|---|
| 183 | |
|---|
| 184 | ! Loop index for model level and column |
|---|
| 185 | integer :: jlev, jcol |
|---|
| 186 | |
|---|
| 187 | real(jprb) :: hook_handle |
|---|
| 188 | |
|---|
| 189 | if (lhook) call dr_hook('radiation_adding_ica_lw:fast_adding_ica_lw',0,hook_handle) |
|---|
| 190 | |
|---|
| 191 | ! Copy over downwelling fluxes above cloud from clear sky |
|---|
| 192 | flux_dn(:,1:i_cloud_top) = flux_dn_clear(:,1:i_cloud_top) |
|---|
| 193 | |
|---|
| 194 | albedo(:,nlev+1) = albedo_surf |
|---|
| 195 | |
|---|
| 196 | ! At the surface, the source is thermal emission |
|---|
| 197 | source(:,nlev+1) = emission_surf |
|---|
| 198 | |
|---|
| 199 | ! Work back up through the atmosphere and compute the albedo of |
|---|
| 200 | ! the entire earth/atmosphere system below that half-level, and |
|---|
| 201 | ! also the "source", which is the upwelling flux due to emission |
|---|
| 202 | ! below that level |
|---|
| 203 | do jlev = nlev,i_cloud_top,-1 |
|---|
| 204 | if (is_clear_sky_layer(jlev)) then |
|---|
| 205 | ! Reflectance of this layer is zero, simplifying the expression |
|---|
| 206 | do jcol = 1,ncol |
|---|
| 207 | albedo(jcol,jlev) = transmittance(jcol,jlev)*transmittance(jcol,jlev)*albedo(jcol,jlev+1) |
|---|
| 208 | source(jcol,jlev) = source_up(jcol,jlev) & |
|---|
| 209 | & + transmittance(jcol,jlev) * (source(jcol,jlev+1) & |
|---|
| 210 | & + albedo(jcol,jlev+1)*source_dn(jcol,jlev)) |
|---|
| 211 | end do |
|---|
| 212 | else |
|---|
| 213 | ! Loop over columns; explicit loop seems to be faster |
|---|
| 214 | do jcol = 1,ncol |
|---|
| 215 | ! Lacis and Hansen (1974) Eq 33, Shonk & Hogan (2008) Eq 10: |
|---|
| 216 | inv_denominator(jcol,jlev) = 1.0_jprb & |
|---|
| 217 | & / (1.0_jprb-albedo(jcol,jlev+1)*reflectance(jcol,jlev)) |
|---|
| 218 | ! Shonk & Hogan (2008) Eq 9, Petty (2006) Eq 13.81: |
|---|
| 219 | albedo(jcol,jlev) = reflectance(jcol,jlev) + transmittance(jcol,jlev)*transmittance(jcol,jlev) & |
|---|
| 220 | & * albedo(jcol,jlev+1) * inv_denominator(jcol,jlev) |
|---|
| 221 | ! Shonk & Hogan (2008) Eq 11: |
|---|
| 222 | source(jcol,jlev) = source_up(jcol,jlev) & |
|---|
| 223 | & + transmittance(jcol,jlev) * (source(jcol,jlev+1) & |
|---|
| 224 | & + albedo(jcol,jlev+1)*source_dn(jcol,jlev)) & |
|---|
| 225 | & * inv_denominator(jcol,jlev) |
|---|
| 226 | end do |
|---|
| 227 | end if |
|---|
| 228 | end do |
|---|
| 229 | |
|---|
| 230 | ! Compute the fluxes above the highest cloud |
|---|
| 231 | flux_up(:,i_cloud_top) = source(:,i_cloud_top) & |
|---|
| 232 | & + albedo(:,i_cloud_top)*flux_dn(:,i_cloud_top) |
|---|
| 233 | do jlev = i_cloud_top-1,1,-1 |
|---|
| 234 | flux_up(:,jlev) = transmittance(:,jlev)*flux_up(:,jlev+1) + source_up(:,jlev) |
|---|
| 235 | end do |
|---|
| 236 | |
|---|
| 237 | ! Work back down through the atmosphere from cloud top computing |
|---|
| 238 | ! the fluxes at each half-level |
|---|
| 239 | do jlev = i_cloud_top,nlev |
|---|
| 240 | if (is_clear_sky_layer(jlev)) then |
|---|
| 241 | do jcol = 1,ncol |
|---|
| 242 | flux_dn(jcol,jlev+1) = transmittance(jcol,jlev)*flux_dn(jcol,jlev) & |
|---|
| 243 | & + source_dn(jcol,jlev) |
|---|
| 244 | flux_up(jcol,jlev+1) = albedo(jcol,jlev+1)*flux_dn(jcol,jlev+1) & |
|---|
| 245 | & + source(jcol,jlev+1) |
|---|
| 246 | end do |
|---|
| 247 | else |
|---|
| 248 | do jcol = 1,ncol |
|---|
| 249 | ! Shonk & Hogan (2008) Eq 14 (after simplification): |
|---|
| 250 | flux_dn(jcol,jlev+1) & |
|---|
| 251 | & = (transmittance(jcol,jlev)*flux_dn(jcol,jlev) & |
|---|
| 252 | & + reflectance(jcol,jlev)*source(jcol,jlev+1) & |
|---|
| 253 | & + source_dn(jcol,jlev)) * inv_denominator(jcol,jlev) |
|---|
| 254 | ! Shonk & Hogan (2008) Eq 12: |
|---|
| 255 | flux_up(jcol,jlev+1) = albedo(jcol,jlev+1)*flux_dn(jcol,jlev+1) & |
|---|
| 256 | & + source(jcol,jlev+1) |
|---|
| 257 | end do |
|---|
| 258 | end if |
|---|
| 259 | end do |
|---|
| 260 | |
|---|
| 261 | if (lhook) call dr_hook('radiation_adding_ica_lw:fast_adding_ica_lw',1,hook_handle) |
|---|
| 262 | |
|---|
| 263 | end subroutine fast_adding_ica_lw |
|---|
| 264 | |
|---|
| 265 | |
|---|
| 266 | !--------------------------------------------------------------------- |
|---|
| 267 | ! If there is no scattering then fluxes may be computed simply by |
|---|
| 268 | ! passing down through the atmosphere computing the downwelling |
|---|
| 269 | ! fluxes from the transmission and emission of each layer, and then |
|---|
| 270 | ! passing back up through the atmosphere to compute the upwelling |
|---|
| 271 | ! fluxes in the same way. |
|---|
| 272 | subroutine calc_fluxes_no_scattering_lw(ncol, nlev, & |
|---|
| 273 | & transmittance, source_up, source_dn, emission_surf, albedo_surf, flux_up, flux_dn) |
|---|
| 274 | |
|---|
| 275 | use parkind1, only : jprb |
|---|
| 276 | use yomhook, only : lhook, dr_hook |
|---|
| 277 | |
|---|
| 278 | implicit none |
|---|
| 279 | |
|---|
| 280 | ! Inputs |
|---|
| 281 | integer, intent(in) :: ncol ! number of columns (may be spectral intervals) |
|---|
| 282 | integer, intent(in) :: nlev ! number of levels |
|---|
| 283 | |
|---|
| 284 | ! Surface emission (W m-2) and albedo |
|---|
| 285 | real(jprb), intent(in), dimension(ncol) :: emission_surf, albedo_surf |
|---|
| 286 | |
|---|
| 287 | ! Diffuse reflectance and transmittance of each layer |
|---|
| 288 | real(jprb), intent(in), dimension(ncol, nlev) :: transmittance |
|---|
| 289 | |
|---|
| 290 | ! Emission from each layer in an upward and downward direction |
|---|
| 291 | real(jprb), intent(in), dimension(ncol, nlev) :: source_up, source_dn |
|---|
| 292 | |
|---|
| 293 | ! Resulting fluxes (W m-2) at half-levels: diffuse upwelling and |
|---|
| 294 | ! downwelling |
|---|
| 295 | real(jprb), intent(out), dimension(ncol, nlev+1) :: flux_up, flux_dn |
|---|
| 296 | |
|---|
| 297 | ! Loop index for model level |
|---|
| 298 | integer :: jlev, jcol |
|---|
| 299 | |
|---|
| 300 | real(jprb) :: hook_handle |
|---|
| 301 | |
|---|
| 302 | if (lhook) call dr_hook('radiation_adding_ica_lw:calc_fluxes_no_scattering_lw',0,hook_handle) |
|---|
| 303 | |
|---|
| 304 | ! At top-of-atmosphere there is no diffuse downwelling radiation |
|---|
| 305 | flux_dn(:,1) = 0.0_jprb |
|---|
| 306 | |
|---|
| 307 | ! Work down through the atmosphere computing the downward fluxes |
|---|
| 308 | ! at each half-level |
|---|
| 309 | ! Added for DWD (2020) |
|---|
| 310 | !NEC$ outerloop_unroll(8) |
|---|
| 311 | do jlev = 1,nlev |
|---|
| 312 | do jcol = 1,ncol |
|---|
| 313 | flux_dn(jcol,jlev+1) = transmittance(jcol,jlev)*flux_dn(jcol,jlev) + source_dn(jcol,jlev) |
|---|
| 314 | end do |
|---|
| 315 | end do |
|---|
| 316 | |
|---|
| 317 | ! Surface reflection and emission |
|---|
| 318 | flux_up(:,nlev+1) = emission_surf + albedo_surf * flux_dn(:,nlev+1) |
|---|
| 319 | |
|---|
| 320 | ! Work back up through the atmosphere computing the upward fluxes |
|---|
| 321 | ! at each half-level |
|---|
| 322 | ! Added for DWD (2020) |
|---|
| 323 | !NEC$ outerloop_unroll(8) |
|---|
| 324 | do jlev = nlev,1,-1 |
|---|
| 325 | do jcol = 1,ncol |
|---|
| 326 | flux_up(jcol,jlev) = transmittance(jcol,jlev)*flux_up(jcol,jlev+1) + source_up(jcol,jlev) |
|---|
| 327 | end do |
|---|
| 328 | end do |
|---|
| 329 | |
|---|
| 330 | if (lhook) call dr_hook('radiation_adding_ica_lw:calc_fluxes_no_scattering_lw',1,hook_handle) |
|---|
| 331 | |
|---|
| 332 | end subroutine calc_fluxes_no_scattering_lw |
|---|
| 333 | |
|---|
| 334 | end module radiation_adding_ica_lw |
|---|
