source: LMDZ6/branches/blowing_snow/libf/phylmd/ecrad/radiation_spartacus_lw.F90 @ 5411

Last change on this file since 5411 was 3908, checked in by idelkadi, 4 years ago

Online implementation of the radiative transfer code ECRAD in the LMDZ model.

  • Inclusion of the ecrad directory containing the sources of the ECRAD code
    • interface routine : radiation_scheme.F90
  • Adaptation of compilation scripts :
    • compilation under CPP key CPP_ECRAD
    • compilation with option "-rad ecard" or "-ecard true"
    • The "-rad old/rtm/ecran" build option will need to replace the "-rrtm true" and "-ecrad true" options in the future.
  • Runing LMDZ simulations with ecrad, you need :
    • logical key iflag_rrtm = 2 in physiq.def
    • namelist_ecrad (DefLists?)
    • the directory "data" containing the configuration files is temporarily placed in ../libfphylmd/ecrad/
  • Compilation and execution are tested in the 1D case. The repository under svn would allow to continue the implementation work: tests, verification of the results, ...
File size: 48.7 KB
Line 
1! radiation_spartacus_lw.F90 - SPARTACUS longwave 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 emission/albedo rather than planck/emissivity
17!   2017-04-22  R. Hogan  Store surface fluxes at all g-points
18!   2018-09-03  R. Hogan  Security via min_cloud_effective_size
19!   2018-10-08  R. Hogan  Call calc_region_properties
20!   2019-01-12  R. Hogan  Use inv_inhom_effective_size if allocated
21
22module radiation_spartacus_lw
23
24  public
25
26contains
27
28  ! Small routine for scaling cloud optical depth in the cloudy
29  ! regions
30#include "radiation_optical_depth_scaling.h"
31
32  ! This module contains just one subroutine, the longwave solver
33  ! using the Speedy Algorithm for Radiative Transfer through Cloud
34  ! Sides (SPARTACUS), which can represent 3D effects using a matrix
35  ! form of the two-stream equations.
36  !
37  ! Sections:
38  !   1: Prepare general variables and arrays
39  !   2: Prepare column-specific variables and arrays
40  !   3: First loop over layers
41  !     3.1: Layer-specific initialization
42  !     3.2: Compute gamma variables
43  !       3.2a: Clear-sky case
44  !       3.2b: Cloudy case
45  !     3.3: Compute reflection, transmission and emission
46  !       3.3a: g-points with 3D effects
47  !       3.3b: g-points without 3D effects
48  !   4: Compute total sources and albedos
49  !   5: Compute fluxes
50  subroutine solver_spartacus_lw(nlev,istartcol,iendcol, &
51       &  config, thermodynamics, cloud, &
52       &  od, ssa, g, od_cloud, ssa_cloud, g_cloud, planck_hl, &
53       &  emission, albedo, &
54       &  flux)
55
56    use parkind1,                 only : jprb
57    use yomhook,                  only : lhook, dr_hook
58
59    use radiation_io,             only : nulout
60    use radiation_config,         only : config_type, IPdfShapeGamma
61    use radiation_thermodynamics, only : thermodynamics_type
62    use radiation_cloud,          only : cloud_type
63    use radiation_regions,        only : calc_region_properties
64    use radiation_overlap,        only : calc_overlap_matrices
65    use radiation_flux,           only : flux_type, indexed_sum
66    use radiation_matrix
67    use radiation_two_stream,     only : calc_two_stream_gammas_lw, &
68         calc_reflectance_transmittance_lw, LwDiffusivity
69    use radiation_lw_derivatives, only : calc_lw_derivatives_matrix
70    use radiation_constants,      only : Pi, GasConstantDryAir, &
71         AccelDueToGravity
72
73    implicit none
74
75    ! Inputs
76    integer, intent(in) :: nlev               ! number of model levels
77    integer, intent(in) :: istartcol, iendcol ! range of columns to process
78    type(config_type), intent(in)        :: config
79    type(thermodynamics_type),intent(in) :: thermodynamics
80    type(cloud_type),   intent(in)       :: cloud
81
82    ! Gas and aerosol optical depth of each layer at each longwave
83    ! g-point
84    real(jprb), intent(in), dimension(config%n_g_lw,nlev,istartcol:iendcol) :: od
85
86    ! Gas and aerosol single-scattering albedo and asymmetry factor,
87    ! only if longwave scattering by aerosols is to be represented
88    real(jprb), intent(in), &
89         &  dimension(config%n_g_lw_if_scattering,nlev,istartcol:iendcol) :: ssa, g
90
91    ! Cloud and precipitation optical depth of each layer in each
92    ! longwave band
93    real(jprb), intent(in) :: od_cloud(config%n_bands_lw,nlev,istartcol:iendcol)
94
95    ! Cloud and precipitation single-scattering albedo and asymmetry
96    ! factor, only if longwave scattering by clouds is to be
97    ! represented
98    real(jprb), intent(in), dimension(config%n_bands_lw_if_scattering, &
99         &                            nlev,istartcol:iendcol) :: ssa_cloud, g_cloud
100
101    ! Planck function (emitted flux from a black body) at half levels
102    real(jprb), intent(in), dimension(config%n_g_lw,nlev+1,istartcol:iendcol) &
103         &  :: planck_hl
104
105    ! Emission (Planck*emissivity) and albedo (1-emissivity) at the
106    ! surface at each longwave g-point
107    real(jprb), intent(in), dimension(config%n_g_lw, istartcol:iendcol) &
108         &  :: emission, albedo
109
110    ! Output
111    type(flux_type),          intent(inout):: flux
112
113    integer :: nreg, ng
114    integer :: nregActive ! =1 in clear layer, =nreg in a cloudy layer
115    integer :: jcol, jlev, jg, jreg, iband, jreg2
116    integer :: ng3D ! Number of g-points with small enough gas optical
117                    ! depth that 3D effects need to be represented
118
119    ! Ratio of gas constant for dry air to acceleration due to gravity
120    real(jprb), parameter &
121         &  :: R_over_g = GasConstantDryAir / AccelDueToGravity
122
123    ! Used in computing rates of lateral radiation transfer
124    real(jprb), parameter :: four_over_pi = 4.0_jprb / Pi
125
126    ! The tangent of the effective zenith angle for diffuse radiation
127    ! that is appropriate for 3D transport
128    real(jprb), parameter :: tan_diffuse_angle_3d = Pi * 0.5_jprb
129    ! Incorrect value of tand(53) used by Hogan & Shonk (2013)
130    !    real(jprb), parameter :: tan_diffuse_angle_3d = 1.32704482162041
131
132    ! Optical depth, single scattering albedo and asymmetry factor in
133    ! each region at each g-point
134    real(jprb), dimension(config%n_g_lw, config%nregions) &
135         &  :: od_region, ssa_region, g_region
136
137    ! Scattering optical depths of gases and clouds
138    real(jprb) :: scat_od, scat_od_cloud
139
140    ! The area fractions of each region
141    real(jprb) :: region_fracs(1:config%nregions,nlev,istartcol:iendcol)
142
143    ! The scaling used for the optical depth in the cloudy regions
144    real(jprb) :: od_scaling(2:config%nregions,nlev,istartcol:iendcol)
145
146    ! The length of the interface between regions per unit area of
147    ! gridbox, equal to L_diff^ab in Hogan and Shonk (2013). This is
148    ! actually the effective length oriented to a photon with random
149    ! azimuth angle. The three values are the edge length between
150    ! regions a-b, b-c and a-c.
151    real(jprb) :: edge_length(3)
152
153    ! Element i,j gives the rate of 3D transfer of diffuse
154    ! radiation from region i to region j, multiplied by the thickness
155    ! of the layer in m
156    real(jprb) :: transfer_rate(config%nregions,config%nregions)
157
158    ! Directional overlap matrices defined at all layer interfaces
159    ! including top-of-atmosphere and the surface
160    real(jprb), dimension(config%nregions,config%nregions,nlev+1,istartcol:iendcol) &
161         &  :: u_matrix, v_matrix
162
163    ! Two-stream variables
164    real(jprb), dimension(config%n_g_lw, config%nregions) &
165         &  :: gamma1, gamma2
166
167    ! Matrix Gamma multiplied by the layer thickness z1, so units
168    ! metres.  After calling expm, this array contains the matrix
169    ! exponential of the original.
170    real(jprb), dimension(config%n_g_lw, 2*config%nregions, &
171         &  2*config%nregions) :: Gamma_z1
172
173    ! Diffuse reflection and transmission matrices of each layer
174    real(jprb), dimension(config%n_g_lw, config%nregions, &
175         &  config%nregions, nlev) :: reflectance, transmittance
176
177    ! Clear-sky diffuse reflection and transmission matrices of each
178    ! layer
179    real(jprb), dimension(config%n_g_lw, nlev) :: ref_clear, trans_clear
180
181    ! If the Planck term at the top of a layer in each region is the
182    ! vector b0 then the following is vector [-b0; b0]*dz
183    real(jprb), dimension(config%n_g_lw,2*config%nregions) :: planck_top
184
185    ! The difference between the Planck terms at the bottom and top of
186    ! a layer, doubled as with planck_top; in terms of b' in the
187    ! paper, planck_diff is [-b'; b']*dz*dz
188    real(jprb), dimension(config%n_g_lw,2*config%nregions) :: planck_diff
189
190    ! Parts of the particular solution associated with the
191    ! inhomogeneous ODE. In terms of quantities from the paper,
192    ! solution0 is [c0;d0] and solution_diff is [c';d']*dz.
193    real(jprb), dimension(config%n_g_lw,2*config%nregions) &
194         &  :: solution0, solution_diff
195
196    ! Used for computing the Planck emission per layer
197    real(jprb), dimension(config%n_g_lw,config%nregions) &
198         &  :: tmp_vectors
199
200    ! The fluxes upwelling from the top of a layer (source_up) and
201    ! downwelling from the bottom of the layer (source_dn) due to
202    ! emission
203    real(jprb), dimension(config%n_g_lw, config%nregions, nlev) &
204         &  :: source_up, source_dn
205    ! ...clear-sky equivalents
206    real(jprb), dimension(config%n_g_lw, nlev) &
207         &  :: source_up_clear, source_dn_clear
208
209    ! Upwelling radiation just above a layer interface due to emission
210    ! below that interface, where level index = 1 corresponds to the
211    ! top-of-atmosphere
212    real(jprb), dimension(config%n_g_lw, config%nregions, nlev+1) &
213         &  :: total_source
214    ! ...clear-sky equivalent
215    real(jprb), dimension(config%n_g_lw, nlev+1) :: total_source_clear
216
217    ! As total_source, but just below a layer interface
218    real(jprb), dimension(config%n_g_lw, config%nregions) &
219         &  :: total_source_below
220
221    ! Total albedo of the atmosphere/surface just above a layer
222    ! interface with respect to downwelling diffuse radiation at that
223    ! interface, where level index = 1 corresponds to the
224    ! top-of-atmosphere
225    real(jprb), dimension(config%n_g_lw, config%nregions, &
226         &  config%nregions, nlev+1) :: total_albedo
227    ! ...clear-sky equivalent
228    real(jprb), dimension(config%n_g_lw, nlev+1) :: total_albedo_clear
229
230    ! As total_albedo, but just below a layer interface
231    real(jprb), dimension(config%n_g_lw, config%nregions, config%nregions) &
232         &  :: total_albedo_below
233
234    ! The following is used to store matrices of the form I-A*B that
235    ! are used on the denominator of some expressions
236    real(jprb), dimension(config%n_g_lw, config%nregions, config%nregions) &
237         &  :: denominator
238    ! Clear-sky equivalent, but actually its reciprocal to replace
239    ! some divisions by multiplications
240    real(jprb), dimension(config%n_g_lw) :: inv_denom_scalar
241
242    ! Layer depth (m)
243    real(jprb) :: dz
244
245    ! Upwelling and downwelling fluxes above/below layer interfaces
246    real(jprb), dimension(config%n_g_lw, config%nregions) &
247         &  :: flux_up_above, flux_dn_above, flux_dn_below
248    ! Clear-sky upwelling and downwelling fluxes (which don't
249    ! distinguish between whether they are above/below a layer
250    ! interface)
251    real(jprb), dimension(config%n_g_lw) :: flux_up_clear, flux_dn_clear
252
253    ! Parameterization of internal flux distribution in a cloud that
254    ! can lead to the flux just about to exit a cloud being different
255    ! from the mean in-cloud value.  "lateral_od" is the typical
256    ! absorption optical depth through the cloud in a horizontal
257    ! direction.
258    real(jprb), dimension(config%n_g_lw) :: lateral_od, sqrt_1_minus_ssa
259    real(jprb), dimension(config%n_g_lw) :: side_emiss_thick, side_emiss
260
261    ! In the optically thin limit, the flux near the edge is greater
262    ! than that in the interior
263    real(jprb), parameter :: side_emiss_thin = 1.4107
264
265    real(jprb) :: aspect_ratio
266
267    ! Keep a count of the number of calls to the two ways of computing
268    ! reflectance/transmittance matrices
269    integer :: n_calls_expm, n_calls_meador_weaver
270
271    ! Identify clear-sky layers, with pseudo layers for outer space
272    ! and below the ground, both treated as single-region clear skies
273    logical :: is_clear_sky_layer(0:nlev+1)
274
275    real(jprb) :: hook_handle
276
277    if (lhook) call dr_hook('radiation_spartacus_lw:solver_spartacus_lw',0,hook_handle)
278
279    ! --------------------------------------------------------
280    ! Section 1: Prepare general variables and arrays
281    ! --------------------------------------------------------
282
283    ! Copy array dimensions to local variables for convenience
284    nreg = config%nregions
285    ng = config%n_g_lw
286
287    ! Reset count of number of calls to the two ways to compute
288    ! reflectance/transmission matrices
289    n_calls_expm = 0
290    n_calls_meador_weaver = 0
291
292    ! Initialize dz to avoid compiler warning
293    dz = 1.0_jprb
294
295    ! Compute the wavelength-independent region fractions and
296    ! optical-depth scalings
297    call calc_region_properties(nlev,nreg,istartcol,iendcol, &
298         &  config%i_cloud_pdf_shape == IPdfShapeGamma, &
299         &  cloud%fraction, cloud%fractional_std, region_fracs, &
300         &  od_scaling, config%cloud_fraction_threshold)
301
302    ! Compute wavelength-independent overlap matrices u_matrix and v_matrix
303    call calc_overlap_matrices(nlev, nreg, istartcol, iendcol, &
304         &  region_fracs, cloud%overlap_param, &
305         &  u_matrix, v_matrix, decorrelation_scaling=config%cloud_inhom_decorr_scaling, &
306         &  cloud_fraction_threshold=config%cloud_fraction_threshold, &
307         &  use_beta_overlap=config%use_beta_overlap, &
308         &  cloud_cover=flux%cloud_cover_lw)
309
310    if (config%iverbose >= 3) then
311      write(nulout,'(a)',advance='no') '  Processing columns'
312    end if
313
314    ! Main loop over columns
315    do jcol = istartcol, iendcol
316      ! --------------------------------------------------------
317      ! Section 2: Prepare column-specific variables and arrays
318      ! --------------------------------------------------------
319
320      if (config%iverbose >= 3) then
321        write(nulout,'(a)',advance='no') '.'
322      end if
323
324      ! Define which layers contain cloud; assume that
325      ! cloud%crop_cloud_fraction has already been called
326      is_clear_sky_layer = .true.
327      do jlev = 1,nlev
328        if (cloud%fraction(jcol,jlev) > 0.0_jprb) then
329          is_clear_sky_layer(jlev) = .false.
330        end if
331      end do
332
333      ! --------------------------------------------------------
334      ! Section 3: First loop over layers
335      ! --------------------------------------------------------
336      ! In this section the reflectance, transmittance and sources
337      ! are computed for each layer
338      do jlev = 1,nlev ! Start at top-of-atmosphere
339        ! --------------------------------------------------------
340        ! Section 3.1: Layer-specific initialization
341        ! --------------------------------------------------------
342        ! Array-wise assignments
343        gamma1 = 0.0_jprb
344        gamma2 = 0.0_jprb
345        Gamma_z1= 0.0_jprb
346        planck_top = 0.0_jprb
347        planck_diff = 0.0_jprb
348        transfer_rate(:,:) = 0.0_jprb
349        solution0 = 0.0_jprb
350        solution_diff = 0.0_jprb
351
352        ! Reset clear-sky absorption/scattering properties in each
353        ! region
354        od_region = 0.0_jprb
355        ssa_region = 0.0_jprb ! No Rayleigh scattering in longwave by default
356        g_region = 0.0_jprb
357
358        ! Set optical depth of clear-sky region (region 1) to the
359        ! gas/aerosol optical depth
360        od_region(1:ng,1) = od(1:ng, jlev, jcol)
361
362        ! Can't deal with zero gas optical depths, but this is now
363        ! dealt with in radiation_ifsrrtm
364        ! where (od_region(:,1) < 1.0e-15_jprb) od_region(:,1) = 1.0e-15_jprb
365
366        ! Set clear-sky single-scattering albedo and asymmetry
367        ! factor, if non-zero
368        if (config%do_lw_aerosol_scattering) then
369          ssa_region(1:ng,1) = ssa(1:ng, jlev, jcol)
370          g_region(1:ng,1)   = g(1:ng, jlev, jcol)
371        end if
372
373        ! --------------------------------------------------------
374        ! Section 3.2: Compute gamma variables
375        ! --------------------------------------------------------
376        if (is_clear_sky_layer(jlev)) then
377          ! --- Section 3.2a: Clear-sky case --------------------
378
379          nregActive = 1   ! Only region 1 (clear-sky) is active
380
381          ! Compute two-stream variables gamma1 and gamma2
382          call calc_two_stream_gammas_lw(ng, &
383               &  ssa_region(1:ng,1), g_region(1:ng,1), &
384               &  gamma1(1:ng,1), gamma2(1:ng,1))
385
386          if (config%use_expm_everywhere) then
387            ! Treat cloud-free layer as 3D: the matrix exponential
388            ! "expm" is used as much as possible (provided the
389            ! optical depth is not too high). ng3D is initially
390            ! set to the total number of g-points, then the
391            ! clear-sky optical depths are searched until a value
392            ! exceeds the threshold for treating 3D effects, and
393            ! if found it is reset to that.  Note that we are
394            ! assuming that the g-points have been reordered in
395            ! approximate order of gas optical depth.
396            ng3D = ng
397            do jg = 1, ng
398              if (od_region(jg,1) > config%max_gas_od_3D) then
399                ng3D = jg-1
400                exit
401              end if
402            end do
403          else
404            ! Otherwise treat cloud-free layer using the classical
405            ! Meador & Weaver formulae for all g-points
406            ng3D = 0
407          end if
408
409        else
410          ! --- Section 3.2b: Cloudy case -----------------------
411
412          ! Default number of g-points to treat with
413          ! matrix-exponential scheme
414          if (config%use_expm_everywhere) then
415            ng3D = ng   ! All g-points
416          else
417            ng3D = 0    ! No g-points
418          end if
419
420          ! Do we compute 3D effects; note that if we only have one
421          ! region and the sky is overcast then 3D calculations must
422          ! be turned off as there will be only one region
423          if (config%do_3d_effects .and. &
424               &  allocated(cloud%inv_cloud_effective_size) .and. &
425               &  .not. (nreg == 2 .and. cloud%fraction(jcol,jlev) &
426               &  > 1.0_jprb-config%cloud_fraction_threshold)) then
427            if (cloud%inv_cloud_effective_size(jcol,jlev) &
428                 &  > 0.0_jprb) then
429              ! 3D effects are only simulated if
430              ! inv_cloud_effective_size is defined and greater
431              ! than zero
432              ng3D = ng
433
434              ! The following is from the hydrostatic equation
435              ! and ideal gas law: dz = dp * R * T / (p * g)
436              dz = R_over_g &
437                   &  * (thermodynamics%pressure_hl(jcol,jlev+1) &
438                   &   - thermodynamics%pressure_hl(jcol,jlev)) &
439                   &  * (thermodynamics%temperature_hl(jcol,jlev) &
440                   &   + thermodynamics%temperature_hl(jcol,jlev+1)) &
441                   &  / (thermodynamics%pressure_hl(jcol,jlev) &
442                   &   + thermodynamics%pressure_hl(jcol,jlev+1))
443
444              ! Compute cloud edge length per unit area of gridbox
445              ! from rearranging Hogan & Shonk (2013) Eq. 45, but
446              ! adding a factor of (1-frac) so that a region that
447              ! fully occupies the gridbox (frac=1) has an edge of
448              ! zero. We can therefore use the fraction of clear sky,
449              ! region_fracs(1,jlev,jcol) for convenience instead. The
450              ! pi on the denominator means that this is actually edge
451              ! length with respect to a light ray with a random
452              ! azimuthal direction.
453              edge_length(1) = four_over_pi &
454                   &  * region_fracs(1,jlev,jcol)*(1.0_jprb-region_fracs(1,jlev,jcol)) &
455                   &  * min(cloud%inv_cloud_effective_size(jcol,jlev), &
456                   &        1.0_jprb / config%min_cloud_effective_size)
457              if (nreg > 2) then
458                ! The corresponding edge length between the two cloudy
459                ! regions is computed in the same way but treating the
460                ! optically denser of the two regions (region 3) as
461                ! the cloud; note that the fraction of this region may
462                ! be less than that of the optically less dense of the
463                ! two regions (region 2).  For increased flexibility,
464                ! the user may specify the effective size of
465                ! inhomogeneities separately from the cloud effective
466                ! size.
467                if (allocated(cloud%inv_inhom_effective_size)) then
468                  edge_length(2) = four_over_pi &
469                       &  * region_fracs(3,jlev,jcol)*(1.0_jprb-region_fracs(3,jlev,jcol)) &
470                       &  * min(cloud%inv_inhom_effective_size(jcol,jlev), &
471                       &        1.0_jprb / config%min_cloud_effective_size)
472                else
473                  edge_length(2) = four_over_pi &
474                       &  * region_fracs(3,jlev,jcol)*(1.0_jprb-region_fracs(3,jlev,jcol)) &
475                       &  * min(cloud%inv_cloud_effective_size(jcol,jlev), &
476                       &        1.0_jprb / config%min_cloud_effective_size)
477                end if
478                ! In the case of three regions, some of the cloud
479                ! boundary may go directly to the "thick" region
480                if (config%clear_to_thick_fraction > 0.0_jprb) then
481                  edge_length(3) = config%clear_to_thick_fraction &
482                       &  * min(edge_length(1), edge_length(2))
483                  edge_length(1) = edge_length(1) - edge_length(3)
484                  edge_length(2) = edge_length(2) - edge_length(3)
485                else
486                  edge_length(3) = 0.0_jprb
487                end if
488              end if
489
490              do jreg = 1, nreg-1
491                ! Compute lateral transfer rate from region jreg to
492                ! jreg+1 following Hogan & Shonk (2013) Eq. 47, but
493                ! multiplied by dz because the transfer rate is
494                ! vertically integrated across the depth of the layer
495                if (region_fracs(jreg,jlev,jcol) > epsilon(1.0_jprb)) then
496                  transfer_rate(jreg,jreg+1) = dz &
497                       &  * edge_length(jreg) &
498                       &  * tan_diffuse_angle_3d / region_fracs(jreg,jlev,jcol)
499                end if
500                ! Compute transfer rate from region jreg+1 to jreg
501                if (region_fracs(jreg+1,jlev,jcol) > epsilon(1.0_jprb)) then
502                  transfer_rate(jreg+1,jreg) = dz &
503                       &  * edge_length(jreg) &
504                       &  * tan_diffuse_angle_3d / region_fracs(jreg+1,jlev,jcol)
505                end if
506              end do
507
508              ! Compute transfer rates directly between regions 1 and
509              ! 3
510              if (edge_length(3) > 0.0_jprb) then
511                if (region_fracs(1,jlev,jcol) > epsilon(1.0_jprb)) then
512                  transfer_rate(1,3) = dz &
513                       &  * edge_length(3) &
514                       &  * tan_diffuse_angle_3d / region_fracs(1,jlev,jcol)
515                end if
516                if (region_fracs(3,jlev,jcol) > epsilon(1.0_jprb)) then
517                  transfer_rate(3,1) = dz &
518                       &  * edge_length(3) &
519                       &  * tan_diffuse_angle_3d / region_fracs(3,jlev,jcol)
520                end if
521              end if
522
523              ! Don't allow the transfer rate out of a region to be
524              ! equivalent to a loss of exp(-10) through the layer
525              where (transfer_rate > config%max_3d_transfer_rate)
526                transfer_rate = config%max_3d_transfer_rate
527              end where
528            end if ! Cloud has edge length required for 3D effects
529          end if ! Include 3D effects
530
531          ! In a cloudy layer the number of active regions equals
532          ! the number of regions
533          nregActive = nreg
534
535          ! Compute scattering properties of the regions at each
536          ! g-point, mapping from the cloud properties
537          ! defined in each band.
538          do jg = 1,ng
539            ! Mapping from g-point to band
540            iband = config%i_band_from_reordered_g_lw(jg)
541
542            ! Scattering optical depth of clear-sky region
543            scat_od = od_region(jg,1)*ssa_region(jg,1)
544
545            ! Loop over cloudy regions
546            do jreg = 2,nreg
547              ! Add scaled cloud optical depth to clear-sky value
548              od_region(jg,jreg) = od_region(jg,1) &
549                   &  + od_cloud(iband,jlev,jcol)*od_scaling(jreg,jlev,jcol)
550              if (config%do_lw_cloud_scattering) then
551                ! Compute single-scattering albedo of gas-cloud
552                ! combination
553                scat_od_cloud = od_cloud(iband,jlev,jcol) &
554                     &  * ssa_cloud(iband,jlev,jcol)*od_scaling(jreg,jlev,jcol)
555                ssa_region(jg,jreg) = (scat_od+scat_od_cloud) &
556                     &  / od_region(jg,jreg)
557                if (scat_od + scat_od_cloud > 0.0_jprb) then
558                  ! Compute asymmetry factor of gas-cloud
559                  ! combination
560                  g_region(jg,jreg) &
561                       &  = (scat_od*g_region(jg,1) &
562                       &  + scat_od_cloud*g_cloud(iband,jlev,jcol)) &
563                       &  / (scat_od + scat_od_cloud)
564                  ! else g_region is already set to zero
565                end if
566                ! else ssa_region is already set to zero
567              end if
568
569              ! Apply maximum cloud optical depth for stability in the
570              ! 3D case
571              if (od_region(jg,jreg) > config%max_cloud_od) then
572                od_region(jg,jreg) = config%max_cloud_od
573              end if
574
575            end do
576
577            ! Calculate two-stream variables gamma1 and gamma2 of
578            ! all regions at once
579            call calc_two_stream_gammas_lw(nreg, &
580                 &  ssa_region(jg,:), g_region(jg,:), &
581                 &  gamma1(jg,:), gamma2(jg,:))
582
583            ! Loop is in order of g-points with typically
584            ! increasing optical depth: if optical depth of
585            ! clear-sky region exceeds a threshold then turn off
586            ! 3D effects for any further g-points
587            if (ng3D == ng &
588                 &  .and. od_region(jg,1) > config%max_gas_od_3D) then
589              ng3D = jg-1
590            end if
591          end do ! Loop over g points
592        end if ! Cloudy level
593
594        ! --------------------------------------------------------
595        ! Section 3.3: Compute reflection, transmission and emission
596        ! --------------------------------------------------------
597        if (ng3D > 0) then
598          ! --- Section 3.3a: g-points with 3D effects ----------
599
600          ! 3D effects need to be represented in "ng3D" of the g
601          ! points.  This is done by creating ng3D square matrices
602          ! each of dimension 2*nreg by 2*nreg, computing the matrix
603          ! exponential, then computing the various
604          ! transmission/reflectance matrices from that.
605          do jreg = 1,nregActive
606            ! Write the diagonal elements of -Gamma1*z1
607            Gamma_z1(1:ng3D,jreg,jreg) &
608                 &  = od_region(1:ng3D,jreg)*gamma1(1:ng3D,jreg)
609            ! Write the diagonal elements of +Gamma2*z1
610            Gamma_z1(1:ng3D,jreg+nreg,jreg) &
611                 &  = od_region(1:ng3D,jreg)*gamma2(1:ng3D,jreg)
612
613            ! Write the vectors corresponding to the inhomogeneous
614            ! parts of the matrix ODE
615            planck_top(1:ng3D,nreg+jreg) = od_region(1:ng3D,jreg) &
616                 &  *(1.0_jprb-ssa_region(1:ng3D,jreg))*region_fracs(jreg,jlev,jcol) &
617                 &  *planck_hl(1:ng3D,jlev,jcol)*LwDiffusivity
618            planck_top(1:ng3D,jreg) = -planck_top(1:ng3D,nreg+jreg)
619            planck_diff(1:ng3D,nreg+jreg) = od_region(1:ng3D,jreg) &
620                 &  * (1.0_jprb-ssa_region(1:ng3D,jreg))*region_fracs(jreg,jlev,jcol) &
621                 &  * (planck_hl(1:ng3D,jlev+1,jcol) &
622                 &  -planck_hl(1:ng3D,jlev,jcol))*LwDiffusivity
623            planck_diff(1:ng3D,jreg) = -planck_diff(1:ng3D,nreg+jreg)
624          end do
625
626          if (nregActive < nreg) then
627            ! To avoid NaNs in some situations we need to put
628            ! non-zeros in the cloudy-sky regions even if the cloud
629            ! fraction is zero
630            do jreg = nregActive+1,nreg
631              Gamma_z1(1:ng3D,jreg,jreg) = Gamma_z1(1:ng3D,1,1)
632              Gamma_z1(1:ng3D,nreg+jreg,jreg) = Gamma_z1(1:ng3D,nreg+1,1)
633            end do
634          end if
635
636          ! Parameterization for the effective emissivity of the side
637          ! of the cloud
638          if (config%do_lw_side_emissivity &
639             & .and. region_fracs(1,jlev,jcol) > 0.0_jprb .and. region_fracs(2,jlev,jcol) > 0.0_jprb &
640             & .and. config%do_3d_effects &
641             & .and. cloud%inv_cloud_effective_size(jcol,jlev) > 0.0_jprb) then
642            aspect_ratio = 1.0_jprb / (min(cloud%inv_cloud_effective_size(jcol,jlev), &
643                 &                         1.0_jprb / config%min_cloud_effective_size) &
644                 &                     * region_fracs(1,jlev,jcol) * dz)
645            lateral_od(1:ng3D) = (aspect_ratio / (nreg-1.0_jprb)) &
646                 &  * sum(od_region(1:ng3D,2:nreg)*(1.0_jprb-ssa_region(1:ng3D,2:nreg)),2)
647            sqrt_1_minus_ssa(1:ng3D) = sqrt(1.0_jprb - ssa_region(1:ng3D,2))
648            side_emiss_thick(1:ng3D) = 2.0_jprb * sqrt_1_minus_ssa(1:ng3D) / &
649                 & (sqrt_1_minus_ssa(1:ng3D) &
650                 &  + sqrt(1.0_jprb-ssa_region(1:ng3D,2)*g_region(1:ng3D,2)))
651            side_emiss(1:ng3D) = (side_emiss_thin - side_emiss_thick(1:ng3D)) &
652                 &  / (lateral_od(1:ng3D) + 1.0_jprb) &
653                 &  + side_emiss_thick(1:ng3D)
654          else
655            side_emiss(1:ng3D) = 1.0_jprb
656          end if
657
658          do jreg = 1,nregActive-1
659            ! Add the terms assocated with 3D transport to Gamma1*z1.
660            ! First the rate from region jreg to jreg+1
661            Gamma_z1(1:ng3D,jreg,jreg) = Gamma_z1(1:ng3D,jreg,jreg) &
662                 &  + transfer_rate(jreg,jreg+1)
663            Gamma_z1(1:ng3D,jreg+1,jreg) = -transfer_rate(jreg,jreg+1)
664
665            ! Next the rate from region jreg+1 to jreg
666            if (jreg > 1) then
667              ! Flow between one cloudy region and another
668              Gamma_z1(1:ng3D,jreg+1,jreg+1) = Gamma_z1(1:ng3D,jreg+1,jreg+1) &
669                   &  + transfer_rate(jreg+1,jreg)
670              Gamma_z1(1:ng3D,jreg,jreg+1) = -transfer_rate(jreg+1,jreg)
671            else
672              ! Only for the lateral transfer between cloud and clear
673              ! skies do we account for the effective emissivity of
674              ! the side of the cloud
675              Gamma_z1(1:ng3D,jreg+1,jreg+1) = Gamma_z1(1:ng3D,jreg+1,jreg+1) &
676                   &  + side_emiss(1:ng3D) * transfer_rate(jreg+1,jreg)
677              Gamma_z1(1:ng3D,jreg,jreg+1) = -side_emiss(1:ng3D)*transfer_rate(jreg+1,jreg)
678            end if
679          end do
680
681          ! Possible flow between regions a and c
682          if (edge_length(3) > 0.0_jprb) then
683            Gamma_z1(1:ng3D,1,1) = Gamma_z1(1:ng3D,1,1) &
684                 &  + transfer_rate(1,3)
685            Gamma_z1(1:ng3D,3,1) = -transfer_rate(1,3)
686            Gamma_z1(1:ng3D,3,3) = Gamma_z1(1:ng3D,3,3) &
687                 &  + side_emiss(1:ng3D) * transfer_rate(3,1)
688            Gamma_z1(1:ng3D,1,3) = -side_emiss(1:ng3D)*transfer_rate(3,1)
689          end if
690
691          ! Copy Gamma1*z1
692          Gamma_z1(1:ng3D,nreg+1:2*nreg,nreg+1:2*nreg) &
693               &  = -Gamma_z1(1:ng3D,1:nreg,1:nreg)
694          ! Copy Gamma2*z1
695          Gamma_z1(1:ng3D,1:nreg,nreg+1:2*nreg) &
696               &  = -Gamma_z1(1:ng3D,nreg+1:2*nreg,1:nreg)
697
698          ! Compute the parts of the particular solution
699          solution_diff(1:ng3D,1:2*nreg) &
700               &  = solve_vec(ng,ng3D,2*nreg,Gamma_z1,planck_diff)
701          solution_diff(1:ng3D,1:2*nreg) = - solution_diff(1:ng3D,1:2*nreg)
702          solution0(1:ng3D,1:2*nreg) = solve_vec(ng,ng3D,2*nreg,Gamma_z1, &
703               &  solution_diff-planck_top)
704
705          ! Compute the matrix exponential of Gamma_z1, returning the
706          ! result in-place
707          call expm(ng, ng3D, 2*nreg, Gamma_z1, IMatrixPatternDense)
708
709          ! Update count of expm calls
710          n_calls_expm = n_calls_expm + ng3D
711
712          ! Diffuse reflectance matrix
713          reflectance(1:ng3D,:,:,jlev) = -solve_mat(ng,ng3D,nreg, &
714               &  Gamma_z1(1:ng3D,1:nreg,1:nreg), &
715               &  Gamma_z1(1:ng3D,1:nreg,nreg+1:2*nreg))
716          ! Diffuse transmission matrix
717          transmittance(1:ng3D,:,:,jlev) = mat_x_mat(ng,ng3D,nreg, &
718               &  Gamma_z1(1:ng3D,nreg+1:2*nreg,1:nreg), &
719               &  reflectance(1:ng3D,:,:,jlev)) &
720               &  + Gamma_z1(1:ng3D,nreg+1:2*nreg,nreg+1:2*nreg)
721
722          ! Upward and downward sources due to emission within the
723          ! layer
724          tmp_vectors(1:ng3D,:) = solution0(1:ng3D,1:nreg) &
725               &  + solution_diff(1:ng3D,1:nreg) &
726               &  - mat_x_vec(ng,ng3D,nreg,Gamma_z1(:,1:nreg,nreg+1:2*nreg), &
727               &  solution0(:,nreg+1:2*nreg))
728          source_up(1:ng3D,:,jlev) = solution0(1:ng3D,1:nreg) &
729               &  - solve_vec(ng,ng3D,nreg,Gamma_z1(:,1:nreg,1:nreg), &
730               &  tmp_vectors)
731          tmp_vectors(1:ng3D,:) = source_up(1:ng3D,:,jlev) &
732               &  - solution0(1:ng3D,1:nreg)
733          source_dn(1:ng3D,:,jlev) = mat_x_vec(ng,ng3D,nreg,&
734               &  Gamma_z1(:,nreg+1:2*nreg,1:nreg),tmp_vectors) &
735               &  + solution0(1:ng3D,nreg+1:2*nreg) &
736               &  - mat_x_vec(ng,ng3D,nreg, &
737               &  Gamma_z1(:,nreg+1:2*nreg,nreg+1:2*nreg), &
738               &  solution0(:,nreg+1:2*nreg)) &
739               &  + solution_diff(1:ng3D,nreg+1:2*nreg)
740        end if ! we are treating 3D effects for some g points
741
742        ! --- Section 3.3b: g-points without 3D effects ----------
743
744        ! Compute reflectance, transmittance and upward and downward
745        ! sources for clear skies, using the Meador-Weaver formulas
746        ! for reflectance and transmittance, and equivalent solutions
747        ! to the coupled ODEs for the sources.
748        call calc_reflectance_transmittance_lw(ng, &
749             &  od_region(1:ng,1), &
750             &  gamma1(1:ng,1), gamma2(1:ng,1), &
751             &  planck_hl(1:ng,jlev,jcol), planck_hl(1:ng,jlev+1,jcol), &
752             &  ref_clear(1:ng,jlev), trans_clear(1:ng,jlev), &
753             &  source_up_clear(1:ng,jlev), source_dn_clear(1:ng,jlev))
754
755        n_calls_meador_weaver = n_calls_meador_weaver + ng
756
757        if (ng3D < ng) then
758          ! Some of the g points are to be treated using the
759          ! conventional plane-parallel method.  First zero the
760          ! relevant parts of the arrays
761          reflectance  (ng3D+1:ng,:,:,jlev) = 0.0_jprb
762          transmittance(ng3D+1:ng,:,:,jlev) = 0.0_jprb
763          source_up    (ng3D+1:ng,:,jlev) = 0.0_jprb
764          source_dn    (ng3D+1:ng,:,jlev) = 0.0_jprb
765
766          ! Since there is no lateral transport, the clear-sky parts
767          ! of the arrays can be copied from the clear-sky arrays
768          reflectance(ng3D+1:ng,1,1,jlev) = ref_clear(ng3D+1:ng,jlev)
769          transmittance(ng3D+1:ng,1,1,jlev) = trans_clear(ng3D+1:ng,jlev)
770          source_up(ng3D+1:ng,1,jlev) = region_fracs(1,jlev,jcol)*source_up_clear(ng3D+1:ng,jlev)
771          source_dn(ng3D+1:ng,1,jlev) = region_fracs(1,jlev,jcol)*source_dn_clear(ng3D+1:ng,jlev)
772
773          ! Compute reflectance, transmittance and upward and downward
774          ! sources, for each cloudy region, using the Meador-Weaver
775          ! formulas for reflectance and transmittance, and equivalent
776          ! solutions to the coupled ODEs for the sources.
777          do jreg = 2, nregActive
778            call calc_reflectance_transmittance_lw(ng-ng3D, &
779                 &  od_region(ng3D+1:ng,jreg), &
780                 &  gamma1(ng3D+1:ng,jreg), gamma2(ng3D+1:ng,jreg), &
781                 &  region_fracs(jreg,jlev,jcol)*planck_hl(ng3D+1:ng,jlev,jcol), &
782                 &  region_fracs(jreg,jlev,jcol)*planck_hl(ng3D+1:ng,jlev+1,jcol), &
783                 &  reflectance(ng3D+1:ng,jreg,jreg,jlev), &
784                 &  transmittance(ng3D+1:ng,jreg,jreg,jlev), &
785                 &  source_up(ng3D+1:ng,jreg,jlev), &
786                 &  source_dn(ng3D+1:ng,jreg,jlev))
787          end do
788          n_calls_meador_weaver &
789               &  = n_calls_meador_weaver + (ng-ng3D)*(nregActive-1)
790        end if
791
792      end do ! Loop over levels
793
794      ! --------------------------------------------------------
795      ! Section 4: Compute total sources and albedos
796      ! --------------------------------------------------------
797
798      total_albedo(:,:,:,:) = 0.0_jprb
799      total_source(:,:,:) = 0.0_jprb
800
801      if (config%do_clear) then
802        total_albedo_clear(:,:) = 0.0_jprb
803        total_source_clear(:,:) = 0.0_jprb
804      end if
805
806      ! Calculate the upwelling radiation emitted by the surface, and
807      ! copy the surface albedo into total_albedo
808      do jreg = 1,nreg
809        do jg = 1,ng
810          ! region_fracs(jreg,nlev,jcol) is the fraction of each
811          ! region in the lowest model level
812          total_source(jg,jreg,nlev+1) = region_fracs(jreg,nlev,jcol)*emission(jg,jcol)
813          total_albedo(jg,jreg,jreg,nlev+1) = albedo(jg,jcol)
814        end do
815      end do
816      ! Equivalent surface values for computing clear-sky fluxes
817      if (config%do_clear) then
818        do jg = 1,ng
819          total_source_clear(jg,nlev+1) = emission(jg,jcol)
820        end do
821        ! In the case of surface albedo there is no dependence on
822        ! cloud fraction so we can copy the all-sky value
823        total_albedo_clear(1:ng,nlev+1) = total_albedo(1:ng,1,1,nlev+1)
824      end if
825
826      ! Loop back up through the atmosphere computing the total albedo
827      ! and the total upwelling due to emission below each level
828      do jlev = nlev,1,-1
829        if (config%do_clear) then
830          ! Use adding method for clear-sky arrays; note that there
831          ! is no need to consider "above" and "below" quantities
832          ! since with no cloud overlap to worry about, these are
833          ! the same
834          inv_denom_scalar(:) = 1.0_jprb &
835               &  / (1.0_jprb - total_albedo_clear(:,jlev+1)*ref_clear(:,jlev))
836          total_albedo_clear(:,jlev) = ref_clear(:,jlev) &
837               &  + trans_clear(:,jlev)*trans_clear(:,jlev)*total_albedo_clear(:,jlev+1) &
838               &  * inv_denom_scalar(:)
839          total_source_clear(:,jlev) = source_up_clear(:,jlev) &
840               &  + trans_clear(:,jlev)*(total_source_clear(:,jlev+1) &
841               &  + total_albedo_clear(:,jlev+1)*source_dn_clear(:,jlev)) &
842               &  * inv_denom_scalar(:)
843        end if
844
845        if (is_clear_sky_layer(jlev)) then
846          ! Clear-sky layer: use scalar adding method
847          inv_denom_scalar(:) = 1.0_jprb &
848               &  / (1.0_jprb - total_albedo(:,1,1,jlev+1)*reflectance(:,1,1,jlev))
849          total_albedo_below = 0.0_jprb
850          total_albedo_below(:,1,1) = reflectance(:,1,1,jlev) &
851               &  + transmittance(:,1,1,jlev)  * transmittance(:,1,1,jlev) &
852               &  * total_albedo(:,1,1,jlev+1) * inv_denom_scalar(:)
853          total_source_below = 0.0_jprb
854          total_source_below(:,1) = source_up(:,1,jlev) &
855               &  + transmittance(:,1,1,jlev)*(total_source(:,1,jlev+1) &
856               &  + total_albedo(:,1,1,jlev+1)*source_dn(:,1,jlev)) &
857               &  * inv_denom_scalar(:)
858        else if (config%do_3d_effects .or. &
859             &   config%do_3d_lw_multilayer_effects) then
860          ! Cloudy layer: use matrix adding method
861          denominator = identity_minus_mat_x_mat(ng,ng,nreg,&
862               &  total_albedo(:,:,:,jlev+1), reflectance(:,:,:,jlev))
863          total_albedo_below = reflectance(:,:,:,jlev) &
864               &  + mat_x_mat(ng,ng,nreg,transmittance(:,:,:,jlev), &
865               &  solve_mat(ng,ng,nreg,denominator, &
866               &  mat_x_mat(ng,ng,nreg,total_albedo(:,:,:,jlev+1), &
867               &  transmittance(:,:,:,jlev))))
868          total_source_below = source_up(:,:,jlev) &
869               &  + mat_x_vec(ng,ng,nreg,transmittance(:,:,:,jlev), &
870               &  solve_vec(ng,ng,nreg,denominator, &
871               &  total_source(:,:,jlev+1) + mat_x_vec(ng,ng,nreg, &
872               &  total_albedo(:,:,:,jlev+1),source_dn(:,:,jlev))))
873        else
874          ! Cloudy layer for which reflectance, transmittance and
875          ! total_albedo matrices are diagonal
876          total_albedo_below = 0.0_jprb
877          total_source_below = 0.0_jprb
878          do jreg = 1,nreg
879            inv_denom_scalar(:) = 1.0_jprb / (1.0_jprb &
880                 &  - total_albedo(:,jreg,jreg,jlev+1)*reflectance(:,jreg,jreg,jlev))
881            total_albedo_below(:,jreg,jreg) = reflectance(:,jreg,jreg,jlev) &
882                 &  + transmittance(:,jreg,jreg,jlev)*transmittance(:,jreg,jreg,jlev) &
883                 &  * total_albedo(:,jreg,jreg,jlev+1) &
884                 &  * inv_denom_scalar(:)
885            total_source_below(:,jreg) = source_up(:,jreg,jlev) &
886                 &  + transmittance(:,jreg,jreg,jlev)*(total_source(:,jreg,jlev+1) &
887                 &  + total_albedo(:,jreg,jreg,jlev+1)*source_dn(:,jreg,jlev)) &
888                 &  * inv_denom_scalar(:)
889          end do
890
891        end if
892
893        ! Account for cloud overlap when converting albedo and
894        ! source below a layer interface to the equivalent values
895        ! just above
896        if (is_clear_sky_layer(jlev) .and. is_clear_sky_layer(jlev-1)) then
897          ! If both layers are cloud free, this is trivial...
898          total_albedo(:,:,:,jlev) = 0.0_jprb
899          total_albedo(:,1,1,jlev) = total_albedo_below(:,1,1)
900          total_source(:,:,jlev) = 0.0_jprb
901          total_source(:,1,jlev) = total_source_below(:,1)
902        else
903          total_source(:,:,jlev) = singlemat_x_vec(ng,ng,nreg,&
904               &  u_matrix(:,:,jlev,jcol), total_source_below)
905
906!          if (config%do_3d_effects .or. config%do_3d_lw_multilayer_effects) then
907          if (config%do_3d_lw_multilayer_effects) then
908            ! Use the overlap matrices u_matrix and v_matrix
909            total_albedo(:,:,:,jlev) = singlemat_x_mat(ng,ng,nreg,&
910                 &  u_matrix(:,:,jlev,jcol), &
911                 &  mat_x_singlemat(ng,ng,nreg,total_albedo_below,&
912                 &  v_matrix(:,:,jlev,jcol)))
913          else
914            total_albedo(:,:,:,jlev) = 0.0_jprb
915            ! "total_albedo" is diagonal and we wish to exclude
916            ! anomalous horizontal transport described by Shonk &
917            ! Hogan (2008).  Therefore, the operation we perform is
918            ! essentially diag(total_albedo) = matmul(transpose(v_matrix),
919            ! diag(total_albedo_below)).
920            do jreg = 1,nreg
921              do jreg2 = 1,nreg
922                total_albedo(:,jreg,jreg,jlev) &
923                     &  = total_albedo(:,jreg,jreg,jlev) &
924                     &  + total_albedo_below(:,jreg2,jreg2) &
925                     &  * v_matrix(jreg2,jreg,jlev,jcol)
926              end do
927            end do
928          end if
929        end if
930
931      end do ! Reverse loop over levels
932
933      ! --------------------------------------------------------
934      ! Section 5: Compute fluxes
935      ! --------------------------------------------------------
936
937      ! Top-of-atmosphere fluxes into the regions of the top-most
938      ! layer: zero since we assume no extraterrestrial longwave
939      ! radiation (the shortwave scheme accounts for solar radiation
940      ! in the "longwave" part of the spectrum)
941      flux_dn_below = 0.0_jprb
942      flux%lw_dn(jcol,1) = 0.0_jprb
943      if (config%do_clear) then
944        flux_dn_clear = 0.0_jprb
945        flux%lw_dn_clear(jcol,1) = 0.0_jprb
946      end if
947
948      ! Store the outgoing longwave radiation at top-of-atmosphere
949      flux%lw_up(jcol,1) = sum(sum(total_source(:,:,1),1))
950      if (config%do_clear) then
951        flux%lw_up_clear(jcol,1) = sum(total_source_clear(:,1))
952      end if
953
954      if (config%do_save_spectral_flux) then
955        call indexed_sum(sum(total_source(:,:,1),2), &
956             &           config%i_spec_from_reordered_g_lw, &
957             &           flux%lw_up_band(:,jcol,1))
958        flux%lw_dn_band(:,jcol,1) = 0.0_jprb
959        if (config%do_clear) then
960          call indexed_sum(total_source_clear(:,1), &
961               &           config%i_spec_from_reordered_g_lw, &
962               &           flux%lw_up_clear_band(:,jcol,1))
963          flux%lw_dn_clear_band(:,jcol,1) = 0.0_jprb
964        end if
965      end if
966
967      ! Final loop back down through the atmosphere to compute fluxes
968      do jlev = 1,nlev
969        if (config%do_clear) then
970          ! Scalar operations for clear-sky fluxes
971          flux_dn_clear(:) = (trans_clear(:,jlev)*flux_dn_clear(:) &
972               + ref_clear(:,jlev)*total_source_clear(:,jlev+1) &
973               + source_dn_clear(:,jlev)) &
974               / (1.0_jprb - ref_clear(:,jlev)*total_albedo_clear(:,jlev+1))
975          flux_up_clear(:) = total_source_clear(:,jlev+1) &
976               + total_albedo_clear(:,jlev+1)*flux_dn_clear(:)
977        end if
978
979        if (is_clear_sky_layer(jlev)) then
980          ! Scalar operations for clear-sky layer
981          flux_dn_above(:,1) = (transmittance(:,1,1,jlev)*flux_dn_below(:,1) &
982               &  + reflectance(:,1,1,jlev)*total_source(:,1,jlev+1) &
983               &  + source_dn(:,1,jlev)) &
984               &  / (1.0_jprb - reflectance(:,1,1,jlev)*total_albedo(:,1,1,jlev+1))
985          flux_dn_above(:,2:nreg) = 0.0_jprb
986          flux_up_above(:,1) = total_source(:,1,jlev+1) &
987               &  + total_albedo(:,1,1,jlev+1)*flux_dn_above(:,1)
988          flux_up_above(:,2:nreg) = 0.0_jprb
989        else if (config%do_3d_effects .or. config%do_3d_lw_multilayer_effects) then
990          ! Matrix operations for cloudy layer
991          denominator = identity_minus_mat_x_mat(ng,ng,nreg,reflectance(:,:,:,jlev), &
992               &  total_albedo(:,:,:,jlev+1))
993          flux_dn_above = solve_vec(ng,ng,nreg,denominator, &
994               &  mat_x_vec(ng,ng,nreg,transmittance(:,:,:,jlev),flux_dn_below) &
995               &  + mat_x_vec(ng,ng,nreg,reflectance(:,:,:,jlev), &
996               &  total_source(:,:,jlev+1)) &
997               &  + source_dn(:,:,jlev))
998          flux_up_above = mat_x_vec(ng,ng,nreg,total_albedo(:,:,:,jlev+1), &
999               &  flux_dn_above) + total_source(:,:,jlev+1)
1000        else
1001          do jreg = 1,nreg
1002            ! Scalar operations for all regions, requiring that
1003            ! reflectance, transmittance and total_albedo are diagonal
1004            flux_dn_above(:,jreg) = (transmittance(:,jreg,jreg,jlev)*flux_dn_below(:,jreg) &
1005                 &  + reflectance(:,jreg,jreg,jlev)*total_source(:,jreg,jlev+1) &
1006                 &  + source_dn(:,jreg,jlev)) &
1007                 &  / (1.0_jprb - reflectance(:,jreg,jreg,jlev) &
1008                 &              * total_albedo(:,jreg,jreg,jlev+1))
1009            flux_up_above(:,jreg) = total_source(:,jreg,jlev+1) &
1010                 &  + total_albedo(:,jreg,jreg,jlev+1)*flux_dn_above(:,jreg)
1011          end do
1012        end if
1013
1014        ! Account for overlap rules in translating fluxes just above
1015        ! a layer interface to the values just below
1016        if (is_clear_sky_layer(jlev) .and. is_clear_sky_layer(jlev+1)) then
1017          flux_dn_below = flux_dn_above
1018        else
1019          flux_dn_below = singlemat_x_vec(ng,ng,nreg,v_matrix(:,:,jlev+1,jcol), &
1020               &    flux_dn_above)
1021        end if
1022
1023        ! Store the broadband fluxes
1024        flux%lw_up(jcol,jlev+1) = sum(sum(flux_up_above,1))
1025        flux%lw_dn(jcol,jlev+1) = sum(sum(flux_dn_above,1))
1026        if (config%do_clear) then
1027          flux%lw_up_clear(jcol,jlev+1) = sum(flux_up_clear)
1028          flux%lw_dn_clear(jcol,jlev+1) = sum(flux_dn_clear)
1029        end if
1030
1031        if (config%do_save_spectral_flux) then
1032          call indexed_sum(sum(flux_up_above,2), &
1033               &           config%i_spec_from_reordered_g_lw, &
1034               &           flux%lw_up_band(:,jcol,jlev+1))
1035          call indexed_sum(sum(flux_dn_above,2), &
1036               &           config%i_spec_from_reordered_g_lw, &
1037               &           flux%lw_dn_band(:,jcol,jlev+1))
1038          if (config%do_clear) then
1039            call indexed_sum(flux_up_clear, &
1040                 &           config%i_spec_from_reordered_g_lw, &
1041                 &           flux%lw_up_clear_band(:,jcol,jlev+1))
1042            call indexed_sum(flux_dn_clear, &
1043                 &           config%i_spec_from_reordered_g_lw, &
1044                 &           flux%lw_dn_clear_band(:,jcol,jlev+1))
1045          end if
1046        end if
1047
1048      end do ! Final loop over levels
1049
1050      ! Store surface spectral downwelling fluxes, which at this point
1051      ! are at the surface
1052      flux%lw_dn_surf_g(:,jcol) = sum(flux_dn_above,2)
1053      if (config%do_clear) then
1054        flux%lw_dn_surf_clear_g(:,jcol) = flux_dn_clear
1055      end if
1056
1057      ! Compute the longwave derivatives needed by Hogan and Bozzo
1058      ! (2015) approximate radiation update scheme
1059      if (config%do_lw_derivatives) then
1060        ! Note that at this point flux_up_above contains the spectral
1061        ! fluxes into the regions of the lowest layer; we sum over
1062        ! regions first to provide a simple spectral flux upwelling
1063        ! from the surface
1064        call calc_lw_derivatives_matrix(ng, nlev, nreg, jcol, transmittance, &
1065             &  u_matrix(:,:,:,jcol), sum(flux_up_above,2), flux%lw_derivatives)
1066      end if
1067       
1068    end do ! Loop over columns
1069
1070    if (config%iverbose >= 3) then
1071      write(nulout,*)
1072    end if
1073
1074    ! Report number of calls to each method of solving single-layer
1075    ! two-stream equations
1076    if (config%iverbose >= 4) then
1077      write(nulout,'(a,i0)') '  Matrix-exponential calls: ', n_calls_expm
1078      write(nulout,'(a,i0)') '  Meador-Weaver calls: ', n_calls_meador_weaver
1079    end if
1080
1081    if (lhook) call dr_hook('radiation_spartacus_lw:solver_spartacus_lw',1,hook_handle)
1082
1083  end subroutine solver_spartacus_lw
1084
1085end module radiation_spartacus_lw
Note: See TracBrowser for help on using the repository browser.