| 1 | ! radiation_spartacus_sw.F90 - SPARTACUS shortwave solver |
|---|
| 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-04-11 R. Hogan Receive albedos at g-points |
|---|
| 17 | ! 2017-04-22 R. Hogan Store surface fluxes at all g-points |
|---|
| 18 | ! 2017-06-30 R. Hogan Reformulate to use total_albedo_direct not total_source |
|---|
| 19 | ! 2017-07-03 R. Hogan Explicit calculation of encroachment |
|---|
| 20 | ! 2017-10-23 R. Hogan Renamed single-character variables |
|---|
| 21 | ! 2018-02-20 R. Hogan Corrected "computed" encroachment |
|---|
| 22 | ! 2018-03-09 R. Hogan Security on computed encroachment, transmittance and reflectance |
|---|
| 23 | ! 2018-08-29 R. Hogan Reformulate horizontal migration distances in step_migrations |
|---|
| 24 | ! 2018-09-02 R. Hogan Bug fix in x_direct in step_migrations |
|---|
| 25 | ! 2018-09-03 R. Hogan Security via min_cloud_effective_size |
|---|
| 26 | ! 2018-09-04 R. Hogan Use encroachment_scaling and read encroachment edge_length from upper layer |
|---|
| 27 | ! 2018-09-13 R. Hogan Added "Fractal" encroachment option |
|---|
| 28 | ! 2018-09-14 R. Hogan Added "Zero" encroachment option |
|---|
| 29 | ! 2018-10-08 R. Hogan Call calc_region_properties |
|---|
| 30 | ! 2018-10-15 R. Hogan Added call to fast_expm_exchange instead of expm |
|---|
| 31 | ! 2019-01-12 R. Hogan Use inv_inhom_effective_size if allocated |
|---|
| 32 | ! 2019-02-10 R. Hogan Renamed "encroachment" to "entrapment" |
|---|
| 33 | |
|---|
| 34 | module radiation_spartacus_sw |
|---|
| 35 | |
|---|
| 36 | public |
|---|
| 37 | |
|---|
| 38 | contains |
|---|
| 39 | |
|---|
| 40 | ! Small routine for scaling cloud optical depth in the cloudy |
|---|
| 41 | ! regions |
|---|
| 42 | #include "radiation_optical_depth_scaling.h" |
|---|
| 43 | |
|---|
| 44 | ! This module contains just one exported subroutine, the shortwave |
|---|
| 45 | ! solver using the Speedy Algorithm for Radiative Transfer through |
|---|
| 46 | ! Cloud Sides (SPARTACUS), which can represent 3D effects using a |
|---|
| 47 | ! matrix form of the two-stream equations. |
|---|
| 48 | ! |
|---|
| 49 | ! Sections: |
|---|
| 50 | ! 1: Prepare general variables and arrays |
|---|
| 51 | ! 2: Prepare column-specific variables and arrays |
|---|
| 52 | ! 3: First loop over layers |
|---|
| 53 | ! 3.1: Layer-specific initialization |
|---|
| 54 | ! 3.2: Compute gamma variables |
|---|
| 55 | ! 3.2a: Clear-sky case |
|---|
| 56 | ! 3.2b: Cloudy case |
|---|
| 57 | ! 3.3: Compute reflection, transmission and sources |
|---|
| 58 | ! 3.3a: g-points with 3D effects |
|---|
| 59 | ! 3.3b: g-points without 3D effects |
|---|
| 60 | ! 4: Compute total albedos |
|---|
| 61 | ! 4.1: Adding method |
|---|
| 62 | ! 4.2: Overlap and entrapment |
|---|
| 63 | ! 5: Compute fluxes |
|---|
| 64 | subroutine solver_spartacus_sw(nlev,istartcol,iendcol, & |
|---|
| 65 | & config, single_level, thermodynamics, cloud, & |
|---|
| 66 | & od, ssa, g, od_cloud, ssa_cloud, g_cloud, & |
|---|
| 67 | & albedo_direct, albedo_diffuse, incoming_sw, & |
|---|
| 68 | & flux) |
|---|
| 69 | |
|---|
| 70 | use parkind1, only : jprb |
|---|
| 71 | use yomhook, only : lhook, dr_hook, jphook |
|---|
| 72 | |
|---|
| 73 | use radiation_io, only : nulout |
|---|
| 74 | use radiation_config, only : config_type, IPdfShapeGamma, & |
|---|
| 75 | & IEntrapmentZero, IEntrapmentEdgeOnly, IEntrapmentExplicit, & |
|---|
| 76 | & IEntrapmentExplicitNonFractal, IEntrapmentMaximum |
|---|
| 77 | use radiation_single_level, only : single_level_type |
|---|
| 78 | use radiation_thermodynamics, only : thermodynamics_type |
|---|
| 79 | use radiation_cloud, only : cloud_type |
|---|
| 80 | use radiation_regions, only : calc_region_properties |
|---|
| 81 | use radiation_overlap, only : calc_overlap_matrices |
|---|
| 82 | use radiation_flux, only : flux_type, & |
|---|
| 83 | & indexed_sum, add_indexed_sum |
|---|
| 84 | use radiation_matrix |
|---|
| 85 | use radiation_two_stream, only : calc_two_stream_gammas_sw, & |
|---|
| 86 | & calc_reflectance_transmittance_sw, calc_frac_scattered_diffuse_sw |
|---|
| 87 | use radiation_constants, only : Pi, GasConstantDryAir, & |
|---|
| 88 | & AccelDueToGravity |
|---|
| 89 | USE mod_phys_lmdz_para |
|---|
| 90 | |
|---|
| 91 | implicit none |
|---|
| 92 | |
|---|
| 93 | ! Inputs |
|---|
| 94 | integer, intent(in) :: nlev ! number of model levels |
|---|
| 95 | integer, intent(in) :: istartcol, iendcol ! range of columns to process |
|---|
| 96 | type(config_type), intent(in) :: config |
|---|
| 97 | type(single_level_type), intent(in) :: single_level |
|---|
| 98 | type(thermodynamics_type),intent(in) :: thermodynamics |
|---|
| 99 | type(cloud_type), intent(in) :: cloud |
|---|
| 100 | |
|---|
| 101 | ! Gas and aerosol optical depth, single-scattering albedo and |
|---|
| 102 | ! asymmetry factor at each shortwave g-point |
|---|
| 103 | real(jprb), intent(in), dimension(config%n_g_sw,nlev,istartcol:iendcol) :: & |
|---|
| 104 | & od, ssa, g |
|---|
| 105 | |
|---|
| 106 | ! Cloud and precipitation optical depth, single-scattering albedo and |
|---|
| 107 | ! asymmetry factor in each shortwave band |
|---|
| 108 | real(jprb), intent(in), dimension(config%n_bands_sw,nlev,istartcol:iendcol) :: & |
|---|
| 109 | & od_cloud, ssa_cloud, g_cloud |
|---|
| 110 | |
|---|
| 111 | ! Direct and diffuse surface albedos, and the incoming shortwave |
|---|
| 112 | ! flux into a plane perpendicular to the incoming radiation at |
|---|
| 113 | ! top-of-atmosphere in each of the shortwave g points |
|---|
| 114 | real(jprb), intent(in), dimension(config%n_g_sw,istartcol:iendcol) :: & |
|---|
| 115 | & albedo_direct, albedo_diffuse, incoming_sw |
|---|
| 116 | |
|---|
| 117 | ! Output |
|---|
| 118 | type(flux_type), intent(inout):: flux |
|---|
| 119 | |
|---|
| 120 | integer :: nreg, ng |
|---|
| 121 | integer :: nregactive ! =1 in clear layer, =nreg in a cloudy layer |
|---|
| 122 | integer :: jcol, jlev, jg, jreg, iband, jreg2,jreg3 |
|---|
| 123 | #ifdef EXPLICIT_EDGE_ENTRAPMENT |
|---|
| 124 | integer :: jreg4 |
|---|
| 125 | #endif |
|---|
| 126 | integer :: ng3D ! Number of g-points with small enough gas optical |
|---|
| 127 | ! depth that 3D effects need to be represented |
|---|
| 128 | |
|---|
| 129 | ! Ratio of gas constant for dry air to acceleration due to gravity |
|---|
| 130 | real(jprb), parameter :: R_over_g = GasConstantDryAir / AccelDueToGravity |
|---|
| 131 | |
|---|
| 132 | real(jprb) :: mu0, one_over_mu0, tan_sza, transfer_scaling |
|---|
| 133 | |
|---|
| 134 | ! The tangent of the effective zenith angle for diffuse radiation |
|---|
| 135 | ! that is appropriate for 3D transport |
|---|
| 136 | real(jprb), parameter :: tan_diffuse_angle_3d = Pi * 0.5_jprb |
|---|
| 137 | |
|---|
| 138 | ! The minimum cosine of solar zenith angle to consider for 3D |
|---|
| 139 | ! effects, equivalent to one solar radius (0.2615 degrees) |
|---|
| 140 | real(jprb), parameter :: min_mu0_3d = 0.004625_jprb |
|---|
| 141 | |
|---|
| 142 | ! Optical depth, single scattering albedo and asymmetry factor in |
|---|
| 143 | ! each region and (except for asymmetry) at each g-point |
|---|
| 144 | real(jprb), dimension(config%n_g_sw, config%nregions) & |
|---|
| 145 | & :: od_region, ssa_region |
|---|
| 146 | real(jprb), dimension(config%nregions) :: g_region |
|---|
| 147 | |
|---|
| 148 | ! Scattering optical depths of gases and clouds |
|---|
| 149 | real(jprb) :: scat_od, scat_od_cloud |
|---|
| 150 | |
|---|
| 151 | ! The area fractions of each region |
|---|
| 152 | real(jprb) :: region_fracs(1:config%nregions,nlev,istartcol:iendcol) |
|---|
| 153 | |
|---|
| 154 | ! The scaling used for the optical depth in the cloudy regions |
|---|
| 155 | real(jprb) :: od_scaling(2:config%nregions,nlev,istartcol:iendcol) |
|---|
| 156 | |
|---|
| 157 | ! The length of the interface between regions jreg and jreg+1 per |
|---|
| 158 | ! unit area of gridbox, equal to L_diff^ab in Hogan and Shonk |
|---|
| 159 | ! (2013). When the first index is 3, this is the length of the |
|---|
| 160 | ! interface between regions 3 and 1. This is actually the |
|---|
| 161 | ! effective length oriented to a photon with random azimuth angle, |
|---|
| 162 | ! so is the true edge length divided by pi. |
|---|
| 163 | real(jprb) :: edge_length(3,nlev) |
|---|
| 164 | |
|---|
| 165 | ! Element i,j gives the rate of 3D transfer of diffuse/direct |
|---|
| 166 | ! radiation from region i to region j, multiplied by the thickness |
|---|
| 167 | ! of the layer in m |
|---|
| 168 | real(jprb) :: transfer_rate_diffuse(config%nregions,config%nregions) |
|---|
| 169 | real(jprb) :: transfer_rate_direct(config%nregions,config%nregions) |
|---|
| 170 | |
|---|
| 171 | ! Directional overlap matrices defined at all layer interfaces |
|---|
| 172 | ! including top-of-atmosphere and the surface |
|---|
| 173 | real(jprb), dimension(config%nregions,config%nregions,nlev+1, & |
|---|
| 174 | & istartcol:iendcol) :: u_matrix, v_matrix |
|---|
| 175 | |
|---|
| 176 | ! Two-stream variables |
|---|
| 177 | real(jprb), dimension(config%n_g_sw, config%nregions) & |
|---|
| 178 | & :: gamma1, gamma2, gamma3 |
|---|
| 179 | |
|---|
| 180 | ! Matrix Gamma multiplied by the layer thickness z1, so units |
|---|
| 181 | ! metres. After calling expm, this array contains the matrix |
|---|
| 182 | ! exponential of the original. |
|---|
| 183 | real(jprb) :: Gamma_z1(config%n_g_sw,3*config%nregions,3*config%nregions) |
|---|
| 184 | |
|---|
| 185 | ! Diffuse reflection and transmission matrices of each layer |
|---|
| 186 | real(jprb), dimension(config%n_g_sw, config%nregions, & |
|---|
| 187 | & config%nregions, nlev) :: reflectance, transmittance |
|---|
| 188 | |
|---|
| 189 | ! Clear-sky diffuse reflection and transmission matrices of each |
|---|
| 190 | ! layer |
|---|
| 191 | real(jprb), dimension(config%n_g_sw, nlev) :: ref_clear, trans_clear |
|---|
| 192 | |
|---|
| 193 | ! Matrices translating the direct flux entering the layer from |
|---|
| 194 | ! above to the reflected radiation exiting upwards (ref_dir) and |
|---|
| 195 | ! the scattered radiation exiting downwards (trans_dir_diff), |
|---|
| 196 | ! along with the direct unscattered transmission matrix |
|---|
| 197 | ! (trans_dir_dir). |
|---|
| 198 | real(jprb), dimension(config%n_g_sw, config%nregions, config%nregions, nlev) & |
|---|
| 199 | & :: ref_dir, trans_dir_diff, trans_dir_dir |
|---|
| 200 | ! ...clear-sky equivalents |
|---|
| 201 | real(jprb), dimension(config%n_g_sw, nlev) & |
|---|
| 202 | & :: ref_dir_clear, trans_dir_diff_clear, trans_dir_dir_clear |
|---|
| 203 | |
|---|
| 204 | ! The fluxes downwelling from the bottom of the layer due to |
|---|
| 205 | ! scattering by the direct beam within the layer |
|---|
| 206 | real(jprb), dimension(config%n_g_sw, config%nregions) :: source_dn |
|---|
| 207 | ! ...clear-sky equivalent |
|---|
| 208 | real(jprb), dimension(config%n_g_sw) :: source_dn_clear |
|---|
| 209 | |
|---|
| 210 | ! The fluxes upwelling just above the base of a layer due to |
|---|
| 211 | ! reflection of the direct downwelling beam; this is just used as |
|---|
| 212 | ! a temporary variable |
|---|
| 213 | real(jprb), dimension(config%n_g_sw, config%nregions) :: total_source |
|---|
| 214 | |
|---|
| 215 | ! Direct downwelling flux below and above an interface between |
|---|
| 216 | ! layers into a plane perpendicular to the direction of the sun |
|---|
| 217 | real(jprb), dimension(config%n_g_sw, config%nregions) & |
|---|
| 218 | & :: direct_dn_below, direct_dn_above |
|---|
| 219 | ! ...clear-sky equivalent (no distinction between "above/below") |
|---|
| 220 | real(jprb), dimension(config%n_g_sw) :: direct_dn_clear |
|---|
| 221 | |
|---|
| 222 | ! Total albedo of the atmosphere/surface just above a layer |
|---|
| 223 | ! interface with respect to downwelling diffuse and direct |
|---|
| 224 | ! radiation at that interface, where level index = 1 corresponds |
|---|
| 225 | ! to the top-of-atmosphere |
|---|
| 226 | real(jprb), dimension(config%n_g_sw, config%nregions, & |
|---|
| 227 | & config%nregions, nlev+1) :: total_albedo, total_albedo_direct |
|---|
| 228 | ! ...clear-sky equivalent |
|---|
| 229 | real(jprb), dimension(config%n_g_sw, nlev+1) & |
|---|
| 230 | & :: total_albedo_clear, total_albedo_clear_direct |
|---|
| 231 | |
|---|
| 232 | ! As total_albedo, but just below a layer interface |
|---|
| 233 | real(jprb), dimension(config%n_g_sw, config%nregions, config%nregions) & |
|---|
| 234 | & :: total_albedo_below, total_albedo_below_direct |
|---|
| 235 | |
|---|
| 236 | ! Temporary array for applying adding method with entrapment to |
|---|
| 237 | ! albedo matrices |
|---|
| 238 | real(jprb), dimension(config%n_g_sw, config%nregions, config%nregions) & |
|---|
| 239 | & :: albedo_part |
|---|
| 240 | |
|---|
| 241 | ! Horizontal migration distance (m) of reflected light |
|---|
| 242 | real(jprb), dimension(config%n_g_sw, config%nregions) & |
|---|
| 243 | & :: x_diffuse, x_direct |
|---|
| 244 | ! Temporary variables when applying overlap rules |
|---|
| 245 | real(jprb), dimension(config%n_g_sw, config%nregions) & |
|---|
| 246 | & :: x_diffuse_above, x_direct_above |
|---|
| 247 | |
|---|
| 248 | real(jprb), dimension(config%n_g_sw, config%nregions, config%nregions) & |
|---|
| 249 | & :: entrapment |
|---|
| 250 | |
|---|
| 251 | ! The following is used to store matrices of the form I-A*B that |
|---|
| 252 | ! are used on the denominator of some expressions |
|---|
| 253 | real(jprb) :: denominator(config%n_g_sw,config%nregions,config%nregions) |
|---|
| 254 | |
|---|
| 255 | ! Clear-sky equivalent, but actually its reciprocal to replace |
|---|
| 256 | ! some divisions by multiplications |
|---|
| 257 | real(jprb), dimension(config%n_g_sw) :: inv_denom_scalar |
|---|
| 258 | |
|---|
| 259 | ! Final step in working out how much transport between regions |
|---|
| 260 | ! above occurs |
|---|
| 261 | real(jprb), dimension(config%n_g_sw) :: fractal_factor |
|---|
| 262 | |
|---|
| 263 | ! Inverse of cloud effective size (m^-1) |
|---|
| 264 | real(jprb) :: inv_effective_size |
|---|
| 265 | |
|---|
| 266 | ! Layer depth (m) |
|---|
| 267 | real(jprb) :: dz, layer_depth(nlev) |
|---|
| 268 | |
|---|
| 269 | ! Upwelling and downwelling fluxes above/below layer interfaces |
|---|
| 270 | real(jprb), dimension(config%n_g_sw, config%nregions) & |
|---|
| 271 | & :: flux_up_above, flux_dn_above, flux_dn_below |
|---|
| 272 | ! Clear-sky upwelling and downwelling fluxes (which don't |
|---|
| 273 | ! distinguish between whether they are above/below a layer |
|---|
| 274 | ! interface) |
|---|
| 275 | real(jprb), dimension(config%n_g_sw) :: flux_up_clear, flux_dn_clear |
|---|
| 276 | |
|---|
| 277 | ! Index of top-most cloudy layer, or nlev+1 if no cloud |
|---|
| 278 | integer :: i_cloud_top |
|---|
| 279 | |
|---|
| 280 | ! Keep a count of the number of calls to the two ways of computing |
|---|
| 281 | ! reflectance/transmittance matrices |
|---|
| 282 | integer :: n_calls_expm, n_calls_meador_weaver |
|---|
| 283 | |
|---|
| 284 | ! Identify clear-sky layers, with pseudo layers for outer space |
|---|
| 285 | ! and below the ground, both treated as single-region clear skies |
|---|
| 286 | logical :: is_clear_sky_layer(0:nlev+1) |
|---|
| 287 | |
|---|
| 288 | ! Used in computing rates of lateral radiation transfer |
|---|
| 289 | real(jprb), parameter :: four_over_pi = 4.0_jprb / Pi |
|---|
| 290 | |
|---|
| 291 | ! Maximum entrapment coefficient |
|---|
| 292 | real(jprb) :: max_entr |
|---|
| 293 | |
|---|
| 294 | real(jphook) :: hook_handle |
|---|
| 295 | |
|---|
| 296 | if (lhook) call dr_hook('radiation_spartacus_sw:solver_spartacus_sw',0,hook_handle) |
|---|
| 297 | |
|---|
| 298 | ! -------------------------------------------------------- |
|---|
| 299 | ! Section 1: Prepare general variables and arrays |
|---|
| 300 | ! -------------------------------------------------------- |
|---|
| 301 | |
|---|
| 302 | ! Copy array dimensions to local variables for convenience |
|---|
| 303 | nreg = config%nregions |
|---|
| 304 | ng = config%n_g_sw |
|---|
| 305 | |
|---|
| 306 | ! Reset count of number of calls to the two ways to compute |
|---|
| 307 | ! reflectance/transmission matrices |
|---|
| 308 | n_calls_expm = 0 |
|---|
| 309 | n_calls_meador_weaver = 0 |
|---|
| 310 | |
|---|
| 311 | ! Compute the wavelength-independent region fractions and |
|---|
| 312 | ! optical-depth scalings |
|---|
| 313 | call calc_region_properties(nlev,nreg,istartcol,iendcol, & |
|---|
| 314 | & config%i_cloud_pdf_shape == IPdfShapeGamma, & |
|---|
| 315 | & cloud%fraction, cloud%fractional_std, region_fracs, & |
|---|
| 316 | & od_scaling, config%cloud_fraction_threshold) |
|---|
| 317 | |
|---|
| 318 | ! Compute wavelength-independent overlap matrices u_matrix and v_matrix |
|---|
| 319 | call calc_overlap_matrices(nlev,nreg,istartcol,iendcol, & |
|---|
| 320 | & region_fracs, cloud%overlap_param, & |
|---|
| 321 | & u_matrix, v_matrix, decorrelation_scaling=config%cloud_inhom_decorr_scaling, & |
|---|
| 322 | & cloud_fraction_threshold=config%cloud_fraction_threshold, & |
|---|
| 323 | & use_beta_overlap=config%use_beta_overlap, & |
|---|
| 324 | & cloud_cover=flux%cloud_cover_sw) |
|---|
| 325 | |
|---|
| 326 | if (config%iverbose >= 3) then |
|---|
| 327 | write(nulout,'(a)',advance='no') ' Processing columns' |
|---|
| 328 | end if |
|---|
| 329 | |
|---|
| 330 | print*,'Dans radiation_spartacus, mpi_rank, omp_rank, & |
|---|
| 331 | size, chape inv_cloud = ',& |
|---|
| 332 | mpi_rank, omp_rank, & |
|---|
| 333 | shape(cloud%inv_cloud_effective_size), & |
|---|
| 334 | size(cloud%inv_cloud_effective_size) |
|---|
| 335 | |
|---|
| 336 | ! Main loop over columns |
|---|
| 337 | do jcol = istartcol, iendcol |
|---|
| 338 | ! -------------------------------------------------------- |
|---|
| 339 | ! Section 2: Prepare column-specific variables and arrays |
|---|
| 340 | ! -------------------------------------------------------- |
|---|
| 341 | |
|---|
| 342 | if (config%iverbose >= 3) then |
|---|
| 343 | write(nulout,'(a)',advance='no') '.' |
|---|
| 344 | end if |
|---|
| 345 | |
|---|
| 346 | ! Copy local cosine of the solar zenith angle |
|---|
| 347 | mu0 = single_level%cos_sza(jcol) |
|---|
| 348 | |
|---|
| 349 | ! Skip profile if sun is too low in the sky |
|---|
| 350 | if (mu0 < 1.0e-10_jprb) then |
|---|
| 351 | flux%sw_dn(jcol,:) = 0.0_jprb |
|---|
| 352 | flux%sw_up(jcol,:) = 0.0_jprb |
|---|
| 353 | if (allocated(flux%sw_dn_direct)) then |
|---|
| 354 | flux%sw_dn_direct(jcol,:) = 0.0_jprb |
|---|
| 355 | end if |
|---|
| 356 | if (config%do_clear) then |
|---|
| 357 | flux%sw_dn_clear(jcol,:) = 0.0_jprb |
|---|
| 358 | flux%sw_up_clear(jcol,:) = 0.0_jprb |
|---|
| 359 | if (allocated(flux%sw_dn_direct_clear)) then |
|---|
| 360 | flux%sw_dn_direct_clear(jcol,:) = 0.0_jprb |
|---|
| 361 | end if |
|---|
| 362 | end if |
|---|
| 363 | |
|---|
| 364 | if (config%do_save_spectral_flux) then |
|---|
| 365 | flux%sw_dn_band(:,jcol,:) = 0.0_jprb |
|---|
| 366 | flux%sw_up_band(:,jcol,:) = 0.0_jprb |
|---|
| 367 | if (allocated(flux%sw_dn_direct_band)) then |
|---|
| 368 | flux%sw_dn_direct_band(:,jcol,:) = 0.0_jprb |
|---|
| 369 | end if |
|---|
| 370 | if (config%do_clear) then |
|---|
| 371 | flux%sw_dn_clear_band(:,jcol,:) = 0.0_jprb |
|---|
| 372 | flux%sw_up_clear_band(:,jcol,:) = 0.0_jprb |
|---|
| 373 | if (allocated(flux%sw_dn_direct_clear_band)) then |
|---|
| 374 | flux%sw_dn_direct_clear_band(:,jcol,:) = 0.0_jprb |
|---|
| 375 | end if |
|---|
| 376 | end if |
|---|
| 377 | end if |
|---|
| 378 | |
|---|
| 379 | flux%sw_dn_diffuse_surf_g(:,jcol) = 0.0_jprb |
|---|
| 380 | flux%sw_dn_direct_surf_g(:,jcol) = 0.0_jprb |
|---|
| 381 | if (config%do_clear) then |
|---|
| 382 | flux%sw_dn_diffuse_surf_clear_g(:,jcol) = 0.0_jprb |
|---|
| 383 | flux%sw_dn_direct_surf_clear_g(:,jcol) = 0.0_jprb |
|---|
| 384 | end if |
|---|
| 385 | |
|---|
| 386 | cycle |
|---|
| 387 | end if ! sun is below the horizon |
|---|
| 388 | |
|---|
| 389 | ! At this point mu0 >= 1.0e-10 |
|---|
| 390 | |
|---|
| 391 | ! Used to compute rate of attenuation of direct solar beam |
|---|
| 392 | ! through the atmosphere |
|---|
| 393 | one_over_mu0 = 1.0_jprb / mu0 |
|---|
| 394 | |
|---|
| 395 | ! The rate at which direct radiation enters cloud sides is |
|---|
| 396 | ! proportional to the tangent of the solar zenith angle, but |
|---|
| 397 | ! this gets very large when the sun is low in the sky, in which |
|---|
| 398 | ! case we limit it to one solar radius above the horizon |
|---|
| 399 | if (mu0 < min_mu0_3d) then |
|---|
| 400 | tan_sza = sqrt(1.0_jprb/(min_mu0_3d*min_mu0_3d) - 1.0_jprb) |
|---|
| 401 | else if (one_over_mu0 > 1.0_jprb) then |
|---|
| 402 | tan_sza = sqrt(one_over_mu0*one_over_mu0 - 1.0_jprb & |
|---|
| 403 | & + config%overhead_sun_factor) |
|---|
| 404 | else |
|---|
| 405 | ! Just in case we get mu0 > 1... |
|---|
| 406 | tan_sza = sqrt(config%overhead_sun_factor) |
|---|
| 407 | end if |
|---|
| 408 | |
|---|
| 409 | ! Define which layers contain cloud; assume that |
|---|
| 410 | ! cloud%crop_cloud_fraction has already been called |
|---|
| 411 | is_clear_sky_layer = .true. |
|---|
| 412 | i_cloud_top = nlev+1 |
|---|
| 413 | do jlev = nlev,1,-1 |
|---|
| 414 | if (cloud%fraction(jcol,jlev) > 0.0_jprb) then |
|---|
| 415 | is_clear_sky_layer(jlev) = .false. |
|---|
| 416 | i_cloud_top = jlev |
|---|
| 417 | end if |
|---|
| 418 | end do |
|---|
| 419 | |
|---|
| 420 | ! -------------------------------------------------------- |
|---|
| 421 | ! Section 3: First loop over layers |
|---|
| 422 | ! -------------------------------------------------------- |
|---|
| 423 | ! In this section the reflectance, transmittance and sources |
|---|
| 424 | ! are computed for each layer |
|---|
| 425 | do jlev = 1,nlev ! Start at top-of-atmosphere |
|---|
| 426 | ! -------------------------------------------------------- |
|---|
| 427 | ! Section 3.1: Layer-specific initialization |
|---|
| 428 | ! -------------------------------------------------------- |
|---|
| 429 | |
|---|
| 430 | ! Array-wise assignments |
|---|
| 431 | gamma1 = 0.0_jprb |
|---|
| 432 | gamma2 = 0.0_jprb |
|---|
| 433 | gamma3 = 0.0_jprb |
|---|
| 434 | Gamma_z1= 0.0_jprb |
|---|
| 435 | transfer_rate_direct(:,:) = 0.0_jprb |
|---|
| 436 | transfer_rate_diffuse(:,:) = 0.0_jprb |
|---|
| 437 | edge_length(:,jlev) = 0.0_jprb |
|---|
| 438 | |
|---|
| 439 | ! The following is from the hydrostatic equation |
|---|
| 440 | ! and ideal gas law: dz = dp * R * T / (p * g) |
|---|
| 441 | layer_depth(jlev) = R_over_g & |
|---|
| 442 | & * (thermodynamics%pressure_hl(jcol,jlev+1) & |
|---|
| 443 | & - thermodynamics%pressure_hl(jcol,jlev)) & |
|---|
| 444 | & * (thermodynamics%temperature_hl(jcol,jlev) & |
|---|
| 445 | & + thermodynamics%temperature_hl(jcol,jlev+1)) & |
|---|
| 446 | & / (thermodynamics%pressure_hl(jcol,jlev) & |
|---|
| 447 | & + thermodynamics%pressure_hl(jcol,jlev+1)) |
|---|
| 448 | |
|---|
| 449 | ! -------------------------------------------------------- |
|---|
| 450 | ! Section 3.2: Compute gamma variables |
|---|
| 451 | ! -------------------------------------------------------- |
|---|
| 452 | if (is_clear_sky_layer(jlev)) then |
|---|
| 453 | ! --- Section 3.2a: Clear-sky case -------------------- |
|---|
| 454 | |
|---|
| 455 | nregactive = 1 ! Only region 1 (clear-sky) is active |
|---|
| 456 | |
|---|
| 457 | ! Copy optical depth and single-scattering albedo of |
|---|
| 458 | ! clear-sky region |
|---|
| 459 | od_region(1:ng,1) = od(1:ng,jlev,jcol) |
|---|
| 460 | ssa_region(1:ng,1) = ssa(1:ng,jlev,jcol) |
|---|
| 461 | |
|---|
| 462 | ! Compute two-stream variables gamma1-gamma3 (gamma4=1-gamma3) |
|---|
| 463 | call calc_two_stream_gammas_sw(ng, & |
|---|
| 464 | & mu0, ssa(1:ng,jlev,jcol), g(1:ng,jlev,jcol), & |
|---|
| 465 | & gamma1(1:ng,1), gamma2(1:ng,1), gamma3(1:ng,1)) |
|---|
| 466 | |
|---|
| 467 | if (config%use_expm_everywhere) then |
|---|
| 468 | ! Treat cloud-free layer as 3D: the matrix exponential |
|---|
| 469 | ! "expm" is used as much as possible (provided the |
|---|
| 470 | ! optical depth is not too high). ng3D is initially |
|---|
| 471 | ! set to the total number of g-points, then the |
|---|
| 472 | ! clear-sky optical depths are searched until a value |
|---|
| 473 | ! exceeds the threshold for treating 3D effects, and |
|---|
| 474 | ! if found it is reset to that. Note that we are |
|---|
| 475 | ! assuming that the g-points have been reordered in |
|---|
| 476 | ! approximate order of gas optical depth. |
|---|
| 477 | ng3D = ng |
|---|
| 478 | do jg = 1, ng |
|---|
| 479 | if (od_region(jg,1) > config%max_gas_od_3D) then |
|---|
| 480 | ng3D = jg-1 |
|---|
| 481 | exit |
|---|
| 482 | end if |
|---|
| 483 | end do |
|---|
| 484 | else |
|---|
| 485 | ! Otherwise treat cloud-free layer using the classical |
|---|
| 486 | ! Meador & Weaver formulae for all g-points |
|---|
| 487 | ng3D = 0 |
|---|
| 488 | end if |
|---|
| 489 | |
|---|
| 490 | |
|---|
| 491 | else |
|---|
| 492 | ! --- Section 3.2b: Cloudy case ----------------------- |
|---|
| 493 | |
|---|
| 494 | ! Default number of g-points to treat with |
|---|
| 495 | ! matrix-exponential scheme |
|---|
| 496 | if (config%use_expm_everywhere) then |
|---|
| 497 | ng3D = ng ! All g-points |
|---|
| 498 | else |
|---|
| 499 | ng3D = 0 ! No g-points |
|---|
| 500 | end if |
|---|
| 501 | |
|---|
| 502 | if (config%do_3d_effects .and. & |
|---|
| 503 | & allocated(cloud%inv_cloud_effective_size) .and. & |
|---|
| 504 | & .not. (nreg == 2 .and. cloud%fraction(jcol,jlev) & |
|---|
| 505 | & > 1.0-config%cloud_fraction_threshold)) then |
|---|
| 506 | ! print*,' Dans radiation_spartacus mpi_rank, omp_rank, jcol, jlev, & |
|---|
| 507 | ! & cloud%inv_cloud_effective_size =', mpi_rank, & |
|---|
| 508 | ! & omp_rank, jcol, jlev, & |
|---|
| 509 | ! & cloud%inv_cloud_effective_size(jcol,jlev) |
|---|
| 510 | if (cloud%inv_cloud_effective_size(jcol,jlev) > 0.0_jprb) then |
|---|
| 511 | ! 3D effects are only simulated if |
|---|
| 512 | ! inv_cloud_effective_size is defined and greater |
|---|
| 513 | ! than zero |
|---|
| 514 | ng3D = ng |
|---|
| 515 | |
|---|
| 516 | ! Depth of current layer |
|---|
| 517 | dz = layer_depth(jlev) |
|---|
| 518 | |
|---|
| 519 | ! Compute cloud edge length per unit area of gridbox |
|---|
| 520 | ! from rearranging Hogan & Shonk (2013) Eq. 45, but |
|---|
| 521 | ! adding a factor of (1-frac) so that a region that |
|---|
| 522 | ! fully occupies the gridbox (frac=1) has an edge of |
|---|
| 523 | ! zero. We can therefore use the fraction of clear sky, |
|---|
| 524 | ! region_fracs(1,jlev,jcol) for convenience instead. The |
|---|
| 525 | ! pi on the denominator means that this is actually edge |
|---|
| 526 | ! length with respect to a light ray with a random |
|---|
| 527 | ! azimuthal direction. |
|---|
| 528 | edge_length(1,jlev) = four_over_pi & |
|---|
| 529 | & * region_fracs(1,jlev,jcol)*(1.0_jprb-region_fracs(1,jlev,jcol)) & |
|---|
| 530 | & * min(cloud%inv_cloud_effective_size(jcol,jlev), & |
|---|
| 531 | & 1.0_jprb / config%min_cloud_effective_size) |
|---|
| 532 | if (nreg > 2) then |
|---|
| 533 | ! The corresponding edge length between the two cloudy |
|---|
| 534 | ! regions is computed in the same way but treating the |
|---|
| 535 | ! optically denser of the two regions (region 3) as |
|---|
| 536 | ! the cloud; note that the fraction of this region may |
|---|
| 537 | ! be less than that of the optically less dense of the |
|---|
| 538 | ! two regions (region 2). For increased flexibility, |
|---|
| 539 | ! the user may specify the effective size of |
|---|
| 540 | ! inhomogeneities separately from the cloud effective |
|---|
| 541 | ! size. |
|---|
| 542 | if (allocated(cloud%inv_inhom_effective_size)) then |
|---|
| 543 | edge_length(2,jlev) = four_over_pi & |
|---|
| 544 | & * region_fracs(3,jlev,jcol)*(1.0_jprb-region_fracs(3,jlev,jcol)) & |
|---|
| 545 | & * min(cloud%inv_inhom_effective_size(jcol,jlev), & |
|---|
| 546 | & 1.0_jprb / config%min_cloud_effective_size) |
|---|
| 547 | else |
|---|
| 548 | edge_length(2,jlev) = four_over_pi & |
|---|
| 549 | & * region_fracs(3,jlev,jcol)*(1.0_jprb-region_fracs(3,jlev,jcol)) & |
|---|
| 550 | & * min(cloud%inv_cloud_effective_size(jcol,jlev), & |
|---|
| 551 | & 1.0_jprb / config%min_cloud_effective_size) |
|---|
| 552 | end if |
|---|
| 553 | |
|---|
| 554 | ! In the case of three regions, some of the cloud |
|---|
| 555 | ! boundary may go directly to the "thick" region |
|---|
| 556 | if (config%clear_to_thick_fraction > 0.0_jprb) then |
|---|
| 557 | edge_length(3,jlev) = config%clear_to_thick_fraction & |
|---|
| 558 | & * min(edge_length(1,jlev), edge_length(2,jlev)) |
|---|
| 559 | edge_length(1,jlev) = edge_length(1,jlev) - edge_length(3,jlev) |
|---|
| 560 | edge_length(2,jlev) = edge_length(2,jlev) - edge_length(3,jlev) |
|---|
| 561 | else |
|---|
| 562 | edge_length(3,jlev) = 0.0_jprb |
|---|
| 563 | end if |
|---|
| 564 | end if |
|---|
| 565 | |
|---|
| 566 | do jreg = 1, nreg-1 |
|---|
| 567 | ! Compute lateral transfer rates from region jreg to |
|---|
| 568 | ! jreg+1 following Hogan & Shonk (2013) Eq. 47, but |
|---|
| 569 | ! multiplied by dz because the transfer rate is |
|---|
| 570 | ! vertically integrated across the depth of the layer |
|---|
| 571 | if (region_fracs(jreg,jlev,jcol) > epsilon(1.0_jprb)) then |
|---|
| 572 | transfer_rate_direct(jreg,jreg+1) = dz & |
|---|
| 573 | & * edge_length(jreg,jlev) * tan_sza / region_fracs(jreg,jlev,jcol) |
|---|
| 574 | transfer_rate_diffuse(jreg,jreg+1) = dz & |
|---|
| 575 | & * edge_length(jreg,jlev) & |
|---|
| 576 | & * tan_diffuse_angle_3d / region_fracs(jreg,jlev,jcol) |
|---|
| 577 | end if |
|---|
| 578 | ! Compute transfer rates from region jreg+1 to |
|---|
| 579 | ! jreg |
|---|
| 580 | if (region_fracs(jreg+1,jlev,jcol) > epsilon(1.0_jprb)) then |
|---|
| 581 | transfer_rate_direct(jreg+1,jreg) = dz & |
|---|
| 582 | & * edge_length(jreg,jlev) & |
|---|
| 583 | & * tan_sza / region_fracs(jreg+1,jlev,jcol) |
|---|
| 584 | transfer_rate_diffuse(jreg+1,jreg) = dz & |
|---|
| 585 | & * edge_length(jreg,jlev) & |
|---|
| 586 | & * tan_diffuse_angle_3d / region_fracs(jreg+1,jlev,jcol) |
|---|
| 587 | end if |
|---|
| 588 | end do |
|---|
| 589 | |
|---|
| 590 | ! Compute transfer rates directly between regions 1 and |
|---|
| 591 | ! 3 |
|---|
| 592 | if (edge_length(3,jlev) > 0.0_jprb) then |
|---|
| 593 | if (region_fracs(1,jlev,jcol) > epsilon(1.0_jprb)) then |
|---|
| 594 | transfer_rate_direct(1,3) = dz & |
|---|
| 595 | & * edge_length(3,jlev) * tan_sza / region_fracs(1,jlev,jcol) |
|---|
| 596 | transfer_rate_diffuse(1,3) = dz & |
|---|
| 597 | & * edge_length(3,jlev) & |
|---|
| 598 | & * tan_diffuse_angle_3d / region_fracs(1,jlev,jcol) |
|---|
| 599 | end if |
|---|
| 600 | if (region_fracs(3,jlev,jcol) > epsilon(1.0_jprb)) then |
|---|
| 601 | transfer_rate_direct(3,1) = dz & |
|---|
| 602 | & * edge_length(3,jlev) * tan_sza / region_fracs(3,jlev,jcol) |
|---|
| 603 | transfer_rate_diffuse(3,1) = dz & |
|---|
| 604 | & * edge_length(3,jlev) & |
|---|
| 605 | & * tan_diffuse_angle_3d / region_fracs(3,jlev,jcol) |
|---|
| 606 | end if |
|---|
| 607 | end if |
|---|
| 608 | |
|---|
| 609 | ! Don't allow the transfer rate out of a region to be |
|---|
| 610 | ! equivalent to a loss of exp(-10) through the layer |
|---|
| 611 | where (transfer_rate_direct > config%max_3d_transfer_rate) |
|---|
| 612 | transfer_rate_direct = config%max_3d_transfer_rate |
|---|
| 613 | end where |
|---|
| 614 | where (transfer_rate_diffuse > config%max_3d_transfer_rate) |
|---|
| 615 | transfer_rate_diffuse = config%max_3d_transfer_rate |
|---|
| 616 | end where |
|---|
| 617 | |
|---|
| 618 | end if ! Cloud has edge length required for 3D effects |
|---|
| 619 | end if ! Include 3D effects |
|---|
| 620 | |
|---|
| 621 | ! In a cloudy layer the number of active regions equals |
|---|
| 622 | ! the number of regions |
|---|
| 623 | nregactive = nreg |
|---|
| 624 | |
|---|
| 625 | ! Compute scattering properties of the regions at each |
|---|
| 626 | ! g-point, mapping from the cloud properties |
|---|
| 627 | ! defined in each band. |
|---|
| 628 | do jg = 1,ng |
|---|
| 629 | ! Mapping from g-point to band |
|---|
| 630 | iband = config%i_band_from_reordered_g_sw(jg) |
|---|
| 631 | |
|---|
| 632 | ! Scattering optical depth of clear-sky region |
|---|
| 633 | scat_od = od(jg,jlev,jcol)*ssa(jg,jlev,jcol) |
|---|
| 634 | |
|---|
| 635 | ! Scattering properties of clear-sky regions copied |
|---|
| 636 | ! over |
|---|
| 637 | od_region(jg,1) = od(jg, jlev, jcol) |
|---|
| 638 | ssa_region(jg,1) = ssa(jg, jlev, jcol) |
|---|
| 639 | g_region(1) = g(jg, jlev, jcol) |
|---|
| 640 | |
|---|
| 641 | ! Loop over cloudy regions |
|---|
| 642 | do jreg = 2,nreg |
|---|
| 643 | scat_od_cloud = od_cloud(iband,jlev,jcol) & |
|---|
| 644 | & * ssa_cloud(iband,jlev,jcol)*od_scaling(jreg,jlev,jcol) |
|---|
| 645 | ! Add scaled cloud optical depth to clear-sky value |
|---|
| 646 | od_region(jg,jreg) = od(jg,jlev,jcol) & |
|---|
| 647 | & + od_cloud(iband,jlev,jcol)*od_scaling(jreg,jlev,jcol) |
|---|
| 648 | ! Compute single-scattering albedo and asymmetry |
|---|
| 649 | ! factor of gas-cloud combination |
|---|
| 650 | ssa_region(jg,jreg) = (scat_od+scat_od_cloud) & |
|---|
| 651 | & / od_region(jg,jreg) |
|---|
| 652 | g_region(jreg) = (scat_od*g(jg,jlev,jcol) & |
|---|
| 653 | & + scat_od_cloud * g_cloud(iband,jlev,jcol)) & |
|---|
| 654 | & / (scat_od + scat_od_cloud) |
|---|
| 655 | |
|---|
| 656 | ! Apply maximum cloud optical depth for stability in the |
|---|
| 657 | ! 3D case |
|---|
| 658 | if (od_region(jg,jreg) > config%max_cloud_od) then |
|---|
| 659 | od_region(jg,jreg) = config%max_cloud_od |
|---|
| 660 | end if |
|---|
| 661 | |
|---|
| 662 | end do |
|---|
| 663 | |
|---|
| 664 | ! Calculate two-stream variables gamma1-gamma3 of all |
|---|
| 665 | ! regions at once |
|---|
| 666 | call calc_two_stream_gammas_sw(nreg, & |
|---|
| 667 | & mu0, ssa_region(jg,:), g_region, & |
|---|
| 668 | & gamma1(jg,:), gamma2(jg,:), gamma3(jg,:)) |
|---|
| 669 | |
|---|
| 670 | ! Loop is in order of g-points with typically |
|---|
| 671 | ! increasing optical depth: if optical depth of |
|---|
| 672 | ! clear-sky region exceeds a threshold then turn off |
|---|
| 673 | ! 3D effects for any further g-points |
|---|
| 674 | if (ng3D == ng & |
|---|
| 675 | & .and. od_region(jg,1) > config%max_gas_od_3D) then |
|---|
| 676 | ng3D = jg-1 |
|---|
| 677 | end if |
|---|
| 678 | end do ! Loop over g points |
|---|
| 679 | end if ! Cloudy level |
|---|
| 680 | |
|---|
| 681 | ! -------------------------------------------------------- |
|---|
| 682 | ! Section 3.3: Compute reflection, transmission and emission |
|---|
| 683 | ! -------------------------------------------------------- |
|---|
| 684 | if (ng3D > 0) then |
|---|
| 685 | ! --- Section 3.3a: g-points with 3D effects ---------- |
|---|
| 686 | |
|---|
| 687 | ! 3D effects need to be represented in "ng3D" of the g |
|---|
| 688 | ! points. This is done by creating ng3D square matrices |
|---|
| 689 | ! each of dimension 3*nreg by 3*nreg, computing the matrix |
|---|
| 690 | ! exponential, then computing the various |
|---|
| 691 | ! transmission/reflectance matrices from that. |
|---|
| 692 | do jreg = 1,nregactive |
|---|
| 693 | ! Write the diagonal elements of -Gamma1*z1 |
|---|
| 694 | Gamma_z1(1:ng3D,jreg,jreg) & |
|---|
| 695 | & = od_region(1:ng3D,jreg)*gamma1(1:ng3D,jreg) |
|---|
| 696 | ! Write the diagonal elements of +Gamma2*z1 |
|---|
| 697 | Gamma_z1(1:ng3D,jreg+nreg,jreg) & |
|---|
| 698 | & = od_region(1:ng3D,jreg)*gamma2(1:ng3D,jreg) |
|---|
| 699 | ! Write the diagonal elements of -Gamma3*z1 |
|---|
| 700 | Gamma_z1(1:ng3D,jreg,jreg+2*nreg) & |
|---|
| 701 | & = -od_region(1:ng3D,jreg)*ssa_region(1:ng3D,jreg) & |
|---|
| 702 | & * gamma3(1:ng3D,jreg) |
|---|
| 703 | |
|---|
| 704 | ! Write the diagonal elements of +Gamma4*z1 |
|---|
| 705 | Gamma_z1(1:ng3D,jreg+nreg,jreg+2*nreg) & |
|---|
| 706 | & = od_region(1:ng3D,jreg)*ssa_region(1:ng3D,jreg) & |
|---|
| 707 | & * (1.0_jprb - gamma3(1:ng3D,jreg)) |
|---|
| 708 | |
|---|
| 709 | ! Write the diagonal elements of +Gamma0*z1 |
|---|
| 710 | Gamma_z1(1:ng3D,jreg+2*nreg,jreg+2*nreg) & |
|---|
| 711 | & = -od_region(1:ng3D,jreg)*one_over_mu0 |
|---|
| 712 | end do |
|---|
| 713 | |
|---|
| 714 | do jreg = 1,nregactive-1 |
|---|
| 715 | ! Write the elements of -Gamma1*z1 concerned with 3D |
|---|
| 716 | ! transport |
|---|
| 717 | Gamma_z1(1:ng3D,jreg,jreg) = Gamma_z1(1:ng3D,jreg,jreg) & |
|---|
| 718 | & + transfer_rate_diffuse(jreg,jreg+1) |
|---|
| 719 | Gamma_z1(1:ng3D,jreg+1,jreg+1) = Gamma_z1(1:ng3D,jreg+1,jreg+1) & |
|---|
| 720 | & + transfer_rate_diffuse(jreg+1,jreg) |
|---|
| 721 | Gamma_z1(1:ng3D,jreg+1,jreg) = -transfer_rate_diffuse(jreg,jreg+1) |
|---|
| 722 | Gamma_z1(1:ng3D,jreg,jreg+1) = -transfer_rate_diffuse(jreg+1,jreg) |
|---|
| 723 | ! Write the elements of +Gamma0*z1 concerned with 3D |
|---|
| 724 | ! transport |
|---|
| 725 | Gamma_z1(1:ng3D,jreg+2*nreg,jreg+2*nreg) & |
|---|
| 726 | & = Gamma_z1(1:ng3D,jreg+2*nreg,jreg+2*nreg) & |
|---|
| 727 | & - transfer_rate_direct(jreg,jreg+1) |
|---|
| 728 | Gamma_z1(1:ng3D,jreg+2*nreg+1,jreg+2*nreg+1) & |
|---|
| 729 | & = Gamma_z1(1:ng3D,jreg+2*nreg+1,jreg+2*nreg+1) & |
|---|
| 730 | & - transfer_rate_direct(jreg+1,jreg) |
|---|
| 731 | Gamma_z1(1:ng3D,jreg+2*nreg+1,jreg+2*nreg) & |
|---|
| 732 | & = transfer_rate_direct(jreg,jreg+1) |
|---|
| 733 | Gamma_z1(1:ng3D,jreg+2*nreg,jreg+2*nreg+1) & |
|---|
| 734 | & = transfer_rate_direct(jreg+1,jreg) |
|---|
| 735 | end do |
|---|
| 736 | |
|---|
| 737 | ! Possible flow between regions a and c |
|---|
| 738 | if (edge_length(3,jlev) > 0.0_jprb) then |
|---|
| 739 | ! Diffuse transport |
|---|
| 740 | Gamma_z1(1:ng3D,1,1) = Gamma_z1(1:ng3D,1,1) & |
|---|
| 741 | & + transfer_rate_diffuse(1,3) |
|---|
| 742 | Gamma_z1(1:ng3D,3,3) = Gamma_z1(1:ng3D,3,3) & |
|---|
| 743 | & + transfer_rate_diffuse(3,1) |
|---|
| 744 | Gamma_z1(1:ng3D,3,1) = -transfer_rate_diffuse(1,3) |
|---|
| 745 | Gamma_z1(1:ng3D,1,3) = -transfer_rate_diffuse(3,1) |
|---|
| 746 | ! Direct transport |
|---|
| 747 | Gamma_z1(1:ng3D,1+2*nreg,1+2*nreg) = Gamma_z1(1:ng3D,1+2*nreg,1+2*nreg) & |
|---|
| 748 | & - transfer_rate_direct(1,3) |
|---|
| 749 | Gamma_z1(1:ng3D,3+2*nreg,3+2*nreg) = Gamma_z1(1:ng3D,3+2*nreg,3+2*nreg) & |
|---|
| 750 | & - transfer_rate_direct(3,1) |
|---|
| 751 | Gamma_z1(1:ng3D,3+2*nreg,1+2*nreg) = transfer_rate_direct(1,3) |
|---|
| 752 | Gamma_z1(1:ng3D,1+2*nreg,3+2*nreg) = transfer_rate_direct(3,1) |
|---|
| 753 | end if |
|---|
| 754 | |
|---|
| 755 | ! Copy Gamma1*z1 |
|---|
| 756 | Gamma_z1(1:ng3D,nreg+1:nreg+nregactive,nreg+1:nreg+nregactive) & |
|---|
| 757 | & = -Gamma_z1(1:ng3D,1:nregactive,1:nregactive) |
|---|
| 758 | ! Copy Gamma2*z1 |
|---|
| 759 | Gamma_z1(1:ng3D,1:nregactive,nreg+1:nreg+nregactive) & |
|---|
| 760 | & = -Gamma_z1(1:ng3D,nreg+1:nreg+nregactive,1:nregactive) |
|---|
| 761 | |
|---|
| 762 | ! Compute the matrix exponential of Gamma_z1, returning the |
|---|
| 763 | ! result in-place |
|---|
| 764 | call expm(ng, ng3D, 3*nreg, Gamma_z1, IMatrixPatternShortwave) |
|---|
| 765 | |
|---|
| 766 | ! Update count of expm calls |
|---|
| 767 | n_calls_expm = n_calls_expm + ng3D |
|---|
| 768 | |
|---|
| 769 | ! Direct transmission matrix |
|---|
| 770 | trans_dir_dir(1:ng3D,:,:,jlev) = min(1.0_jprb,max(0.0_jprb, & |
|---|
| 771 | & Gamma_z1(1:ng3D,2*nreg+1:3*nreg, 2*nreg+1:3*nreg))) |
|---|
| 772 | ! Diffuse reflectance matrix; security on negative values |
|---|
| 773 | ! necessary occasionally for very low cloud fraction and very high |
|---|
| 774 | ! in-cloud optical depth |
|---|
| 775 | reflectance(1:ng3D,:,:,jlev) = min(1.0_jprb,max(0.0_jprb, & |
|---|
| 776 | & -solve_mat(ng,ng3D,nreg,Gamma_z1(1:ng3D,1:nreg,1:nreg), & |
|---|
| 777 | & Gamma_z1(1:ng3D,1:nreg,nreg+1:2*nreg)))) |
|---|
| 778 | ! Diffuse transmission matrix |
|---|
| 779 | transmittance(1:ng3D,:,:,jlev) = min(1.0_jprb,max(0.0_jprb, & |
|---|
| 780 | & mat_x_mat(ng,ng3D,nreg,Gamma_z1(1:ng3D,nreg+1:2*nreg,1:nreg), & |
|---|
| 781 | & reflectance(1:ng3D,:,:,jlev)) & |
|---|
| 782 | & + Gamma_z1(1:ng3D,nreg+1:2*nreg,nreg+1:2*nreg))) |
|---|
| 783 | ! Transfer matrix between downward direct and upward |
|---|
| 784 | ! diffuse |
|---|
| 785 | ref_dir(1:ng3D,:,:,jlev) = min(mu0,max(0.0_jprb, & |
|---|
| 786 | & -solve_mat(ng,ng3D,nreg,Gamma_z1(1:ng3D,1:nreg,1:nreg), & |
|---|
| 787 | & Gamma_z1(1:ng3D,1:nreg,2*nreg+1:3*nreg)))) |
|---|
| 788 | ! Transfer matrix between downward direct and downward |
|---|
| 789 | ! diffuse in layer interface below. Include correction for |
|---|
| 790 | ! trans_dir_diff out of plausible bounds (note that Meador & |
|---|
| 791 | ! Weaver has the same correction in radiation_two_stream.F90 |
|---|
| 792 | ! - this is not just an expm thing) |
|---|
| 793 | trans_dir_diff(1:ng3D,:,:,jlev) = min(mu0,max(0.0_jprb, & |
|---|
| 794 | & mat_x_mat(ng,ng3D,nreg,Gamma_z1(1:ng3D,nreg+1:2*nreg,1:nreg), & |
|---|
| 795 | & ref_dir(1:ng3D,:,:,jlev)) & |
|---|
| 796 | & + Gamma_z1(1:ng3D,nreg+1:2*nreg,2*nreg+1:3*nreg))) |
|---|
| 797 | |
|---|
| 798 | end if ! we are treating 3D effects for some g points |
|---|
| 799 | |
|---|
| 800 | ! --- Section 3.3b: g-points without 3D effects ---------- |
|---|
| 801 | |
|---|
| 802 | ! Compute reflectance, transmittance and associated terms for |
|---|
| 803 | ! clear skies, using the Meador-Weaver formulas |
|---|
| 804 | call calc_reflectance_transmittance_sw(ng, & |
|---|
| 805 | & mu0, od_region(1:ng,1), ssa_region(1:ng,1), & |
|---|
| 806 | & gamma1(1:ng,1), gamma2(1:ng,1), gamma3(1:ng,1), & |
|---|
| 807 | & ref_clear(1:ng,jlev), trans_clear(1:ng,jlev), & |
|---|
| 808 | & ref_dir_clear(1:ng,jlev), trans_dir_diff_clear(1:ng,jlev), & |
|---|
| 809 | & trans_dir_dir_clear(1:ng,jlev) ) |
|---|
| 810 | |
|---|
| 811 | n_calls_meador_weaver = n_calls_meador_weaver + ng |
|---|
| 812 | |
|---|
| 813 | if (ng3D < ng) then |
|---|
| 814 | ! Some of the g points are to be treated using the |
|---|
| 815 | ! conventional plane-parallel method. First zero the |
|---|
| 816 | ! relevant parts of the matrices |
|---|
| 817 | trans_dir_dir (ng3D+1:ng,:,:,jlev) = 0.0_jprb |
|---|
| 818 | reflectance (ng3D+1:ng,:,:,jlev) = 0.0_jprb |
|---|
| 819 | transmittance (ng3D+1:ng,:,:,jlev) = 0.0_jprb |
|---|
| 820 | ref_dir (ng3D+1:ng,:,:,jlev) = 0.0_jprb |
|---|
| 821 | trans_dir_diff(ng3D+1:ng,:,:,jlev) = 0.0_jprb |
|---|
| 822 | |
|---|
| 823 | ! Since there is no lateral transport, the clear-sky parts |
|---|
| 824 | ! of the arrays can be copied from the clear-sky arrays |
|---|
| 825 | trans_dir_dir (ng3D+1:ng,1,1,jlev) = trans_dir_dir_clear (ng3D+1:ng,jlev) |
|---|
| 826 | reflectance (ng3D+1:ng,1,1,jlev) = ref_clear (ng3D+1:ng,jlev) |
|---|
| 827 | transmittance (ng3D+1:ng,1,1,jlev) = trans_clear (ng3D+1:ng,jlev) |
|---|
| 828 | ref_dir (ng3D+1:ng,1,1,jlev) = ref_dir_clear (ng3D+1:ng,jlev) |
|---|
| 829 | trans_dir_diff(ng3D+1:ng,1,1,jlev) = trans_dir_diff_clear(ng3D+1:ng,jlev) |
|---|
| 830 | |
|---|
| 831 | ! Compute reflectance, transmittance and associated terms |
|---|
| 832 | ! for each cloudy region, using the Meador-Weaver formulas |
|---|
| 833 | do jreg = 2, nregactive |
|---|
| 834 | call calc_reflectance_transmittance_sw(ng-ng3D, & |
|---|
| 835 | & mu0, & |
|---|
| 836 | & od_region(ng3D+1:ng,jreg), ssa_region(ng3D+1:ng,jreg), & |
|---|
| 837 | & gamma1(ng3D+1:ng,jreg), gamma2(ng3D+1:ng,jreg), & |
|---|
| 838 | & gamma3(ng3D+1:ng,jreg), & |
|---|
| 839 | & reflectance(ng3D+1:ng,jreg,jreg,jlev), & |
|---|
| 840 | & transmittance(ng3D+1:ng,jreg,jreg,jlev), & |
|---|
| 841 | & ref_dir(ng3D+1:ng,jreg,jreg,jlev), & |
|---|
| 842 | & trans_dir_diff(ng3D+1:ng,jreg,jreg,jlev), & |
|---|
| 843 | & trans_dir_dir(ng3D+1:ng,jreg,jreg,jlev) ) |
|---|
| 844 | end do |
|---|
| 845 | n_calls_meador_weaver & |
|---|
| 846 | & = n_calls_meador_weaver + (ng-ng3D)*(nregactive-1) |
|---|
| 847 | end if |
|---|
| 848 | |
|---|
| 849 | end do ! Loop over levels |
|---|
| 850 | |
|---|
| 851 | ! -------------------------------------------------------- |
|---|
| 852 | ! Section 4: Compute total albedos |
|---|
| 853 | ! -------------------------------------------------------- |
|---|
| 854 | |
|---|
| 855 | total_albedo(:,:,:,:) = 0.0_jprb |
|---|
| 856 | total_albedo_direct(:,:,:,:) = 0.0_jprb |
|---|
| 857 | |
|---|
| 858 | if (config%do_clear) then |
|---|
| 859 | total_albedo_clear(:,:) = 0.0_jprb |
|---|
| 860 | total_albedo_clear_direct(:,:) = 0.0_jprb |
|---|
| 861 | end if |
|---|
| 862 | |
|---|
| 863 | ! Calculate the upwelling radiation scattered from the direct |
|---|
| 864 | ! beam incident on the surface, and copy the surface albedo |
|---|
| 865 | ! into total_albedo |
|---|
| 866 | do jreg = 1,nreg |
|---|
| 867 | do jg = 1,ng |
|---|
| 868 | total_albedo(jg,jreg,jreg,nlev+1) = albedo_diffuse(jg,jcol) |
|---|
| 869 | total_albedo_direct(jg,jreg,jreg,nlev+1) & |
|---|
| 870 | & = mu0 * albedo_direct(jg,jcol) |
|---|
| 871 | end do |
|---|
| 872 | end do |
|---|
| 873 | |
|---|
| 874 | if (config%do_clear) then |
|---|
| 875 | ! Surface albedo is the same |
|---|
| 876 | total_albedo_clear(1:ng,nlev+1) = total_albedo(1:ng,1,1,nlev+1) |
|---|
| 877 | total_albedo_clear_direct(1:ng,nlev+1) & |
|---|
| 878 | & = total_albedo_direct(1:ng,1,1,nlev+1) |
|---|
| 879 | end if |
|---|
| 880 | |
|---|
| 881 | ! Horizontal migration distances of reflected radiation at the |
|---|
| 882 | ! surface are zero |
|---|
| 883 | x_diffuse = 0.0_jprb |
|---|
| 884 | x_direct = 0.0_jprb |
|---|
| 885 | |
|---|
| 886 | ! Work back up through the atmosphere computing the total albedo |
|---|
| 887 | ! of the atmosphere below that point using the adding method |
|---|
| 888 | do jlev = nlev,1,-1 |
|---|
| 889 | |
|---|
| 890 | ! -------------------------------------------------------- |
|---|
| 891 | ! Section 4.1: Adding method |
|---|
| 892 | ! -------------------------------------------------------- |
|---|
| 893 | |
|---|
| 894 | if (config%do_clear) then |
|---|
| 895 | ! Use adding method for clear-sky arrays; note that there |
|---|
| 896 | ! is no need to consider "above" and "below" quantities |
|---|
| 897 | ! since with no cloud overlap to worry about, these are |
|---|
| 898 | ! the same |
|---|
| 899 | inv_denom_scalar(:) = 1.0_jprb & |
|---|
| 900 | & / (1.0_jprb - total_albedo_clear(:,jlev+1)*ref_clear(:,jlev)) |
|---|
| 901 | total_albedo_clear(:,jlev) = ref_clear(:,jlev) & |
|---|
| 902 | & + trans_clear(:,jlev)*trans_clear(:,jlev)*total_albedo_clear(:,jlev+1) & |
|---|
| 903 | & * inv_denom_scalar(:) |
|---|
| 904 | total_albedo_clear_direct(:,jlev) = ref_dir_clear(:,jlev) & |
|---|
| 905 | & + (trans_dir_dir_clear(:,jlev) * total_albedo_clear_direct(:,jlev+1) & |
|---|
| 906 | & +trans_dir_diff_clear(:,jlev) * total_albedo_clear(:,jlev+1)) & |
|---|
| 907 | & * trans_clear(:,jlev) * inv_denom_scalar(:) |
|---|
| 908 | end if |
|---|
| 909 | |
|---|
| 910 | if (is_clear_sky_layer(jlev)) then |
|---|
| 911 | ! Clear-sky layer: use scalar adding method |
|---|
| 912 | inv_denom_scalar(:) = 1.0_jprb & |
|---|
| 913 | & / (1.0_jprb - total_albedo(:,1,1,jlev+1)*reflectance(:,1,1,jlev)) |
|---|
| 914 | total_albedo_below = 0.0_jprb |
|---|
| 915 | total_albedo_below(:,1,1) = reflectance(:,1,1,jlev) & |
|---|
| 916 | & + transmittance(:,1,1,jlev) * transmittance(:,1,1,jlev) & |
|---|
| 917 | & * total_albedo(:,1,1,jlev+1) * inv_denom_scalar(:) |
|---|
| 918 | total_albedo_below_direct = 0.0_jprb |
|---|
| 919 | total_albedo_below_direct(:,1,1) = ref_dir(:,1,1,jlev) & |
|---|
| 920 | & + (trans_dir_dir(:,1,1,jlev)*total_albedo_direct(:,1,1,jlev+1) & |
|---|
| 921 | & +trans_dir_diff(:,1,1,jlev)*total_albedo(:,1,1,jlev+1)) & |
|---|
| 922 | & * transmittance(:,1,1,jlev) * inv_denom_scalar(:) |
|---|
| 923 | else |
|---|
| 924 | ! Cloudy layer: use matrix adding method |
|---|
| 925 | denominator = identity_minus_mat_x_mat(ng,ng,nreg, & |
|---|
| 926 | & total_albedo(:,:,:,jlev+1), reflectance(:,:,:,jlev)) |
|---|
| 927 | total_albedo_below = reflectance(:,:,:,jlev) & |
|---|
| 928 | & + mat_x_mat(ng,ng,nreg,transmittance(:,:,:,jlev), & |
|---|
| 929 | & solve_mat(ng,ng,nreg,denominator, & |
|---|
| 930 | & mat_x_mat(ng,ng,nreg,total_albedo(:,:,:,jlev+1), & |
|---|
| 931 | & transmittance(:,:,:,jlev)))) |
|---|
| 932 | total_albedo_below_direct = ref_dir(:,:,:,jlev) & |
|---|
| 933 | & + mat_x_mat(ng,ng,nreg,transmittance(:,:,:,jlev), & |
|---|
| 934 | & solve_mat(ng,ng,nreg,denominator, & |
|---|
| 935 | & mat_x_mat(ng,ng,nreg,total_albedo_direct(:,:,:,jlev+1), & |
|---|
| 936 | & trans_dir_dir(:,:,:,jlev)) & |
|---|
| 937 | & +mat_x_mat(ng,ng,nreg,total_albedo(:,:,:,jlev+1), & |
|---|
| 938 | & trans_dir_diff(:,:,:,jlev)))) |
|---|
| 939 | end if |
|---|
| 940 | |
|---|
| 941 | ! -------------------------------------------------------- |
|---|
| 942 | ! Section 4.2: Overlap and entrapment |
|---|
| 943 | ! -------------------------------------------------------- |
|---|
| 944 | |
|---|
| 945 | #ifndef PRINT_ENTRAPMENT_DATA |
|---|
| 946 | if ((config%i_3d_sw_entrapment == IEntrapmentExplicitNonFractal & |
|---|
| 947 | & .or. config%i_3d_sw_entrapment == IEntrapmentExplicit) & |
|---|
| 948 | & .and. jlev >= i_cloud_top) then |
|---|
| 949 | #else |
|---|
| 950 | if (config%i_3d_sw_entrapment == IEntrapmentExplicitNonFractal & |
|---|
| 951 | & .or. config%i_3d_sw_entrapment == IEntrapmentExplicit) then |
|---|
| 952 | #endif |
|---|
| 953 | ! "Explicit entrapment": we have the horizontal migration |
|---|
| 954 | ! distances just above the base of the layer, and need to |
|---|
| 955 | ! step them to just below the top of the same layer |
|---|
| 956 | call step_migrations(ng, nreg, cloud%fraction(jcol,jlev), & |
|---|
| 957 | & layer_depth(jlev), tan_diffuse_angle_3d, tan_sza, & |
|---|
| 958 | & reflectance(:,:,:,jlev), transmittance(:,:,:,jlev), & |
|---|
| 959 | & ref_dir(:,:,:,jlev), trans_dir_dir(:,:,:,jlev), & |
|---|
| 960 | & trans_dir_diff(:,:,:,jlev), total_albedo(:,:,:,jlev+1), & |
|---|
| 961 | & total_albedo_direct(:,:,:,jlev+1), & |
|---|
| 962 | & x_diffuse, x_direct) |
|---|
| 963 | |
|---|
| 964 | #ifdef PRINT_ENTRAPMENT_DATA |
|---|
| 965 | ! Write out for later analysis: these are the entrapment |
|---|
| 966 | ! statistics at the top of layer "jlev" |
|---|
| 967 | ! Note that number of scattering events is now not computed, |
|---|
| 968 | ! so print "1.0" |
|---|
| 969 | if (nreg == 2) then |
|---|
| 970 | write(101,'(i4,i4,6e14.6)') jcol, jlev, & |
|---|
| 971 | & x_direct(1,:), x_diffuse(1,:), x_direct(1,:)*0.0_jprb+1.0_jprb |
|---|
| 972 | else |
|---|
| 973 | write(101,'(i4,i4,9e14.6)') jcol, jlev, & |
|---|
| 974 | & x_direct(1,1:3), x_diffuse(1,1:3), 1.0_jprb,1.0_jprb,1.0_jprb |
|---|
| 975 | end if |
|---|
| 976 | #endif |
|---|
| 977 | |
|---|
| 978 | end if |
|---|
| 979 | |
|---|
| 980 | ! Account for cloud overlap when converting albedo and source |
|---|
| 981 | ! below a layer interface to the equivalent values just above |
|---|
| 982 | if (is_clear_sky_layer(jlev) .and. is_clear_sky_layer(jlev-1)) then |
|---|
| 983 | ! If both layers are cloud free, this is trivial... |
|---|
| 984 | total_albedo(:,:,:,jlev) = 0.0_jprb |
|---|
| 985 | total_albedo(:,1,1,jlev) = total_albedo_below(:,1,1) |
|---|
| 986 | total_albedo_direct(:,:,:,jlev) = 0.0_jprb |
|---|
| 987 | total_albedo_direct(:,1,1,jlev) = total_albedo_below_direct(:,1,1) |
|---|
| 988 | |
|---|
| 989 | else if (config%i_3d_sw_entrapment == IEntrapmentMaximum & |
|---|
| 990 | & .or. is_clear_sky_layer(jlev-1)) then |
|---|
| 991 | ! "Maximum entrapment": use the overlap matrices u_matrix and v_matrix |
|---|
| 992 | ! (this is the original SPARTACUS method) |
|---|
| 993 | total_albedo(:,:,:,jlev) = singlemat_x_mat(ng,ng,nreg,& |
|---|
| 994 | & u_matrix(:,:,jlev,jcol), & |
|---|
| 995 | & mat_x_singlemat(ng,ng,nreg,total_albedo_below,& |
|---|
| 996 | & v_matrix(:,:,jlev,jcol))) |
|---|
| 997 | total_albedo_direct(:,:,:,jlev) = singlemat_x_mat(ng,ng,nreg,& |
|---|
| 998 | & u_matrix(:,:,jlev,jcol), & |
|---|
| 999 | & mat_x_singlemat(ng,ng,nreg,total_albedo_below_direct,& |
|---|
| 1000 | & v_matrix(:,:,jlev,jcol))) |
|---|
| 1001 | |
|---|
| 1002 | else if (config%i_3d_sw_entrapment == IEntrapmentZero) then |
|---|
| 1003 | ! "Zero entrapment": even radiation transported |
|---|
| 1004 | ! laterally between regions in the layers below is |
|---|
| 1005 | ! reflected back up into the same region. First diffuse |
|---|
| 1006 | ! radiation: |
|---|
| 1007 | total_albedo(:,:,:,jlev) = 0.0_jprb |
|---|
| 1008 | do jreg = 1,nreg ! Target layer (jlev-1) |
|---|
| 1009 | do jreg2 = 1,nreg ! Current layer (jlev) |
|---|
| 1010 | total_albedo(:,jreg,jreg,jlev) = total_albedo(:,jreg,jreg,jlev) & |
|---|
| 1011 | & + sum(total_albedo_below(:,:,jreg2),2) & |
|---|
| 1012 | & * v_matrix(jreg2,jreg,jlev,jcol) |
|---|
| 1013 | end do |
|---|
| 1014 | end do |
|---|
| 1015 | ! ...then direct radiation: |
|---|
| 1016 | total_albedo_direct(:,:,:,jlev) = 0.0_jprb |
|---|
| 1017 | do jreg = 1,nreg ! Target layer (jlev-1) |
|---|
| 1018 | do jreg2 = 1,nreg ! Current layer (jlev) |
|---|
| 1019 | total_albedo_direct(:,jreg,jreg,jlev) = total_albedo_direct(:,jreg,jreg,jlev) & |
|---|
| 1020 | & + sum(total_albedo_below_direct(:,:,jreg2),2) & |
|---|
| 1021 | & * v_matrix(jreg2,jreg,jlev,jcol) |
|---|
| 1022 | end do |
|---|
| 1023 | end do |
|---|
| 1024 | |
|---|
| 1025 | else |
|---|
| 1026 | ! Controlled entrapment |
|---|
| 1027 | |
|---|
| 1028 | #ifdef EXPLICIT_EDGE_ENTRAPMENT |
|---|
| 1029 | ! If "EXPLICIT_EDGE_ENTRAPMENT" is defined then we use the |
|---|
| 1030 | ! explicit entrapment approach for both horizontal transport |
|---|
| 1031 | ! within regions, and horizontal transport between regions |
|---|
| 1032 | ! (otherwise, horizontal transport between regions is |
|---|
| 1033 | ! automatically treated using maximum entrapment). This is |
|---|
| 1034 | ! experimental, which is why it is not a run-time option. |
|---|
| 1035 | |
|---|
| 1036 | if (config%i_3d_sw_entrapment == IEntrapmentEdgeOnly) then |
|---|
| 1037 | #endif |
|---|
| 1038 | ! Add the contribution from off-diagonal elements of the |
|---|
| 1039 | ! albedo matrix in the lower layer, i.e. radiation that |
|---|
| 1040 | ! flows between regions... |
|---|
| 1041 | |
|---|
| 1042 | ! First diffuse radiation: |
|---|
| 1043 | albedo_part = total_albedo_below |
|---|
| 1044 | do jreg = 1,nreg |
|---|
| 1045 | albedo_part(:,jreg,jreg) = 0.0_jprb |
|---|
| 1046 | end do |
|---|
| 1047 | total_albedo(:,:,:,jlev) = singlemat_x_mat(ng,ng,nreg,& |
|---|
| 1048 | & u_matrix(:,:,jlev,jcol), & |
|---|
| 1049 | & mat_x_singlemat(ng,ng,nreg,albedo_part,& |
|---|
| 1050 | & v_matrix(:,:,jlev,jcol))) |
|---|
| 1051 | ! ...then direct radiation: |
|---|
| 1052 | albedo_part = total_albedo_below_direct |
|---|
| 1053 | do jreg = 1,nreg |
|---|
| 1054 | albedo_part(:,jreg,jreg) = 0.0_jprb |
|---|
| 1055 | end do |
|---|
| 1056 | total_albedo_direct(:,:,:,jlev) = singlemat_x_mat(ng,ng,nreg,& |
|---|
| 1057 | & u_matrix(:,:,jlev,jcol), & |
|---|
| 1058 | & mat_x_singlemat(ng,ng,nreg,albedo_part,& |
|---|
| 1059 | & v_matrix(:,:,jlev,jcol))) |
|---|
| 1060 | |
|---|
| 1061 | #ifdef EXPLICIT_EDGE_ENTRAPMENT |
|---|
| 1062 | end if |
|---|
| 1063 | #endif |
|---|
| 1064 | |
|---|
| 1065 | ! Now the contribution from the diagonals of the albedo |
|---|
| 1066 | ! matrix in the lower layer |
|---|
| 1067 | if (config%i_3d_sw_entrapment == IEntrapmentEdgeOnly & |
|---|
| 1068 | & .or. (.not. config%do_3d_effects)) then |
|---|
| 1069 | ! "Edge-only entrapment": the operation we perform is |
|---|
| 1070 | ! essentially diag(total_albedo) += matmul(transpose(v_matrix), |
|---|
| 1071 | ! diag(total_albedo_below)). |
|---|
| 1072 | do jreg = 1,nreg |
|---|
| 1073 | do jreg2 = 1,nreg |
|---|
| 1074 | total_albedo(:,jreg,jreg,jlev) & |
|---|
| 1075 | & = total_albedo(:,jreg,jreg,jlev) & |
|---|
| 1076 | & + total_albedo_below(:,jreg2,jreg2) & |
|---|
| 1077 | & * v_matrix(jreg2,jreg,jlev,jcol) |
|---|
| 1078 | total_albedo_direct(:,jreg,jreg,jlev) & |
|---|
| 1079 | & = total_albedo_direct(:,jreg,jreg,jlev) & |
|---|
| 1080 | & + total_albedo_below_direct(:,jreg2,jreg2) & |
|---|
| 1081 | & * v_matrix(jreg2,jreg,jlev,jcol) |
|---|
| 1082 | end do |
|---|
| 1083 | end do |
|---|
| 1084 | |
|---|
| 1085 | else |
|---|
| 1086 | ! "Explicit entrapment" |
|---|
| 1087 | |
|---|
| 1088 | do jreg2 = 1,nreg |
|---|
| 1089 | ! Loop through each region in the lower layer. For one |
|---|
| 1090 | ! of the regions in the lower layer, we are imagining it |
|---|
| 1091 | ! to be divided into "nreg" subregions that map onto the |
|---|
| 1092 | ! regions in the upper layer. The rate of exchange |
|---|
| 1093 | ! between these subregions is computed via a coupled |
|---|
| 1094 | ! differential equation written in terms of a singular |
|---|
| 1095 | ! exchange matrix (there are only terms for the exchange |
|---|
| 1096 | ! between subregions, but no propagation effects since |
|---|
| 1097 | ! we already know the albedo of this region). This is |
|---|
| 1098 | ! solved using the matrix-exponential method. |
|---|
| 1099 | |
|---|
| 1100 | ! Use the following array for the transfer of either |
|---|
| 1101 | ! diffuse or direct radiation (despite the name), per |
|---|
| 1102 | ! unit horizontal distance travelled |
|---|
| 1103 | transfer_rate_diffuse = 0.0_jprb |
|---|
| 1104 | |
|---|
| 1105 | ! As we need to reference the layer above the interface, |
|---|
| 1106 | ! don't do the following on the highest layer |
|---|
| 1107 | if (jlev > 1) then |
|---|
| 1108 | |
|---|
| 1109 | ! Given a horizontal migration distance, there is |
|---|
| 1110 | ! still uncertainty about how much entrapment occurs |
|---|
| 1111 | ! associated with how one assumes cloud boundaries |
|---|
| 1112 | ! line up in adjacent layers. "overhang_factor" |
|---|
| 1113 | ! can be varied between 0.0 (the boundaries line up to |
|---|
| 1114 | ! the greatest extent possible given the overlap |
|---|
| 1115 | ! parameter) and 1.0 (the boundaries line up to the |
|---|
| 1116 | ! minimum extent possible); here this is used to |
|---|
| 1117 | ! produce a scaling factor for the transfer rate. |
|---|
| 1118 | transfer_scaling = 1.0_jprb - (1.0_jprb - config%overhang_factor) & |
|---|
| 1119 | & * cloud%overlap_param(jcol,jlev-1) & |
|---|
| 1120 | & * min(region_fracs(jreg2,jlev,jcol),region_fracs(jreg2,jlev-1,jcol)) & |
|---|
| 1121 | & / max(config%cloud_fraction_threshold, region_fracs(jreg2,jlev,jcol)) |
|---|
| 1122 | |
|---|
| 1123 | do jreg = 1, nreg-1 |
|---|
| 1124 | ! Compute lateral transfer rates from region jreg to |
|---|
| 1125 | ! jreg+1 as before, but without length scale which |
|---|
| 1126 | ! is wavelength dependent. |
|---|
| 1127 | |
|---|
| 1128 | ! Recall that overlap indexing is |
|---|
| 1129 | ! u_matrix(upper_region, lower_region, level, |
|---|
| 1130 | ! column). |
|---|
| 1131 | transfer_rate_diffuse(jreg,jreg+1) = transfer_scaling & |
|---|
| 1132 | & * edge_length(jreg,jlev-1) / max(u_matrix(jreg,jreg2,jlev,jcol),1.0e-5_jprb) |
|---|
| 1133 | ! Compute transfer rates from region jreg+1 to jreg |
|---|
| 1134 | transfer_rate_diffuse(jreg+1,jreg) = transfer_scaling & |
|---|
| 1135 | & * edge_length(jreg,jlev-1) / max(u_matrix(jreg+1,jreg2,jlev,jcol),1.0e-5_jprb) |
|---|
| 1136 | end do |
|---|
| 1137 | |
|---|
| 1138 | ! Compute transfer rates directly between regions 1 |
|---|
| 1139 | ! and 3 (not used below) |
|---|
| 1140 | if (edge_length(3,jlev) > 0.0_jprb) then |
|---|
| 1141 | transfer_rate_diffuse(1,3) = transfer_scaling & |
|---|
| 1142 | & * edge_length(3,jlev-1) / max(u_matrix(1,jreg2,jlev,jcol),1.0e-5_jprb) |
|---|
| 1143 | transfer_rate_diffuse(3,1) = transfer_scaling & |
|---|
| 1144 | & * edge_length(3,jlev-1) / max(u_matrix(3,jreg2,jlev,jcol),1.0e-5_jprb) |
|---|
| 1145 | end if |
|---|
| 1146 | end if |
|---|
| 1147 | |
|---|
| 1148 | ! Compute matrix of exchange coefficients |
|---|
| 1149 | entrapment = 0.0_jprb |
|---|
| 1150 | inv_effective_size = min(cloud%inv_cloud_effective_size(jcol,jlev-1), & |
|---|
| 1151 | & 1.0_jprb/config%min_cloud_effective_size) |
|---|
| 1152 | do jreg = 1,nreg-1 |
|---|
| 1153 | ! Diffuse transport down and up with one random |
|---|
| 1154 | ! scattering event |
|---|
| 1155 | if (config%i_3d_sw_entrapment == IEntrapmentExplicit) then |
|---|
| 1156 | fractal_factor = 1.0_jprb / sqrt(max(1.0_jprb, 2.5_jprb*x_diffuse(:,jreg2) & |
|---|
| 1157 | & * inv_effective_size)) |
|---|
| 1158 | entrapment(:,jreg+1,jreg) = entrapment(:,jreg+1,jreg) & |
|---|
| 1159 | & + transfer_rate_diffuse(jreg,jreg+1)*x_diffuse(:,jreg2) & |
|---|
| 1160 | & * fractal_factor |
|---|
| 1161 | entrapment(:,jreg,jreg+1) = entrapment(:,jreg,jreg+1) & |
|---|
| 1162 | & + transfer_rate_diffuse(jreg+1,jreg)*x_diffuse(:,jreg2) & |
|---|
| 1163 | & * fractal_factor |
|---|
| 1164 | else |
|---|
| 1165 | entrapment(:,jreg+1,jreg) = entrapment(:,jreg+1,jreg) & |
|---|
| 1166 | & + transfer_rate_diffuse(jreg,jreg+1)*x_diffuse(:,jreg2) |
|---|
| 1167 | entrapment(:,jreg,jreg+1) = entrapment(:,jreg,jreg+1) & |
|---|
| 1168 | & + transfer_rate_diffuse(jreg+1,jreg)*x_diffuse(:,jreg2) |
|---|
| 1169 | end if |
|---|
| 1170 | entrapment(:,jreg,jreg) = entrapment(:,jreg,jreg) & |
|---|
| 1171 | & - entrapment(:,jreg+1,jreg) |
|---|
| 1172 | entrapment(:,jreg+1,jreg+1) = entrapment(:,jreg+1,jreg+1) & |
|---|
| 1173 | & - entrapment(:,jreg,jreg+1) |
|---|
| 1174 | end do |
|---|
| 1175 | |
|---|
| 1176 | ! If rate of exchange is excessive the expm can throw a |
|---|
| 1177 | ! floating point exception, even if it tends towards a |
|---|
| 1178 | ! trival limit, so we cap the maximum input to expm by |
|---|
| 1179 | ! scaling down if necessary |
|---|
| 1180 | do jg = 1,ng |
|---|
| 1181 | max_entr = -min(entrapment(jg,1,1),entrapment(jg,2,2)) |
|---|
| 1182 | if (max_entr > config%max_cloud_od) then |
|---|
| 1183 | ! Scale down all inputs for this g point |
|---|
| 1184 | entrapment(jg,:,:) = entrapment(jg,:,:) * (config%max_cloud_od/max_entr) |
|---|
| 1185 | end if |
|---|
| 1186 | end do |
|---|
| 1187 | |
|---|
| 1188 | ! Since the matrix to be exponentiated has a simple |
|---|
| 1189 | ! structure we may use a faster method described in the |
|---|
| 1190 | ! appendix of Hogan et al. (GMD 2018) |
|---|
| 1191 | #define USE_FAST_EXPM_EXCHANGE 1 |
|---|
| 1192 | #ifdef USE_FAST_EXPM_EXCHANGE |
|---|
| 1193 | if (nreg == 2) then |
|---|
| 1194 | call fast_expm_exchange(ng, ng, entrapment(:,2,1), entrapment(:,1,2), & |
|---|
| 1195 | & albedo_part) |
|---|
| 1196 | else |
|---|
| 1197 | call fast_expm_exchange(ng, ng, entrapment(:,2,1), entrapment(:,1,2), & |
|---|
| 1198 | & entrapment(:,3,2), entrapment(:,2,3), & |
|---|
| 1199 | & albedo_part) |
|---|
| 1200 | end if |
|---|
| 1201 | #else |
|---|
| 1202 | ! Use matrix exponential to compute rate of exchange |
|---|
| 1203 | albedo_part = entrapment |
|---|
| 1204 | call expm(ng, ng, nreg, albedo_part, IMatrixPatternDense) |
|---|
| 1205 | n_calls_expm = n_calls_expm + ng |
|---|
| 1206 | #endif |
|---|
| 1207 | |
|---|
| 1208 | #ifndef EXPLICIT_EDGE_ENTRAPMENT |
|---|
| 1209 | ! Scale to get the contribution to the diffuse albedo |
|---|
| 1210 | do jreg3 = 1,nreg |
|---|
| 1211 | do jreg = 1,nreg |
|---|
| 1212 | albedo_part(:,jreg3,jreg) = albedo_part(:,jreg3,jreg) & |
|---|
| 1213 | & * v_matrix(jreg2,jreg,jlev,jcol) * total_albedo_below(:,jreg2,jreg2) |
|---|
| 1214 | end do |
|---|
| 1215 | end do |
|---|
| 1216 | #else |
|---|
| 1217 | ! The following is an experimental treatment that tries |
|---|
| 1218 | ! to explicitly account for the horizontal distance |
|---|
| 1219 | ! traveled by radiation that passes through cloud sides |
|---|
| 1220 | entrapment = albedo_part |
|---|
| 1221 | albedo_part = 0.0_jprb |
|---|
| 1222 | ! Scale to get the contribution to the diffuse albedo |
|---|
| 1223 | do jreg3 = 1,nreg ! TO upper region |
|---|
| 1224 | do jreg = 1,nreg ! FROM upper region |
|---|
| 1225 | transfer_scaling = 1.0_jprb - (1.0_jprb - config%overhang_factor) & |
|---|
| 1226 | & * cloud%overlap_param(jcol,jlev-1) & |
|---|
| 1227 | & * min(region_fracs(jreg,jlev,jcol), region_fracs(jreg,jlev,jcol)) & |
|---|
| 1228 | & / max(config%cloud_fraction_threshold, region_fracs(jreg,jlev,jcol)) |
|---|
| 1229 | do jreg4 = 1,nreg ! VIA first lower region (jreg2 is second lower region) |
|---|
| 1230 | if (.not. (jreg4 == jreg .and. jreg4 /= jreg2)) then |
|---|
| 1231 | albedo_part(:,jreg3,jreg) = albedo_part(:,jreg3,jreg) + entrapment(:,jreg3,jreg) & |
|---|
| 1232 | & * v_matrix(jreg4,jreg,jlev,jcol) * total_albedo_below(:,jreg2,jreg4) |
|---|
| 1233 | else |
|---|
| 1234 | albedo_part(:,jreg3,jreg) = albedo_part(:,jreg3,jreg) & |
|---|
| 1235 | & + v_matrix(jreg4,jreg,jlev,jcol) * total_albedo_below(:,jreg2,jreg4) & |
|---|
| 1236 | & * (transfer_scaling * entrapment(:,jreg3,jreg) & |
|---|
| 1237 | & +((1.0_jprb-transfer_scaling) * entrapment(:,jreg3,jreg2))) |
|---|
| 1238 | end if |
|---|
| 1239 | end do |
|---|
| 1240 | end do |
|---|
| 1241 | end do |
|---|
| 1242 | #endif |
|---|
| 1243 | |
|---|
| 1244 | ! Increment diffuse albedo |
|---|
| 1245 | total_albedo(:,:,:,jlev) = total_albedo(:,:,:,jlev) + albedo_part |
|---|
| 1246 | |
|---|
| 1247 | ! Now do the same for the direct albedo |
|---|
| 1248 | entrapment = 0.0_jprb |
|---|
| 1249 | do jreg = 1,nreg-1 |
|---|
| 1250 | ! Direct transport down and diffuse up with one random |
|---|
| 1251 | ! scattering event |
|---|
| 1252 | if (config%i_3d_sw_entrapment == IEntrapmentExplicit) then |
|---|
| 1253 | fractal_factor = 1.0_jprb / sqrt(max(1.0_jprb, 2.5_jprb*x_direct(:,jreg2) & |
|---|
| 1254 | & * inv_effective_size)) |
|---|
| 1255 | entrapment(:,jreg+1,jreg) = entrapment(:,jreg+1,jreg) & |
|---|
| 1256 | & + transfer_rate_diffuse(jreg,jreg+1)*x_direct(:,jreg2) & |
|---|
| 1257 | & * fractal_factor |
|---|
| 1258 | entrapment(:,jreg,jreg+1) = entrapment(:,jreg,jreg+1) & |
|---|
| 1259 | & + transfer_rate_diffuse(jreg+1,jreg)*x_direct(:,jreg2) & |
|---|
| 1260 | & * fractal_factor |
|---|
| 1261 | else |
|---|
| 1262 | entrapment(:,jreg+1,jreg) = entrapment(:,jreg+1,jreg) & |
|---|
| 1263 | & + transfer_rate_diffuse(jreg,jreg+1)*x_direct(:,jreg2) |
|---|
| 1264 | entrapment(:,jreg,jreg+1) = entrapment(:,jreg,jreg+1) & |
|---|
| 1265 | & + transfer_rate_diffuse(jreg+1,jreg)*x_direct(:,jreg2) |
|---|
| 1266 | end if |
|---|
| 1267 | entrapment(:,jreg,jreg) = entrapment(:,jreg,jreg) & |
|---|
| 1268 | & - entrapment(:,jreg+1,jreg) |
|---|
| 1269 | entrapment(:,jreg+1,jreg+1) = entrapment(:,jreg+1,jreg+1) & |
|---|
| 1270 | & - entrapment(:,jreg,jreg+1) |
|---|
| 1271 | end do |
|---|
| 1272 | |
|---|
| 1273 | ! If rate of exchange is excessive the expm can throw a |
|---|
| 1274 | ! floating point exception, even if it tends towards a |
|---|
| 1275 | ! trival limit, so we cap the maximum input to expm by |
|---|
| 1276 | ! scaling down if necessary |
|---|
| 1277 | do jg = 1,ng |
|---|
| 1278 | max_entr = -min(entrapment(jg,1,1),entrapment(jg,2,2)) |
|---|
| 1279 | if (max_entr > config%max_cloud_od) then |
|---|
| 1280 | ! Scale down all inputs for this g point |
|---|
| 1281 | entrapment(jg,:,:) = entrapment(jg,:,:) * (config%max_cloud_od/max_entr) |
|---|
| 1282 | end if |
|---|
| 1283 | end do |
|---|
| 1284 | |
|---|
| 1285 | |
|---|
| 1286 | #ifdef USE_FAST_EXPM_EXCHANGE |
|---|
| 1287 | if (nreg == 2) then |
|---|
| 1288 | call fast_expm_exchange(ng, ng, entrapment(:,2,1), entrapment(:,1,2), & |
|---|
| 1289 | & albedo_part) |
|---|
| 1290 | else |
|---|
| 1291 | call fast_expm_exchange(ng, ng, entrapment(:,2,1), entrapment(:,1,2), & |
|---|
| 1292 | & entrapment(:,3,2), entrapment(:,2,3), & |
|---|
| 1293 | & albedo_part) |
|---|
| 1294 | end if |
|---|
| 1295 | #else |
|---|
| 1296 | albedo_part = entrapment |
|---|
| 1297 | call expm(ng, ng, nreg, albedo_part, IMatrixPatternDense) |
|---|
| 1298 | n_calls_expm = n_calls_expm + ng |
|---|
| 1299 | #endif |
|---|
| 1300 | |
|---|
| 1301 | #ifndef EXPLICIT_EDGE_ENTRAPMENT |
|---|
| 1302 | do jreg3 = 1,nreg |
|---|
| 1303 | do jreg = 1,nreg |
|---|
| 1304 | albedo_part(:,jreg3,jreg) = albedo_part(:,jreg3,jreg) & |
|---|
| 1305 | & * v_matrix(jreg2,jreg,jlev,jcol) * total_albedo_below_direct(:,jreg2,jreg2) |
|---|
| 1306 | end do |
|---|
| 1307 | end do |
|---|
| 1308 | #else |
|---|
| 1309 | entrapment = albedo_part |
|---|
| 1310 | albedo_part = 0.0_jprb |
|---|
| 1311 | do jreg3 = 1,nreg |
|---|
| 1312 | do jreg = 1,nreg |
|---|
| 1313 | transfer_scaling = 1.0_jprb - (1.0_jprb - config%overhang_factor) & |
|---|
| 1314 | & * cloud%overlap_param(jcol,jlev-1) & |
|---|
| 1315 | & * min(region_fracs(jreg,jlev,jcol), region_fracs(jreg,jlev-1,jcol)) & |
|---|
| 1316 | & / max(config%cloud_fraction_threshold, region_fracs(jreg,jlev,jcol)) |
|---|
| 1317 | do jreg4 = 1,nreg |
|---|
| 1318 | if (.not. (jreg4 == jreg .and. jreg4 /= jreg2)) then |
|---|
| 1319 | albedo_part(:,jreg3,jreg) = albedo_part(:,jreg3,jreg) + entrapment(:,jreg3,jreg) & |
|---|
| 1320 | & * v_matrix(jreg4,jreg,jlev,jcol) * total_albedo_below_direct(:,jreg2,jreg4) |
|---|
| 1321 | else |
|---|
| 1322 | albedo_part(:,jreg3,jreg) = albedo_part(:,jreg3,jreg) & |
|---|
| 1323 | & + v_matrix(jreg4,jreg,jlev,jcol) * total_albedo_below_direct(:,jreg2,jreg4) & |
|---|
| 1324 | & * (transfer_scaling * entrapment(:,jreg3,jreg) & |
|---|
| 1325 | & +((1.0_jprb-transfer_scaling) * entrapment(:,jreg3,jreg2))) |
|---|
| 1326 | end if |
|---|
| 1327 | end do |
|---|
| 1328 | end do |
|---|
| 1329 | end do |
|---|
| 1330 | |
|---|
| 1331 | #endif |
|---|
| 1332 | ! Increment direct albedo |
|---|
| 1333 | total_albedo_direct(:,:,:,jlev) = total_albedo_direct(:,:,:,jlev) + albedo_part |
|---|
| 1334 | |
|---|
| 1335 | end do |
|---|
| 1336 | |
|---|
| 1337 | end if |
|---|
| 1338 | end if |
|---|
| 1339 | |
|---|
| 1340 | if ((config%i_3d_sw_entrapment == IEntrapmentExplicitNonFractal & |
|---|
| 1341 | & .or. config%i_3d_sw_entrapment == IEntrapmentExplicit) & |
|---|
| 1342 | & .and. .not. (is_clear_sky_layer(jlev) .and. is_clear_sky_layer(jlev-1))) then |
|---|
| 1343 | ! Horizontal migration distances are averaged when |
|---|
| 1344 | ! applying overlap rules, so equation is |
|---|
| 1345 | ! x_above=matmul(transpose(v_matrix),x_below) |
|---|
| 1346 | |
|---|
| 1347 | ! We do this into temporary arrays... |
|---|
| 1348 | x_direct_above = 0.0_jprb |
|---|
| 1349 | x_diffuse_above = 0.0_jprb |
|---|
| 1350 | |
|---|
| 1351 | nregactive = nreg |
|---|
| 1352 | if (is_clear_sky_layer(jlev)) then |
|---|
| 1353 | nregactive = 1 |
|---|
| 1354 | end if |
|---|
| 1355 | |
|---|
| 1356 | do jreg = 1,nreg ! Target layer (jlev-1) |
|---|
| 1357 | do jreg2 = 1,nregactive ! Current layer (jlev) |
|---|
| 1358 | x_direct_above(:,jreg) = x_direct_above(:,jreg) & |
|---|
| 1359 | & + x_direct(:,jreg2) * v_matrix(jreg2,jreg,jlev,jcol) |
|---|
| 1360 | x_diffuse_above(:,jreg) = x_diffuse_above(:,jreg) & |
|---|
| 1361 | & + x_diffuse(:,jreg2) * v_matrix(jreg2,jreg,jlev,jcol) |
|---|
| 1362 | end do |
|---|
| 1363 | end do |
|---|
| 1364 | |
|---|
| 1365 | !... then copy out of the temporary arrays |
|---|
| 1366 | x_direct = x_direct_above |
|---|
| 1367 | x_diffuse = x_diffuse_above |
|---|
| 1368 | end if |
|---|
| 1369 | |
|---|
| 1370 | end do ! Reverse loop over levels |
|---|
| 1371 | |
|---|
| 1372 | ! -------------------------------------------------------- |
|---|
| 1373 | ! Section 5: Compute fluxes |
|---|
| 1374 | ! -------------------------------------------------------- |
|---|
| 1375 | |
|---|
| 1376 | ! Top-of-atmosphere fluxes into the regions of the top-most |
|---|
| 1377 | ! layer, zero since we assume no diffuse downwelling |
|---|
| 1378 | flux_dn_below = 0.0_jprb |
|---|
| 1379 | ! Direct downwelling flux (into a plane perpendicular to the |
|---|
| 1380 | ! sun) entering the top of each region in the top-most layer |
|---|
| 1381 | do jreg = 1,nreg |
|---|
| 1382 | direct_dn_below(:,jreg) = incoming_sw(:,jcol)*region_fracs(jreg,1,jcol) |
|---|
| 1383 | end do |
|---|
| 1384 | ! We're using flux_up_above as a container; actually its |
|---|
| 1385 | ! interpretation at top of atmosphere here is just 'below' the |
|---|
| 1386 | ! TOA interface, so using the regions of the first model layer |
|---|
| 1387 | flux_up_above = mat_x_vec(ng,ng,nreg,total_albedo_direct(:,:,:,1),direct_dn_below) |
|---|
| 1388 | |
|---|
| 1389 | if (config%do_clear) then |
|---|
| 1390 | flux_dn_clear = 0.0_jprb |
|---|
| 1391 | direct_dn_clear(:) = incoming_sw(:,jcol) |
|---|
| 1392 | flux_up_clear = direct_dn_clear*total_albedo_clear_direct(:,1) |
|---|
| 1393 | end if |
|---|
| 1394 | |
|---|
| 1395 | ! Store the TOA broadband fluxes |
|---|
| 1396 | flux%sw_up(jcol,1) = sum(sum(flux_up_above,1)) |
|---|
| 1397 | flux%sw_dn(jcol,1) = mu0 * sum(direct_dn_clear(:)) |
|---|
| 1398 | if (allocated(flux%sw_dn_direct)) then |
|---|
| 1399 | flux%sw_dn_direct(jcol,1) = flux%sw_dn(jcol,1) |
|---|
| 1400 | end if |
|---|
| 1401 | if (config%do_clear) then |
|---|
| 1402 | flux%sw_up_clear(jcol,1) = sum(flux_up_clear) |
|---|
| 1403 | flux%sw_dn_clear(jcol,1) = flux%sw_dn(jcol,1) |
|---|
| 1404 | if (allocated(flux%sw_dn_direct_clear)) then |
|---|
| 1405 | flux%sw_dn_direct_clear(jcol,1) = flux%sw_dn_clear(jcol,1) |
|---|
| 1406 | end if |
|---|
| 1407 | end if |
|---|
| 1408 | |
|---|
| 1409 | ! Save the spectral fluxes if required |
|---|
| 1410 | if (config%do_save_spectral_flux) then |
|---|
| 1411 | call indexed_sum(sum(flux_up_above(:,:),2), & |
|---|
| 1412 | & config%i_spec_from_reordered_g_sw, & |
|---|
| 1413 | & flux%sw_up_band(:,jcol,1)) |
|---|
| 1414 | call indexed_sum(sum(direct_dn_below(:,:),2), & |
|---|
| 1415 | & config%i_spec_from_reordered_g_sw, & |
|---|
| 1416 | & flux%sw_dn_band(:,jcol,1)) |
|---|
| 1417 | flux%sw_dn_band(:,jcol,1) = mu0 * flux%sw_dn_band(:,jcol,1) |
|---|
| 1418 | if (allocated(flux%sw_dn_direct_band)) then |
|---|
| 1419 | flux%sw_dn_direct_band(:,jcol,1) = flux%sw_dn_band(:,jcol,1) |
|---|
| 1420 | end if |
|---|
| 1421 | if (config%do_clear) then |
|---|
| 1422 | flux%sw_dn_clear_band(:,jcol,1) = flux%sw_dn_band(:,jcol,1) |
|---|
| 1423 | call indexed_sum(flux_up_clear, & |
|---|
| 1424 | & config%i_spec_from_reordered_g_sw, & |
|---|
| 1425 | & flux%sw_up_clear_band(:,jcol,1)) |
|---|
| 1426 | if (allocated(flux%sw_dn_direct_clear_band)) then |
|---|
| 1427 | flux%sw_dn_direct_clear_band(:,jcol,1) & |
|---|
| 1428 | & = flux%sw_dn_clear_band(:,jcol,1) |
|---|
| 1429 | end if |
|---|
| 1430 | end if |
|---|
| 1431 | end if |
|---|
| 1432 | |
|---|
| 1433 | ! Final loop back down through the atmosphere to compute fluxes |
|---|
| 1434 | do jlev = 1,nlev |
|---|
| 1435 | |
|---|
| 1436 | #ifdef PRINT_ENTRAPMENT_DATA |
|---|
| 1437 | if (config%i_3d_sw_entrapment == IEntrapmentExplicitNonFractal & |
|---|
| 1438 | & .or. config%i_3d_sw_entrapment == IEntrapmentExplicit) then |
|---|
| 1439 | ! Save downwelling direct and diffuse fluxes at the top of |
|---|
| 1440 | ! layer "jlev" in each of the regions of layer "jlev" |
|---|
| 1441 | if (nreg == 2) then |
|---|
| 1442 | write(102,'(i4,i4,4e14.6)') jcol, jlev, direct_dn_below(1,:), flux_dn_below(1,:) |
|---|
| 1443 | else |
|---|
| 1444 | write(102,'(i4,i4,6e14.6)') jcol, jlev, direct_dn_below(1,1:3), flux_dn_below(1,1:3) |
|---|
| 1445 | end if |
|---|
| 1446 | end if |
|---|
| 1447 | #endif |
|---|
| 1448 | |
|---|
| 1449 | ! Compute the solar downwelling "source" at the base of the |
|---|
| 1450 | ! layer due to scattering of the direct beam within it |
|---|
| 1451 | if (config%do_clear) then |
|---|
| 1452 | source_dn_clear = trans_dir_diff_clear(:,jlev)*direct_dn_clear |
|---|
| 1453 | end if |
|---|
| 1454 | source_dn(:,:) = mat_x_vec(ng,ng,nreg,trans_dir_diff(:,:,:,jlev),direct_dn_below, & |
|---|
| 1455 | & is_clear_sky_layer(jlev)) |
|---|
| 1456 | |
|---|
| 1457 | ! Compute direct downwelling flux in each region at base of |
|---|
| 1458 | ! current layer |
|---|
| 1459 | if (config%do_clear) then |
|---|
| 1460 | direct_dn_clear = trans_dir_dir_clear(:,jlev)*direct_dn_clear |
|---|
| 1461 | end if |
|---|
| 1462 | direct_dn_above = mat_x_vec(ng,ng,nreg,trans_dir_dir(:,:,:,jlev),direct_dn_below, & |
|---|
| 1463 | & is_clear_sky_layer(jlev)) |
|---|
| 1464 | |
|---|
| 1465 | ! Integrate downwelling direct flux across spectrum and |
|---|
| 1466 | ! regions, and store (the diffuse part will be added later) |
|---|
| 1467 | flux%sw_dn(jcol,jlev+1) = mu0 * sum(sum(direct_dn_above,1)) |
|---|
| 1468 | if (allocated(flux%sw_dn_direct)) then |
|---|
| 1469 | flux%sw_dn_direct(jcol,jlev+1) = flux%sw_dn(jcol,jlev+1) |
|---|
| 1470 | end if |
|---|
| 1471 | if (config%do_clear) then |
|---|
| 1472 | flux%sw_dn_clear(jcol,jlev+1) = mu0 * sum(direct_dn_clear) |
|---|
| 1473 | if (allocated(flux%sw_dn_direct_clear)) then |
|---|
| 1474 | flux%sw_dn_direct_clear(jcol,jlev+1) & |
|---|
| 1475 | & = flux%sw_dn_clear(jcol,jlev+1) |
|---|
| 1476 | end if |
|---|
| 1477 | end if |
|---|
| 1478 | |
|---|
| 1479 | if (config%do_save_spectral_flux) then |
|---|
| 1480 | call indexed_sum(sum(direct_dn_above,2), & |
|---|
| 1481 | & config%i_spec_from_reordered_g_sw, & |
|---|
| 1482 | & flux%sw_dn_band(:,jcol,jlev+1)) |
|---|
| 1483 | flux%sw_dn_band(:,jcol,jlev+1) = mu0 * flux%sw_dn_band(:,jcol,jlev+1) |
|---|
| 1484 | |
|---|
| 1485 | if (allocated(flux%sw_dn_direct_band)) then |
|---|
| 1486 | flux%sw_dn_direct_band(:,jcol,jlev+1) & |
|---|
| 1487 | & = flux%sw_dn_band(:,jcol,jlev+1) |
|---|
| 1488 | end if |
|---|
| 1489 | if (config%do_clear) then |
|---|
| 1490 | call indexed_sum(direct_dn_clear, & |
|---|
| 1491 | & config%i_spec_from_reordered_g_sw, & |
|---|
| 1492 | & flux%sw_dn_clear_band(:,jcol,jlev+1)) |
|---|
| 1493 | flux%sw_dn_clear_band(:,jcol,jlev+1) = mu0 & |
|---|
| 1494 | & * flux%sw_dn_clear_band(:,jcol,jlev+1) |
|---|
| 1495 | if (allocated(flux%sw_dn_direct_clear_band)) then |
|---|
| 1496 | flux%sw_dn_direct_clear_band(:,jcol,jlev+1) & |
|---|
| 1497 | & = flux%sw_dn_clear_band(:,jcol,jlev+1) |
|---|
| 1498 | end if |
|---|
| 1499 | end if |
|---|
| 1500 | end if |
|---|
| 1501 | |
|---|
| 1502 | if (config%do_clear) then |
|---|
| 1503 | ! Scalar operations for clear-sky fluxes |
|---|
| 1504 | flux_dn_clear(:) = (trans_clear(:,jlev)*flux_dn_clear(:) & |
|---|
| 1505 | & + ref_clear(:,jlev)*total_albedo_clear_direct(:,jlev+1)*direct_dn_clear & |
|---|
| 1506 | & + source_dn_clear) & |
|---|
| 1507 | & / (1.0_jprb - ref_clear(:,jlev)*total_albedo_clear(:,jlev+1)) |
|---|
| 1508 | flux_up_clear(:) = total_albedo_clear_direct(:,jlev+1)*direct_dn_clear & |
|---|
| 1509 | & + total_albedo_clear(:,jlev+1)*flux_dn_clear |
|---|
| 1510 | end if |
|---|
| 1511 | |
|---|
| 1512 | if (is_clear_sky_layer(jlev)) then |
|---|
| 1513 | ! Scalar operations for clear-sky layer |
|---|
| 1514 | flux_dn_above(:,1) = (transmittance(:,1,1,jlev)*flux_dn_below(:,1) & |
|---|
| 1515 | & + reflectance(:,1,1,jlev)*total_albedo_direct(:,1,1,jlev+1)*direct_dn_above(:,1) & |
|---|
| 1516 | & + source_dn(:,1)) & |
|---|
| 1517 | & / (1.0_jprb - reflectance(:,1,1,jlev)*total_albedo(:,1,1,jlev+1)) |
|---|
| 1518 | flux_dn_above(:,2:nreg) = 0.0_jprb |
|---|
| 1519 | flux_up_above(:,1) = total_albedo_direct(:,1,1,jlev+1)*direct_dn_above(:,1) & |
|---|
| 1520 | & + total_albedo(:,1,1,jlev+1)*flux_dn_above(:,1) |
|---|
| 1521 | flux_up_above(:,2:nreg) = 0.0_jprb |
|---|
| 1522 | else |
|---|
| 1523 | ! Matrix operations for cloudy layer |
|---|
| 1524 | denominator = identity_minus_mat_x_mat(ng,ng,nreg,reflectance(:,:,:,jlev), & |
|---|
| 1525 | & total_albedo(:,:,:,jlev+1)) |
|---|
| 1526 | total_source = mat_x_vec(ng,ng,nreg,total_albedo_direct(:,:,:,jlev+1),direct_dn_above) |
|---|
| 1527 | |
|---|
| 1528 | flux_dn_above = solve_vec(ng,ng,nreg,denominator, & |
|---|
| 1529 | & mat_x_vec(ng,ng,nreg,transmittance(:,:,:,jlev),flux_dn_below) & |
|---|
| 1530 | & + mat_x_vec(ng,ng,nreg,reflectance(:,:,:,jlev), total_source(:,:)) & |
|---|
| 1531 | & + source_dn(:,:)) |
|---|
| 1532 | flux_up_above = mat_x_vec(ng,ng,nreg,total_albedo(:,:,:,jlev+1), & |
|---|
| 1533 | & flux_dn_above) + total_source(:,:) |
|---|
| 1534 | end if |
|---|
| 1535 | |
|---|
| 1536 | ! Account for overlap rules in translating fluxes just above |
|---|
| 1537 | ! a layer interface to the values just below |
|---|
| 1538 | if (is_clear_sky_layer(jlev) .and. is_clear_sky_layer(jlev+1)) then |
|---|
| 1539 | ! Regions in current layer map directly on to regions in |
|---|
| 1540 | ! layer below |
|---|
| 1541 | flux_dn_below = flux_dn_above |
|---|
| 1542 | direct_dn_below = direct_dn_above |
|---|
| 1543 | else |
|---|
| 1544 | ! Apply downward overlap matrix to compute direct |
|---|
| 1545 | ! downwelling flux entering the top of each region in the |
|---|
| 1546 | ! layer below |
|---|
| 1547 | flux_dn_below = singlemat_x_vec(ng,ng,nreg,v_matrix(:,:,jlev+1,jcol), & |
|---|
| 1548 | & flux_dn_above) |
|---|
| 1549 | direct_dn_below = singlemat_x_vec(ng,ng,nreg,v_matrix(:,:,jlev+1,jcol), & |
|---|
| 1550 | & direct_dn_above) |
|---|
| 1551 | end if |
|---|
| 1552 | |
|---|
| 1553 | ! Store the broadband fluxes |
|---|
| 1554 | flux%sw_up(jcol,jlev+1) = sum(sum(flux_up_above,1)) |
|---|
| 1555 | flux%sw_dn(jcol,jlev+1) & |
|---|
| 1556 | & = flux%sw_dn(jcol,jlev+1) + sum(sum(flux_dn_above,1)) |
|---|
| 1557 | if (config%do_clear) then |
|---|
| 1558 | flux%sw_up_clear(jcol,jlev+1) = sum(flux_up_clear) |
|---|
| 1559 | flux%sw_dn_clear(jcol,jlev+1) & |
|---|
| 1560 | & = flux%sw_dn_clear(jcol,jlev+1) + sum(flux_dn_clear) |
|---|
| 1561 | end if |
|---|
| 1562 | |
|---|
| 1563 | ! Save the spectral fluxes if required |
|---|
| 1564 | if (config%do_save_spectral_flux) then |
|---|
| 1565 | call indexed_sum(sum(flux_up_above,2), & |
|---|
| 1566 | & config%i_spec_from_reordered_g_sw, & |
|---|
| 1567 | & flux%sw_up_band(:,jcol,jlev+1)) |
|---|
| 1568 | call add_indexed_sum(sum(flux_dn_above,2), & |
|---|
| 1569 | & config%i_spec_from_reordered_g_sw, & |
|---|
| 1570 | & flux%sw_dn_band(:,jcol,jlev+1)) |
|---|
| 1571 | if (config%do_clear) then |
|---|
| 1572 | call indexed_sum(flux_up_clear, & |
|---|
| 1573 | & config%i_spec_from_reordered_g_sw, & |
|---|
| 1574 | & flux%sw_up_clear_band(:,jcol,jlev+1)) |
|---|
| 1575 | call add_indexed_sum(flux_dn_clear, & |
|---|
| 1576 | & config%i_spec_from_reordered_g_sw, & |
|---|
| 1577 | & flux%sw_dn_clear_band(:,jcol,jlev+1)) |
|---|
| 1578 | end if |
|---|
| 1579 | end if |
|---|
| 1580 | |
|---|
| 1581 | end do ! Final loop over levels |
|---|
| 1582 | |
|---|
| 1583 | ! Store surface spectral fluxes, if required (after the end of |
|---|
| 1584 | ! the final loop over levels, the current values of these arrays |
|---|
| 1585 | ! will be the surface values) |
|---|
| 1586 | flux%sw_dn_diffuse_surf_g(:,jcol) = sum(flux_dn_above,2) |
|---|
| 1587 | flux%sw_dn_direct_surf_g(:,jcol) = mu0 * sum(direct_dn_above,2) |
|---|
| 1588 | if (config%do_clear) then |
|---|
| 1589 | flux%sw_dn_diffuse_surf_clear_g(:,jcol) = flux_dn_clear |
|---|
| 1590 | flux%sw_dn_direct_surf_clear_g(:,jcol) = mu0 * direct_dn_clear |
|---|
| 1591 | end if |
|---|
| 1592 | |
|---|
| 1593 | end do ! Loop over columns |
|---|
| 1594 | if (config%iverbose >= 3) then |
|---|
| 1595 | write(nulout,*) |
|---|
| 1596 | end if |
|---|
| 1597 | |
|---|
| 1598 | ! Report number of calls to each method of solving single-layer |
|---|
| 1599 | ! two-stream equations |
|---|
| 1600 | if (config%iverbose >= 4) then |
|---|
| 1601 | write(nulout,'(a,i0)') ' Matrix-exponential calls: ', n_calls_expm |
|---|
| 1602 | write(nulout,'(a,i0)') ' Meador-Weaver calls: ', n_calls_meador_weaver |
|---|
| 1603 | end if |
|---|
| 1604 | |
|---|
| 1605 | if (lhook) call dr_hook('radiation_spartacus_sw:solver_spartacus_sw',1,hook_handle) |
|---|
| 1606 | |
|---|
| 1607 | end subroutine solver_spartacus_sw |
|---|
| 1608 | |
|---|
| 1609 | |
|---|
| 1610 | ! Step the horizontal migration distances from the base of a layer |
|---|
| 1611 | ! to the top, accounting for the extra distance travelled within the |
|---|
| 1612 | ! layer |
|---|
| 1613 | subroutine step_migrations(ng, nreg, cloud_frac, & |
|---|
| 1614 | & layer_depth, tan_diffuse_angle_3d, tan_sza, & |
|---|
| 1615 | & reflectance, transmittance, ref_dir, trans_dir_dir, & |
|---|
| 1616 | & trans_dir_diff, total_albedo_diff, total_albedo_dir, & |
|---|
| 1617 | & x_diffuse, x_direct) |
|---|
| 1618 | |
|---|
| 1619 | use parkind1, only : jprb |
|---|
| 1620 | |
|---|
| 1621 | implicit none |
|---|
| 1622 | |
|---|
| 1623 | ! Inputs |
|---|
| 1624 | |
|---|
| 1625 | ! Number of g points and regions |
|---|
| 1626 | integer, intent(in) :: ng, nreg |
|---|
| 1627 | ! Cloud fraction |
|---|
| 1628 | real(jprb), intent(in) :: cloud_frac |
|---|
| 1629 | ! Layer depth (m), tangent of diffuse zenith angle and tangent of |
|---|
| 1630 | ! solar zenith angle |
|---|
| 1631 | real(jprb), intent(in) :: layer_depth, tan_diffuse_angle_3d, tan_sza |
|---|
| 1632 | ! Reflectance and transmittance to diffuse downwelling radiation |
|---|
| 1633 | real(jprb), intent(in), dimension(ng, nreg, nreg) :: reflectance, transmittance |
|---|
| 1634 | ! Reflectance and transmittance to direct downwelling radiation |
|---|
| 1635 | real(jprb), intent(in), dimension(ng, nreg, nreg) :: ref_dir, trans_dir_dir |
|---|
| 1636 | ! Transmittance involving direct entering a layer from the top and |
|---|
| 1637 | ! diffuse leaving from the bottom |
|---|
| 1638 | real(jprb), intent(in), dimension(ng, nreg, nreg) :: trans_dir_diff |
|---|
| 1639 | |
|---|
| 1640 | ! Total albedo of direct and diffuse radiation of the atmosphere |
|---|
| 1641 | ! below the layer in question |
|---|
| 1642 | real(jprb), intent(in), dimension(ng, nreg, nreg) & |
|---|
| 1643 | & :: total_albedo_diff, total_albedo_dir |
|---|
| 1644 | |
|---|
| 1645 | ! Inputs/outputs |
|---|
| 1646 | |
|---|
| 1647 | ! Horizontal migration distance (m) of reflected light |
|---|
| 1648 | real(jprb), intent(inout), dimension(ng, nreg) :: x_diffuse, x_direct |
|---|
| 1649 | |
|---|
| 1650 | ! Local variables |
|---|
| 1651 | |
|---|
| 1652 | ! Top albedo, i.e. the albedo of the top of the layer assuming no |
|---|
| 1653 | ! lateral transport |
|---|
| 1654 | real(jprb), dimension(ng) :: top_albedo |
|---|
| 1655 | |
|---|
| 1656 | ! Multiple-scattering amplitude enhancement |
|---|
| 1657 | real(jprb), dimension(ng) :: ms_enhancement |
|---|
| 1658 | |
|---|
| 1659 | ! Multiple-scattering distance enhancement |
|---|
| 1660 | real(jprb), dimension(ng) :: x_enhancement |
|---|
| 1661 | |
|---|
| 1662 | |
|---|
| 1663 | real(jprb) :: x_layer_diffuse, x_layer_direct |
|---|
| 1664 | integer :: jreg, istartreg, iendreg |
|---|
| 1665 | |
|---|
| 1666 | istartreg = 1 |
|---|
| 1667 | iendreg = nreg |
|---|
| 1668 | |
|---|
| 1669 | if (cloud_frac <= 0.0_jprb) then |
|---|
| 1670 | ! Clear-sky layer: don't waste time on cloudy regions |
|---|
| 1671 | iendreg = 1 |
|---|
| 1672 | else if (cloud_frac >= 1.0_jprb) then |
|---|
| 1673 | ! Overcast layer: don't waste time on clear region |
|---|
| 1674 | istartreg = 2 |
|---|
| 1675 | end if |
|---|
| 1676 | |
|---|
| 1677 | ! This is the mean horizontal distance travelled by diffuse |
|---|
| 1678 | ! radiation that travels from the top of a layer to the centre and |
|---|
| 1679 | ! is then scattered back up and out |
|---|
| 1680 | x_layer_diffuse = layer_depth * tan_diffuse_angle_3d/sqrt(2.0_jprb) |
|---|
| 1681 | |
|---|
| 1682 | ! This is the mean horizontal distance travelled by direct |
|---|
| 1683 | ! radiation that travels from the top of a layer to the centre and |
|---|
| 1684 | ! is then scattered back up and out |
|---|
| 1685 | x_layer_direct = layer_depth * sqrt(tan_sza*tan_sza & |
|---|
| 1686 | & + tan_diffuse_angle_3d*tan_diffuse_angle_3d) * 0.5_jprb |
|---|
| 1687 | |
|---|
| 1688 | do jreg = istartreg,iendreg |
|---|
| 1689 | ! Geometric series enhancement due to multiple scattering: the |
|---|
| 1690 | ! amplitude enhancement is equal to the limit of |
|---|
| 1691 | ! T*[1+RA+(RA)^2+(RA)^3+...] |
|---|
| 1692 | ms_enhancement = transmittance(:,jreg,jreg) & |
|---|
| 1693 | & / (1.0_jprb - reflectance(:,jreg,jreg)*total_albedo_diff(:,jreg,jreg)) |
|---|
| 1694 | ! ...and the distance enhancement is approximately equal to the |
|---|
| 1695 | ! limit of T*[1+sqrt(2)*RA+sqrt(3)*(RA)^2+sqrt(4)*(RA)^3+...] |
|---|
| 1696 | x_enhancement = (1.0_jprb - reflectance(:,jreg,jreg)*total_albedo_diff(:,jreg,jreg))**(-1.5_jprb) |
|---|
| 1697 | |
|---|
| 1698 | ! Horizontal migration of direct downwelling radiation |
|---|
| 1699 | top_albedo = max(1.0e-8_jprb, ref_dir(:,jreg,jreg) + ms_enhancement & |
|---|
| 1700 | & * (trans_dir_diff(:,jreg,jreg)*total_albedo_diff(:,jreg,jreg) & |
|---|
| 1701 | & +trans_dir_dir(:,jreg,jreg)*total_albedo_dir(:,jreg,jreg))) |
|---|
| 1702 | ! The following is approximate and has been found to |
|---|
| 1703 | ! occasionally go negative |
|---|
| 1704 | x_direct(:,jreg) = max(0.0_jprb, x_layer_direct & |
|---|
| 1705 | & + ((trans_dir_diff(:,jreg,jreg)*total_albedo_diff(:,jreg,jreg)*x_enhancement & |
|---|
| 1706 | & +trans_dir_dir(:,jreg,jreg)*total_albedo_dir(:,jreg,jreg)*(x_enhancement-1.0_jprb)) & |
|---|
| 1707 | & *(x_diffuse(:,jreg)+x_layer_diffuse) & |
|---|
| 1708 | & +trans_dir_dir(:,jreg,jreg)*total_albedo_dir(:,jreg,jreg) & |
|---|
| 1709 | & *(x_direct(:,jreg)+x_layer_direct)) & |
|---|
| 1710 | & * transmittance(:,jreg,jreg) / top_albedo) |
|---|
| 1711 | |
|---|
| 1712 | ! Horizontal migration of diffuse downwelling radiation |
|---|
| 1713 | top_albedo = max(1.0e-8_jprb, reflectance(:,jreg,jreg) & |
|---|
| 1714 | & + ms_enhancement*transmittance(:,jreg,jreg)*total_albedo_diff(:,jreg,jreg)) |
|---|
| 1715 | x_diffuse(:,jreg) = x_layer_diffuse + x_enhancement*total_albedo_diff(:,jreg,jreg) & |
|---|
| 1716 | & *(transmittance(:,jreg,jreg)*transmittance(:,jreg,jreg)) & |
|---|
| 1717 | & * (x_diffuse(:,jreg) + x_layer_diffuse) / top_albedo |
|---|
| 1718 | |
|---|
| 1719 | end do |
|---|
| 1720 | if (iendreg < nreg) then |
|---|
| 1721 | x_diffuse(:,iendreg+1:nreg) = 0.0_jprb |
|---|
| 1722 | x_direct(:,iendreg+1:nreg) = 0.0_jprb |
|---|
| 1723 | else if (istartreg == 2) then |
|---|
| 1724 | x_diffuse(:,1) = 0.0_jprb |
|---|
| 1725 | x_direct(:,1) = 0.0_jprb |
|---|
| 1726 | end if |
|---|
| 1727 | |
|---|
| 1728 | end subroutine step_migrations |
|---|
| 1729 | |
|---|
| 1730 | end module radiation_spartacus_sw |
|---|
