source: LMDZ6/trunk/libf/phylmd/ecrad/radiation_tripleclouds_sw.F90 @ 4629

Last change on this file since 4629 was 4489, checked in by lguez, 19 months ago

Merge LMDZ_ECRad branch back into trunk!

File size: 26.3 KB
Line 
1! radiation_tripleclouds_sw.F90 - Shortwave "Tripleclouds" solver
2!
3! (C) Copyright 2016- ECMWF.
4!
5! This software is licensed under the terms of the Apache Licence Version 2.0
6! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
7!
8! In applying this licence, ECMWF does not waive the privileges and immunities
9! granted to it by virtue of its status as an intergovernmental organisation
10! nor does it submit to any jurisdiction.
11!
12! Author:  Robin Hogan
13! Email:   r.j.hogan@ecmwf.int
14!
15! 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-10-23  R. Hogan  Renamed single-character variables
19!   2018-10-08  R. Hogan  Call calc_region_properties
20!   2019-01-02  R. Hogan  Fixed problem of do_save_spectral_flux .and. .not. do_sw_direct
21!   2020-09-18  R. Hogan  Replaced some array expressions with loops for speed
22
23module radiation_tripleclouds_sw
24
25  public
26
27contains
28  ! Provides elemental function "delta_eddington"
29#include "radiation_delta_eddington.h"
30
31  ! Small routine for scaling cloud optical depth in the cloudy
32  ! regions
33#include "radiation_optical_depth_scaling.h"
34
35  !---------------------------------------------------------------------
36  ! This module contains just one subroutine, the shortwave
37  ! "Tripleclouds" solver in which cloud inhomogeneity is treated by
38  ! dividing each model level into three regions, one clear and two
39  ! cloudy (with differing optical depth). This approach was described
40  ! by Shonk and Hogan (2008).
41  subroutine solver_tripleclouds_sw(nlev,istartcol,iendcol, &
42       &  config, single_level, cloud, &
43       &  od, ssa, g, od_cloud, ssa_cloud, g_cloud, &
44       &  albedo_direct, albedo_diffuse, incoming_sw, &
45       &  flux)
46
47    use parkind1, only           : jprb
48    use yomhook,  only           : lhook, dr_hook
49
50!    use radiation_io, only             : nulout
51    use radiation_config, only         : config_type, IPdfShapeGamma
52    use radiation_single_level, only   : single_level_type
53    use radiation_cloud, only          : cloud_type
54    use radiation_regions, only        : calc_region_properties
55    use radiation_overlap, only        : calc_overlap_matrices
56    use radiation_flux, only           : flux_type, &
57         &                               indexed_sum, add_indexed_sum
58    use radiation_matrix, only         : singlemat_x_vec
59    use radiation_two_stream, only     : calc_two_stream_gammas_sw, &
60         &                       calc_reflectance_transmittance_sw
61
62    implicit none
63
64    ! Inputs
65    integer, intent(in) :: nlev               ! number of model levels
66    integer, intent(in) :: istartcol, iendcol ! range of columns to process
67    type(config_type),        intent(in) :: config
68    type(single_level_type),  intent(in) :: single_level
69    type(cloud_type),         intent(in) :: cloud
70
71    ! Gas and aerosol optical depth, single-scattering albedo and
72    ! asymmetry factor at each shortwave g-point
73!    real(jprb), intent(in), dimension(istartcol:iendcol,nlev,config%n_g_sw) :: &
74    real(jprb), intent(in), dimension(config%n_g_sw,nlev,istartcol:iendcol) :: &
75         &  od, ssa, g
76
77    ! Cloud and precipitation optical depth, single-scattering albedo and
78    ! asymmetry factor in each shortwave band
79    real(jprb), intent(in), dimension(config%n_bands_sw,nlev,istartcol:iendcol) :: &
80         &  od_cloud, ssa_cloud, g_cloud
81
82    ! Optical depth, single scattering albedo and asymmetry factor in
83    ! each g-point including gas, aerosol and clouds
84    real(jprb), dimension(config%n_g_sw) :: od_total, ssa_total, g_total
85
86    ! Direct and diffuse surface albedos, and the incoming shortwave
87    ! flux into a plane perpendicular to the incoming radiation at
88    ! top-of-atmosphere in each of the shortwave g points
89    real(jprb), intent(in), dimension(config%n_g_sw,istartcol:iendcol) :: &
90         &  albedo_direct, albedo_diffuse, incoming_sw
91
92    ! Output
93    type(flux_type), intent(inout):: flux
94
95    ! Local constants
96    integer, parameter :: nregions = 3
97
98    ! In a clear-sky layer this will be 1, otherwise equal to nregions
99    integer :: nreg
100
101    ! Local variables
102   
103    ! The area fractions of each region
104    real(jprb) :: region_fracs(1:nregions,nlev,istartcol:iendcol)
105
106    ! The scaling used for the optical depth in the cloudy regions
107    real(jprb) :: od_scaling(2:nregions,nlev,istartcol:iendcol)
108
109    ! Directional overlap matrices defined at all layer interfaces
110    ! including top-of-atmosphere and the surface
111    real(jprb), dimension(nregions,nregions,nlev+1, &
112         &                istartcol:iendcol) :: u_matrix, v_matrix
113
114    ! Two-stream variables
115    real(jprb), dimension(config%n_g_sw) :: gamma1, gamma2, gamma3
116
117    ! Diffuse reflection and transmission matrices of each layer
118    real(jprb), dimension(config%n_g_sw, nregions, nlev) &
119         &  :: reflectance, transmittance
120
121    ! Terms translating the direct flux entering the layer from above
122    ! to the reflected radiation exiting upwards (ref_dir) and the
123    ! scattered radiation exiting downwards (trans_dir_diff), along with the
124    ! direct unscattered transmission matrix (trans_dir_dir).
125    real(jprb), dimension(config%n_g_sw, nregions, nlev) &
126         &  :: ref_dir, trans_dir_diff, trans_dir_dir
127
128    ! Total albedo of the atmosphere/surface just above a layer
129    ! interface with respect to downwelling diffuse and direct
130    ! (respectively) radiation at that interface, where level index =
131    ! 1 corresponds to the top-of-atmosphere
132    real(jprb), dimension(config%n_g_sw, nregions, nlev+1) &
133         &  :: total_albedo, total_albedo_direct
134
135    ! ...equivalent values for clear-skies
136    real(jprb), dimension(config%n_g_sw, nlev+1) &
137         &  :: total_albedo_clear, total_albedo_clear_direct
138
139    ! Total albedo of the atmosphere just below a layer interface
140    real(jprb), dimension(config%n_g_sw, nregions) &
141         &  :: total_albedo_below, total_albedo_below_direct
142
143    ! Direct downwelling flux below and above an interface between
144    ! layers into a plane perpendicular to the direction of the sun
145    real(jprb), dimension(config%n_g_sw, nregions) &
146         &  :: direct_dn, direct_dn_below
147    ! Diffuse equivalents
148    real(jprb), dimension(config%n_g_sw, nregions) &
149         &  :: flux_dn, flux_dn_below, flux_up
150
151    ! ...clear-sky equivalent (no distinction between "above/below")
152    real(jprb), dimension(config%n_g_sw) &
153         &  :: direct_dn_clear, flux_dn_clear, flux_up_clear
154
155    ! Clear-sky equivalent, but actually its reciprocal to replace
156    ! some divisions by multiplications
157    real(jprb), dimension(config%n_g_sw, nregions) :: inv_denom
158
159    ! Identify clear-sky layers, with pseudo layers for outer space
160    ! and below the ground, both treated as single-region clear skies
161    logical :: is_clear_sky_layer(0:nlev+1)
162
163    ! Scattering optical depth of gas+aerosol and of cloud
164    real(jprb) :: scat_od, scat_od_cloud
165
166    real(jprb) :: mu0
167
168    integer :: jcol, jlev, jg, jreg, iband, jreg2, ng
169
170    real(jprb) :: hook_handle
171
172    if (lhook) call dr_hook('radiation_tripleclouds_sw:solver_tripleclouds_sw',0,hook_handle)
173
174    ! --------------------------------------------------------
175    ! Section 1: Prepare general variables and arrays
176    ! --------------------------------------------------------
177    ! Copy array dimensions to local variables for convenience
178    ng   = config%n_g_sw
179
180    ! Compute the wavelength-independent region fractions and
181    ! optical-depth scalings
182    call calc_region_properties(nlev,nregions,istartcol,iendcol, &
183         &  config%i_cloud_pdf_shape == IPdfShapeGamma, &
184         &  cloud%fraction, cloud%fractional_std, region_fracs, &
185         &  od_scaling, config%cloud_fraction_threshold)
186
187    ! Compute wavelength-independent overlap matrices u_matrix and
188    ! v_matrix
189    call calc_overlap_matrices(nlev,nregions,istartcol,iendcol, &
190         &  region_fracs, cloud%overlap_param, &
191         &  u_matrix, v_matrix, &
192         &  decorrelation_scaling=config%cloud_inhom_decorr_scaling, &
193         &  cloud_fraction_threshold=config%cloud_fraction_threshold, &
194         &  use_beta_overlap=config%use_beta_overlap, &
195         &  cloud_cover=flux%cloud_cover_sw)
196
197    ! Main loop over columns
198    do jcol = istartcol, iendcol
199      ! --------------------------------------------------------
200      ! Section 2: Prepare column-specific variables and arrays
201      ! --------------------------------------------------------
202
203      ! Copy local cosine of the solar zenith angle
204      mu0 = single_level%cos_sza(jcol)
205
206      ! Skip profile if sun is too low in the sky
207      if (mu0 < 1.0e-10_jprb) then
208        flux%sw_dn(jcol,:) = 0.0_jprb
209        flux%sw_up(jcol,:) = 0.0_jprb
210        if (allocated(flux%sw_dn_direct)) then
211          flux%sw_dn_direct(jcol,:) = 0.0_jprb
212        end if
213        if (config%do_clear) then
214          flux%sw_dn_clear(jcol,:) = 0.0_jprb
215          flux%sw_up_clear(jcol,:) = 0.0_jprb
216          if (allocated(flux%sw_dn_direct_clear)) then
217            flux%sw_dn_direct_clear(jcol,:) = 0.0_jprb
218          end if
219        end if
220
221        if (config%do_save_spectral_flux) then
222          flux%sw_dn_band(:,jcol,:) = 0.0_jprb
223          flux%sw_up_band(:,jcol,:) = 0.0_jprb
224          if (allocated(flux%sw_dn_direct_band)) then
225            flux%sw_dn_direct_band(:,jcol,:) = 0.0_jprb
226          end if
227          if (config%do_clear) then
228            flux%sw_dn_clear_band(:,jcol,:) = 0.0_jprb
229            flux%sw_up_clear_band(:,jcol,:) = 0.0_jprb
230            if (allocated(flux%sw_dn_direct_clear_band)) then
231              flux%sw_dn_direct_clear_band(:,jcol,:) = 0.0_jprb
232            end if
233          end if
234        end if
235
236        flux%sw_dn_diffuse_surf_g(:,jcol) = 0.0_jprb
237        flux%sw_dn_direct_surf_g(:,jcol)  = 0.0_jprb
238        if (config%do_clear) then
239          flux%sw_dn_diffuse_surf_clear_g(:,jcol) = 0.0_jprb
240          flux%sw_dn_direct_surf_clear_g(:,jcol)  = 0.0_jprb
241        end if
242
243        cycle
244      end if ! sun is below the horizon
245
246      ! At this point mu0 >= 1.0e-10
247
248      ! Define which layers contain cloud; assume that
249      ! cloud%crop_cloud_fraction has already been called
250      is_clear_sky_layer = .true.
251      do jlev = 1,nlev
252        if (cloud%fraction(jcol,jlev) > 0.0_jprb) then
253          is_clear_sky_layer(jlev) = .false.
254        end if
255      end do
256
257      ! --------------------------------------------------------
258      ! Section 3: Loop over layers to compute reflectance and transmittance
259      ! --------------------------------------------------------
260      ! In this section the reflectance, transmittance and sources
261      ! are computed for each layer
262      do jlev = 1,nlev ! Start at top-of-atmosphere
263
264        ! Array-wise assignments
265        gamma1 = 0.0_jprb
266        gamma2 = 0.0_jprb
267        gamma3 = 0.0_jprb
268
269        nreg = nregions
270        if (is_clear_sky_layer(jlev)) then
271          nreg = 1
272        end if
273        do jreg = 1,nreg
274          if (jreg == 1) then
275            od_total  =  od(:,jlev,jcol)
276            ssa_total = ssa(:,jlev,jcol)
277            g_total   =   g(:,jlev,jcol)
278          else
279            do jg = 1,ng
280              ! Mapping from g-point to band
281              iband = config%i_band_from_reordered_g_sw(jg)
282              scat_od = od(jg,jlev,jcol)*ssa(jg,jlev,jcol)
283              scat_od_cloud = od_cloud(iband,jlev,jcol) &
284                   &  * ssa_cloud(iband,jlev,jcol) * od_scaling(jreg,jlev,jcol)
285              ! Add scaled cloud optical depth to clear-sky value
286              od_total(jg) = od(jg,jlev,jcol) &
287                   &  + od_cloud(iband,jlev,jcol)*od_scaling(jreg,jlev,jcol)
288              ! Compute single-scattering albedo and asymmetry
289              ! factor of gas-cloud combination
290              ssa_total(jg) = (scat_od+scat_od_cloud) &
291                   &  / od_total(jg)
292              g_total(jg) = (scat_od*g(jg,jlev,jcol) &
293                   &         + scat_od_cloud * g_cloud(iband,jlev,jcol)) &
294                   &      / (scat_od + scat_od_cloud)
295            end do
296          end if
297
298          if (config%do_sw_delta_scaling_with_gases) then
299            ! Apply delta-Eddington scaling to the aerosol-gas(-cloud)
300            ! mixture
301            call delta_eddington(od_total, ssa_total, g_total)
302          end if
303          call calc_two_stream_gammas_sw(ng, &
304               &  mu0, ssa_total, g_total, &
305               &  gamma1, gamma2, gamma3)
306          call calc_reflectance_transmittance_sw(ng, &
307               &  mu0, od_total, ssa_total, &
308               &  gamma1, gamma2, gamma3, &
309               &  reflectance(:,jreg,jlev), transmittance(:,jreg,jlev), &
310               &  ref_dir(:,jreg,jlev), trans_dir_diff(:,jreg,jlev), &
311               &  trans_dir_dir(:,jreg,jlev) )
312        end do
313      end do
314
315      ! --------------------------------------------------------
316      ! Section 4: Compute total albedos
317      ! --------------------------------------------------------
318
319      total_albedo(:,:,:) = 0.0_jprb
320      total_albedo_direct(:,:,:) = 0.0_jprb
321
322      ! Copy surface albedo in clear-sky region
323      do jg = 1,ng
324        total_albedo(jg,1,nlev+1) = albedo_diffuse(jg,jcol)
325      end do
326
327      ! If direct albedo is available, use it; otherwise copy from
328      ! diffuse albedo
329      do jg = 1,ng
330        total_albedo_direct(jg,1,nlev+1) &
331             &  = mu0 * albedo_direct(jg,jcol)
332      end do
333
334      ! If there is cloud in the lowest layer then we need the albedos
335      ! underneath
336      if (.not. is_clear_sky_layer(nlev)) then
337        do jreg = 2,nregions
338          total_albedo(:,jreg,nlev+1)        = total_albedo(:,1,nlev+1)
339          total_albedo_direct(:,jreg,nlev+1) = total_albedo_direct(:,1,nlev+1)
340        end do
341      end if
342
343      if (config%do_clear) then
344        total_albedo_clear(:,nlev+1)        = total_albedo(:,1,nlev+1)
345        total_albedo_clear_direct(:,nlev+1) = total_albedo_direct(:,1,nlev+1)
346      end if
347
348      ! Work up from the surface computing the total albedo of the
349      ! atmosphere below that point using the adding method
350      do jlev = nlev,1,-1
351
352        total_albedo_below        = 0.0_jprb
353        total_albedo_below_direct = 0.0_jprb
354
355        if (config%do_clear) then
356          ! For clear-skies there is no need to consider "above" and
357          ! "below" quantities since with no cloud overlap to worry
358          ! about, these are the same
359          do jg = 1,ng
360            inv_denom(jg,1) = 1.0_jprb &
361                 &  / (1.0_jprb - total_albedo_clear(jg,jlev+1)*reflectance(jg,1,jlev))
362            total_albedo_clear(jg,jlev) = reflectance(jg,1,jlev) &
363                 &  + transmittance(jg,1,jlev) * transmittance(jg,1,jlev) &
364                 &  * total_albedo_clear(jg,jlev+1) * inv_denom(jg,1)
365            total_albedo_clear_direct(jg,jlev) = ref_dir(jg,1,jlev) &
366                 &  + (trans_dir_dir(jg,1,jlev)*total_albedo_clear_direct(jg,jlev+1) &
367                 &     +trans_dir_diff(jg,1,jlev)*total_albedo_clear(jg,jlev+1)) &
368                 &  * transmittance(jg,1,jlev) * inv_denom(jg,1)
369          end do
370        end if
371
372        if (is_clear_sky_layer(jlev)) then
373          do jg = 1,ng
374            inv_denom(jg,1) = 1.0_jprb &
375                 &  / (1.0_jprb - total_albedo(jg,1,jlev+1)*reflectance(jg,1,jlev))
376            total_albedo_below(jg,1) = reflectance(jg,1,jlev) &
377                 &  + transmittance(jg,1,jlev)  * transmittance(jg,1,jlev) &
378                 &  * total_albedo(jg,1,jlev+1) * inv_denom(jg,1)
379            total_albedo_below_direct(jg,1) = ref_dir(jg,1,jlev) &
380                 &  + (trans_dir_dir(jg,1,jlev)*total_albedo_direct(jg,1,jlev+1) &
381                 &     +trans_dir_diff(jg,1,jlev)*total_albedo(jg,1,jlev+1)) &
382                 &  * transmittance(jg,1,jlev) * inv_denom(jg,1)
383          end do
384        else
385          inv_denom = 1.0_jprb / (1.0_jprb - total_albedo(:,:,jlev+1)*reflectance(:,:,jlev))
386          total_albedo_below = reflectance(:,:,jlev) &
387               &  + transmittance(:,:,jlev)  * transmittance(:,:,jlev) &
388               &  * total_albedo(:,:,jlev+1) * inv_denom
389          total_albedo_below_direct = ref_dir(:,:,jlev) &
390               &  + (trans_dir_dir(:,:,jlev)*total_albedo_direct(:,:,jlev+1) &
391               &     +trans_dir_diff(:,:,jlev)*total_albedo(:,:,jlev+1)) &
392               &  * transmittance(:,:,jlev) * inv_denom
393        end if
394
395        ! Account for cloud overlap when converting albedo below a
396        ! layer interface to the equivalent values just above
397        if (is_clear_sky_layer(jlev) .and. is_clear_sky_layer(jlev-1)) then
398          total_albedo(:,:,jlev)        = total_albedo_below(:,:)
399          total_albedo_direct(:,:,jlev) = total_albedo_below_direct(:,:)
400        else
401          ! Use overlap matrix and exclude "anomalous" horizontal
402          ! transport described by Shonk & Hogan (2008).  Therefore,
403          ! the operation we perform is essentially diag(total_albedo)
404          ! = matmul(transpose(v_matrix)), diag(total_albedo_below)).
405          do jreg = 1,nregions
406            do jreg2 = 1,nregions
407              total_albedo(:,jreg,jlev) &
408                   &  = total_albedo(:,jreg,jlev) &
409                   &  + total_albedo_below(:,jreg2) &
410                   &  * v_matrix(jreg2,jreg,jlev,jcol)
411              total_albedo_direct(:,jreg,jlev) &
412                   &  = total_albedo_direct(:,jreg,jlev) &
413                   &  + total_albedo_below_direct(:,jreg2) &
414                   &  * v_matrix(jreg2,jreg,jlev,jcol)
415            end do
416          end do
417
418        end if
419
420      end do ! Reverse loop over levels
421
422      ! --------------------------------------------------------
423      ! Section 5: Compute fluxes
424      ! --------------------------------------------------------
425
426      ! Top-of-atmosphere fluxes into the regions of the top-most
427      ! layer, zero since we assume no diffuse downwelling
428      flux_dn = 0.0_jprb
429      ! Direct downwelling flux (into a plane perpendicular to the
430      ! sun) entering the top of each region in the top-most layer
431      do jreg = 1,nregions
432        direct_dn(:,jreg) = incoming_sw(:,jcol)*region_fracs(jreg,1,jcol)
433      end do
434      flux_up = direct_dn*total_albedo_direct(:,:,1)
435
436      if (config%do_clear) then
437        flux_dn_clear = 0.0_jprb
438        direct_dn_clear(:) = incoming_sw(:,jcol)
439        flux_up_clear = direct_dn_clear*total_albedo_clear_direct(:,1)
440      end if
441
442      ! Store the TOA broadband fluxes
443      flux%sw_up(jcol,1) = sum(sum(flux_up,1))
444      flux%sw_dn(jcol,1) = mu0 * sum(sum(direct_dn,1))
445      if (allocated(flux%sw_dn_direct)) then
446        flux%sw_dn_direct(jcol,1) = flux%sw_dn(jcol,1)
447      end if
448      if (config%do_clear) then
449        flux%sw_up_clear(jcol,1) = sum(flux_up_clear)
450        flux%sw_dn_clear(jcol,1) = mu0 * sum(direct_dn_clear)
451        if (allocated(flux%sw_dn_direct_clear)) then
452          flux%sw_dn_direct_clear(jcol,1) = flux%sw_dn_clear(jcol,1)
453        end if
454      end if
455
456      ! Save the spectral fluxes if required; some redundancy here as
457      ! the TOA downwelling flux is the same in clear and cloudy skies
458      if (config%do_save_spectral_flux) then
459        call indexed_sum(sum(flux_up,2), &
460             &           config%i_spec_from_reordered_g_sw, &
461             &           flux%sw_up_band(:,jcol,1))
462        call indexed_sum(sum(direct_dn,2), &
463             &           config%i_spec_from_reordered_g_sw, &
464             &           flux%sw_dn_band(:,jcol,1))
465        flux%sw_dn_band(:,jcol,1) = &
466             &  mu0 * flux%sw_dn_band(:,jcol,1)
467        if (allocated(flux%sw_dn_direct_band)) then
468          flux%sw_dn_direct_band(:,jcol,1) = flux%sw_dn_band(:,jcol,1)
469        end if
470        call add_indexed_sum(sum(flux_dn,2), &
471             &           config%i_spec_from_reordered_g_sw, &
472             &           flux%sw_dn_band(:,jcol,1))
473        if (config%do_clear) then
474          call indexed_sum(flux_up_clear, &
475               &           config%i_spec_from_reordered_g_sw, &
476               &           flux%sw_up_clear_band(:,jcol,1))
477          call indexed_sum(direct_dn_clear, &
478               &           config%i_spec_from_reordered_g_sw, &
479               &           flux%sw_dn_clear_band(:,jcol,1))
480          flux%sw_dn_clear_band(:,jcol,1) &
481               &  = mu0 * flux%sw_dn_clear_band(:,jcol,1)
482          if (allocated(flux%sw_dn_direct_clear_band)) then
483            flux%sw_dn_direct_clear_band(:,jcol,1) &
484                 &  = flux%sw_dn_clear_band(:,jcol,1)
485          end if
486          call add_indexed_sum(flux_dn_clear, &
487               &           config%i_spec_from_reordered_g_sw, &
488               &           flux%sw_dn_clear_band(:,jcol,1))
489        end if
490      end if
491
492      ! Final loop back down through the atmosphere to compute fluxes
493      do jlev = 1,nlev
494        if (config%do_clear) then
495          do jg = 1,ng
496            flux_dn_clear(jg) = (transmittance(jg,1,jlev)*flux_dn_clear(jg) + direct_dn_clear(jg) &
497               &  * (trans_dir_dir(jg,1,jlev)*total_albedo_clear_direct(jg,jlev+1)*reflectance(jg,1,jlev) &
498               &     + trans_dir_diff(jg,1,jlev) )) &
499               &  / (1.0_jprb - reflectance(jg,1,jlev)*total_albedo_clear(jg,jlev+1))
500            direct_dn_clear(jg) = trans_dir_dir(jg,1,jlev)*direct_dn_clear(jg)
501            flux_up_clear(jg) = direct_dn_clear(jg)*total_albedo_clear_direct(jg,jlev+1) &
502               &        +   flux_dn_clear(jg)*total_albedo_clear(jg,jlev+1)
503          end do
504        end if
505
506        if (is_clear_sky_layer(jlev)) then
507          do jg = 1,ng
508            flux_dn(jg,1) = (transmittance(jg,1,jlev)*flux_dn(jg,1) + direct_dn(jg,1) &
509                 &  * (trans_dir_dir(jg,1,jlev)*total_albedo_direct(jg,1,jlev+1)*reflectance(jg,1,jlev) &
510                 &     + trans_dir_diff(jg,1,jlev) )) &
511                 &  / (1.0_jprb - reflectance(jg,1,jlev)*total_albedo(jg,1,jlev+1))
512            direct_dn(jg,1) = trans_dir_dir(jg,1,jlev)*direct_dn(jg,1)
513            flux_up(jg,1) = direct_dn(jg,1)*total_albedo_direct(jg,1,jlev+1) &
514                 &  +        flux_dn(jg,1)*total_albedo(jg,1,jlev+1)
515          end do
516
517          flux_dn(:,2:)  = 0.0_jprb
518          flux_up(:,2:)  = 0.0_jprb
519          direct_dn(:,2:)= 0.0_jprb
520        else
521          flux_dn = (transmittance(:,:,jlev)*flux_dn + direct_dn &
522               &  * (trans_dir_dir(:,:,jlev)*total_albedo_direct(:,:,jlev+1)*reflectance(:,:,jlev) &
523               &     + trans_dir_diff(:,:,jlev) )) &
524               &  / (1.0_jprb - reflectance(:,:,jlev)*total_albedo(:,:,jlev+1))
525          direct_dn = trans_dir_dir(:,:,jlev)*direct_dn
526          flux_up = direct_dn*total_albedo_direct(:,:,jlev+1) &
527               &  +   flux_dn*total_albedo(:,:,jlev+1)
528        end if
529
530        if (.not. (is_clear_sky_layer(jlev) &
531             &    .and. is_clear_sky_layer(jlev+1))) then
532          ! Account for overlap rules in translating fluxes just above
533          ! a layer interface to the values just below
534          flux_dn_below = singlemat_x_vec(ng,ng,nregions, &
535               &  v_matrix(:,:,jlev+1,jcol), flux_dn)
536          direct_dn_below = singlemat_x_vec(ng,ng,nregions, &
537               &  v_matrix(:,:,jlev+1,jcol), direct_dn)
538          flux_dn = flux_dn_below
539          direct_dn = direct_dn_below
540        end if ! Otherwise the fluxes in each region are the same so
541               ! nothing to do
542
543        ! Store the broadband fluxes
544        flux%sw_up(jcol,jlev+1) = sum(sum(flux_up,1))
545        if (allocated(flux%sw_dn_direct)) then
546          flux%sw_dn_direct(jcol,jlev+1) = mu0 * sum(sum(direct_dn,1))
547          flux%sw_dn(jcol,jlev+1) &
548               &  = flux%sw_dn_direct(jcol,jlev+1) + sum(sum(flux_dn,1))
549        else
550          flux%sw_dn(jcol,jlev+1) = mu0 * sum(sum(direct_dn,1)) + sum(sum(flux_dn,1))   
551        end if
552        if (config%do_clear) then
553          flux%sw_up_clear(jcol,jlev+1) = sum(flux_up_clear)
554          if (allocated(flux%sw_dn_direct_clear)) then
555            flux%sw_dn_direct_clear(jcol,jlev+1) = mu0 * sum(direct_dn_clear)
556            flux%sw_dn_clear(jcol,jlev+1) &
557                 &  = flux%sw_dn_direct_clear(jcol,jlev+1) + sum(flux_dn_clear)
558          else
559            flux%sw_dn_clear(jcol,jlev+1) = mu0 * sum(direct_dn_clear) &
560                 &  + sum(flux_dn_clear)
561          end if
562        end if
563
564        ! Save the spectral fluxes if required
565        if (config%do_save_spectral_flux) then
566          call indexed_sum(sum(flux_up,2), &
567               &           config%i_spec_from_reordered_g_sw, &
568               &           flux%sw_up_band(:,jcol,jlev+1))
569          call indexed_sum(sum(direct_dn,2), &
570               &           config%i_spec_from_reordered_g_sw, &
571               &           flux%sw_dn_band(:,jcol,jlev+1))
572          flux%sw_dn_band(:,jcol,jlev+1) = &
573               &  mu0 * flux%sw_dn_band(:,jcol,jlev+1)
574          if (allocated(flux%sw_dn_direct_band)) then
575            flux%sw_dn_direct_band(:,jcol,jlev+1) &
576                 &  = flux%sw_dn_band(:,jcol,jlev+1)
577          end if
578          call add_indexed_sum(sum(flux_dn,2), &
579               &           config%i_spec_from_reordered_g_sw, &
580               &           flux%sw_dn_band(:,jcol,jlev+1))
581          if (config%do_clear) then
582            call indexed_sum(flux_up_clear, &
583                 &           config%i_spec_from_reordered_g_sw, &
584                 &           flux%sw_up_clear_band(:,jcol,jlev+1))
585            call indexed_sum(direct_dn_clear, &
586                 &           config%i_spec_from_reordered_g_sw, &
587                 &           flux%sw_dn_clear_band(:,jcol,jlev+1))
588            flux%sw_dn_clear_band(:,jcol,jlev+1) = &
589                 &  mu0 * flux%sw_dn_clear_band(:,jcol,jlev+1)
590            if (allocated(flux%sw_dn_direct_clear_band)) then
591              flux%sw_dn_direct_clear_band(:,jcol,jlev+1) &
592                   &  = flux%sw_dn_clear_band(:,jcol,jlev+1)
593            end if
594            call add_indexed_sum(flux_dn_clear, &
595                 &           config%i_spec_from_reordered_g_sw, &
596                 &           flux%sw_dn_clear_band(:,jcol,jlev+1))
597          end if
598        end if
599
600      end do ! Final loop over levels
601     
602      ! Store surface spectral fluxes, if required (after the end of
603      ! the final loop over levels, the current values of these arrays
604      ! will be the surface values)
605      flux%sw_dn_diffuse_surf_g(:,jcol) = sum(flux_dn,2)
606      flux%sw_dn_direct_surf_g(:,jcol)  = mu0 * sum(direct_dn,2)
607      if (config%do_clear) then
608        flux%sw_dn_diffuse_surf_clear_g(:,jcol) = flux_dn_clear
609        flux%sw_dn_direct_surf_clear_g(:,jcol)  = mu0 * direct_dn_clear
610      end if
611
612    end do ! Loop over columns
613
614    if (lhook) call dr_hook('radiation_tripleclouds_sw:solver_tripleclouds_sw',1,hook_handle)
615
616  end subroutine solver_tripleclouds_sw
617
618end module radiation_tripleclouds_sw
Note: See TracBrowser for help on using the repository browser.