source: LMDZ6/branches/contrails/libf/phylmd/ecrad/radiation/radiation_homogeneous_sw.F90 @ 5450

Last change on this file since 5450 was 4773, checked in by idelkadi, 13 months ago
  • Update of Ecrad in LMDZ The same organization of the Ecrad offline version is retained in order to facilitate the updating of Ecrad in LMDZ and the comparison between online and offline results. version 1.6.1 of Ecrad (https://github.com/lguez/ecrad.git)
  • Implementation of the double call of Ecrad in LMDZ


File size: 15.8 KB
Line 
1! radiation_homogeneous_sw.F90 - Shortwave homogeneous-column (no cloud fraction) 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!   2019-01-14  R. Hogan  Save spectral flux profile if required
20
21module radiation_homogeneous_sw
22
23  public
24
25contains
26
27  ! Provides elemental function "delta_eddington"
28#include "radiation_delta_eddington.h"
29
30  !---------------------------------------------------------------------
31  ! Shortwave homogeneous solver, in which clouds are assumed to fill
32  ! the gridbox horizontally
33  subroutine solver_homogeneous_sw(nlev,istartcol,iendcol, &
34       &  config, single_level, cloud, &
35       &  od, ssa, g, od_cloud, ssa_cloud, g_cloud, &
36       &  albedo_direct, albedo_diffuse, incoming_sw, &
37       &  flux)
38
39    use parkind1, only           : jprb
40    use yomhook,  only           : lhook, dr_hook, jphook
41
42    use radiation_config, only         : config_type
43    use radiation_single_level, only   : single_level_type
44    use radiation_cloud, only          : cloud_type
45    use radiation_flux, only           : flux_type, &
46         &                               indexed_sum_profile, add_indexed_sum_profile
47    use radiation_two_stream, only     : calc_two_stream_gammas_sw, &
48         &                       calc_reflectance_transmittance_sw
49    use radiation_constants, only      : Pi, GasConstantDryAir, &
50         &                               AccelDueToGravity
51    use radiation_adding_ica_sw, only  : adding_ica_sw
52
53    implicit none
54
55    ! Inputs
56    integer, intent(in) :: nlev               ! number of model levels
57    integer, intent(in) :: istartcol, iendcol ! range of columns to process
58    type(config_type),        intent(in) :: config
59    type(single_level_type),  intent(in) :: single_level
60    type(cloud_type),         intent(in) :: cloud
61
62    ! Gas and aerosol optical depth, single-scattering albedo and
63    ! asymmetry factor at each shortwave g-point
64    real(jprb), intent(in), dimension(config%n_g_sw, nlev, istartcol:iendcol) :: &
65         &  od, ssa, g
66
67    ! Cloud and precipitation optical depth, single-scattering albedo and
68    ! asymmetry factor in each shortwave band
69    real(jprb), intent(in), dimension(config%n_bands_sw,nlev,istartcol:iendcol)   :: &
70         &  od_cloud, ssa_cloud, g_cloud
71
72    ! Direct and diffuse surface albedos, and the incoming shortwave
73    ! flux into a plane perpendicular to the incoming radiation at
74    ! top-of-atmosphere in each of the shortwave g points
75    real(jprb), intent(in), dimension(config%n_g_sw,istartcol:iendcol) :: &
76         &  albedo_direct, albedo_diffuse, incoming_sw
77
78    ! Output
79    type(flux_type), intent(inout):: flux
80
81    ! Local variables
82
83    ! Cosine of solar zenith angle
84    real(jprb)                                 :: cos_sza
85
86    ! Diffuse reflectance and transmittance for each layer
87    real(jprb), dimension(config%n_g_sw, nlev) :: reflectance, transmittance
88
89    ! Fraction of direct beam scattered by a layer into the upwelling
90    ! or downwelling diffuse streams
91    real(jprb), dimension(config%n_g_sw, nlev) :: ref_dir, trans_dir_diff
92
93    ! Transmittance for the direct beam in clear and all skies
94    real(jprb), dimension(config%n_g_sw, nlev) :: trans_dir_dir
95
96    ! Fluxes per g point
97    real(jprb), dimension(config%n_g_sw, nlev+1) :: flux_up, flux_dn_diffuse, flux_dn_direct
98
99    ! Combined gas+aerosol+cloud optical depth, single scattering
100    ! albedo and asymmetry factor
101    real(jprb), dimension(config%n_g_sw) :: od_total, ssa_total, g_total
102
103    ! Two-stream coefficients
104    real(jprb), dimension(config%n_g_sw) :: gamma1, gamma2, gamma3
105
106    ! Optical depth of cloud in g-point space
107    real(jprb), dimension(config%n_g_sw) :: od_cloud_g
108
109    ! Is there any cloud in the profile?
110    logical :: is_cloudy_profile
111
112    ! Number of g points
113    integer :: ng
114
115    ! Loop indices for level and column
116    integer :: jlev, jcol
117
118    real(jphook) :: hook_handle
119
120    if (lhook) call dr_hook('radiation_homogeneous_sw:solver_homogeneous_sw',0,hook_handle)
121
122    ng = config%n_g_sw
123
124    ! Loop through columns
125    do jcol = istartcol,iendcol
126      ! Only perform calculation if sun above the horizon
127      if (single_level%cos_sza(jcol) > 0.0_jprb) then
128
129        cos_sza = single_level%cos_sza(jcol)
130       
131        ! Is there any cloud in the profile?
132        is_cloudy_profile = .false.
133        do jlev = 1,nlev
134          if (cloud%fraction(jcol,jlev) >= config%cloud_fraction_threshold) then
135            is_cloudy_profile = .true.
136            exit
137          end if
138        end do
139
140        ! If clear-sky fluxes need to be computed then we first
141        ! compute the reflectance and transmittance of all layers,
142        ! neglecting clouds. If clear-sky fluxes are not required then
143        ! we only do the clear-sky layers since these will be needed
144        ! when we come to do the total-sky fluxes.
145        if (.not. config%do_sw_delta_scaling_with_gases) then
146          ! Delta-Eddington scaling has already been performed to the
147          ! aerosol part of od, ssa and g
148          do jlev = 1,nlev
149            if (config%do_clear .or. cloud%fraction(jcol,jlev) &
150                 &                 < config%cloud_fraction_threshold) then
151              call calc_two_stream_gammas_sw(ng, cos_sza, &
152                   &  ssa(:,jlev,jcol), g(:,jlev,jcol), &
153                   &  gamma1, gamma2, gamma3)
154              call calc_reflectance_transmittance_sw(ng, &
155                   &  cos_sza, &
156                   &  od(:,jlev,jcol), ssa(:,jlev,jcol), &
157                   &  gamma1, gamma2, gamma3, &
158                   &  reflectance(:,jlev), transmittance(:,jlev), &
159                   &  ref_dir(:,jlev), trans_dir_diff(:,jlev), &
160                   &  trans_dir_dir(:,jlev) )
161             
162            end if
163          end do
164        else
165          ! Apply delta-Eddington scaling to the aerosol-gas mixture
166          do jlev = 1,nlev
167            if (config%do_clear .or. cloud%fraction(jcol,jlev) &
168                 &                 < config%cloud_fraction_threshold) then
169              od_total  =  od(:,jlev,jcol)
170              ssa_total = ssa(:,jlev,jcol)
171              g_total   =   g(:,jlev,jcol)
172              call delta_eddington(od_total, ssa_total, g_total)
173              call calc_two_stream_gammas_sw(ng, &
174                   &  cos_sza, ssa_total, g_total, &
175                   &  gamma1, gamma2, gamma3)
176              call calc_reflectance_transmittance_sw(ng, &
177                   &  cos_sza, od_total, ssa_total, &
178                   &  gamma1, gamma2, gamma3, &
179                   &  reflectance(:,jlev), transmittance(:,jlev), &
180                   &  ref_dir(:,jlev), trans_dir_diff(:,jlev), &
181                   &  trans_dir_dir(:,jlev) )
182            end if
183          end do
184        end if
185         
186        if (config%do_clear) then
187          ! Use adding method to compute fluxes
188          call adding_ica_sw(ng, nlev, incoming_sw(:,jcol), &
189               &  albedo_diffuse(:,jcol), albedo_direct(:,jcol), &
190               &  spread(cos_sza,1,ng), reflectance, transmittance, ref_dir, trans_dir_diff, &
191               &  trans_dir_dir, flux_up, flux_dn_diffuse, flux_dn_direct)
192       
193          ! Sum over g-points to compute and save clear-sky broadband
194          ! fluxes
195          flux%sw_up_clear(jcol,:) = sum(flux_up,1)
196          if (allocated(flux%sw_dn_direct_clear)) then
197            flux%sw_dn_direct_clear(jcol,:) &
198                 &  = sum(flux_dn_direct,1)
199            flux%sw_dn_clear(jcol,:) = sum(flux_dn_diffuse,1) &
200                 &  + flux%sw_dn_direct_clear(jcol,:)
201          else
202            flux%sw_dn_clear(jcol,:) = sum(flux_dn_diffuse,1) &
203                 &  + sum(flux_dn_direct,1)
204          end if
205          ! Store spectral downwelling fluxes at surface
206          flux%sw_dn_diffuse_surf_clear_g(:,jcol) = flux_dn_diffuse(:,nlev+1)
207          flux%sw_dn_direct_surf_clear_g(:,jcol)  = flux_dn_direct(:,nlev+1)
208
209          ! Save the spectral fluxes if required
210          if (config%do_save_spectral_flux) then
211            call indexed_sum_profile(flux_up, config%i_spec_from_reordered_g_sw, &
212                 &                   flux%sw_up_clear_band(:,jcol,:))
213            call indexed_sum_profile(flux_dn_direct, config%i_spec_from_reordered_g_sw, &
214                 &                   flux%sw_dn_clear_band(:,jcol,:))
215            if (allocated(flux%sw_dn_direct_clear_band)) then
216              flux%sw_dn_direct_clear_band(:,jcol,:) &
217                   &  = flux%sw_dn_clear_band(:,jcol,:)
218            end if
219            call add_indexed_sum_profile(flux_dn_diffuse, &
220                 &                       config%i_spec_from_reordered_g_sw, &
221                 &                       flux%sw_dn_clear_band(:,jcol,:))
222          end if
223
224        end if ! Do clear-sky calculations
225 
226        ! Now the total-sky calculation.  If this is a clear profile
227        ! and clear-sky fluxes have been calculated then we can simply
228        ! copy over the clear-sky fluxes, otherwise we need to compute
229        ! fluxes now.
230        if (is_cloudy_profile .or. .not. config%do_clear) then
231          do jlev = 1,nlev
232            ! Compute combined gas+aerosol+cloud optical properties;
233            ! note that for clear layers, the reflectance and
234            ! transmittance have already been calculated
235            if (cloud%fraction(jcol,jlev) >= config%cloud_fraction_threshold) then
236              od_cloud_g = od_cloud(config%i_band_from_reordered_g_sw,jlev,jcol)
237              od_total  = od(:,jlev,jcol) + od_cloud_g
238              ssa_total = 0.0_jprb
239              g_total   = 0.0_jprb
240              where (od_total > 0.0_jprb)
241                ssa_total = (ssa(:,jlev,jcol)*od(:,jlev,jcol) &
242                     &     + ssa_cloud(config%i_band_from_reordered_g_sw,jlev,jcol) &
243                     &     *  od_cloud_g) &
244                     &     / od_total
245              end where
246              where (ssa_total > 0.0_jprb .and. od_total > 0.0_jprb)
247                g_total = (g(:,jlev,jcol)*ssa(:,jlev,jcol)*od(:,jlev,jcol) &
248                     &     +   g_cloud(config%i_band_from_reordered_g_sw,jlev,jcol) &
249                     &     * ssa_cloud(config%i_band_from_reordered_g_sw,jlev,jcol) &
250                     &     *  od_cloud_g) &
251                     &     / (ssa_total*od_total)
252              end where
253
254              ! Apply delta-Eddington scaling to the cloud-aerosol-gas
255              ! mixture
256              if (config%do_sw_delta_scaling_with_gases) then
257                call delta_eddington(od_total, ssa_total, g_total)
258              end if
259
260              ! Compute cloudy-sky reflectance, transmittance etc at
261              ! each model level
262              call calc_two_stream_gammas_sw(ng, &
263                   &  cos_sza, ssa_total, g_total, &
264                   &  gamma1, gamma2, gamma3)
265              call calc_reflectance_transmittance_sw(ng, &
266                   &  cos_sza, od_total, ssa_total, &
267                   &  gamma1, gamma2, gamma3, &
268                   &  reflectance(:,jlev), transmittance(:,jlev), &
269                   &  ref_dir(:,jlev), trans_dir_diff(:,jlev), &
270                   &  trans_dir_dir(:,jlev) )
271
272            end if
273          end do
274           
275          ! Use adding method to compute fluxes for an overcast sky
276          call adding_ica_sw(ng, nlev, incoming_sw(:,jcol), &
277               &  albedo_diffuse(:,jcol), albedo_direct(:,jcol), &
278               &  spread(cos_sza,1,ng), reflectance, transmittance, ref_dir, trans_dir_diff, &
279               &  trans_dir_dir, flux_up, flux_dn_diffuse, flux_dn_direct)
280
281          ! Store overcast broadband fluxes
282          flux%sw_up(jcol,:) = sum(flux_up,1)
283          if (allocated(flux%sw_dn_direct)) then
284            flux%sw_dn_direct(jcol,:) = sum(flux_dn_direct,1)
285            flux%sw_dn(jcol,:) = sum(flux_dn_diffuse,1) &
286                 &  + flux%sw_dn_direct(jcol,:)
287          else
288            flux%sw_dn(jcol,:) = sum(flux_dn_diffuse,1) &
289                 &  + sum(flux_dn_direct,1)
290          end if
291
292          ! Likewise for surface spectral fluxes
293          flux%sw_dn_diffuse_surf_g(:,jcol) = flux_dn_diffuse(:,nlev+1)
294          flux%sw_dn_direct_surf_g(:,jcol)  = flux_dn_direct(:,nlev+1)
295
296          ! Save the spectral fluxes if required
297          if (config%do_save_spectral_flux) then
298            call indexed_sum_profile(flux_up, config%i_spec_from_reordered_g_sw, &
299                 &                   flux%sw_up_band(:,jcol,:))
300            call indexed_sum_profile(flux_dn_direct, config%i_spec_from_reordered_g_sw, &
301                 &                   flux%sw_dn_band(:,jcol,:))
302            if (allocated(flux%sw_dn_direct_band)) then
303              flux%sw_dn_direct_band(:,jcol,:) &
304                   &  = flux%sw_dn_band(:,jcol,:)
305            end if
306            call add_indexed_sum_profile(flux_dn_diffuse, &
307                 &                       config%i_spec_from_reordered_g_sw, &
308                 &                       flux%sw_dn_band(:,jcol,:))
309          end if
310
311        else
312          ! No cloud in profile and clear-sky fluxes already
313          ! calculated: copy them over
314          flux%sw_up(jcol,:) = flux%sw_up_clear(jcol,:)
315          flux%sw_dn(jcol,:) = flux%sw_dn_clear(jcol,:)
316          if (allocated(flux%sw_dn_direct)) then
317            flux%sw_dn_direct(jcol,:) = flux%sw_dn_direct_clear(jcol,:)
318          end if
319          flux%sw_dn_diffuse_surf_g(:,jcol) = flux%sw_dn_diffuse_surf_clear_g(:,jcol)
320          flux%sw_dn_direct_surf_g(:,jcol)  = flux%sw_dn_direct_surf_clear_g(:,jcol)
321
322          if (config%do_save_spectral_flux) then
323            flux%sw_up_band(:,jcol,:) = flux%sw_up_clear_band(:,jcol,:)
324            flux%sw_dn_band(:,jcol,:) = flux%sw_dn_clear_band(:,jcol,:)
325            if (allocated(flux%sw_dn_direct_band)) then
326              flux%sw_dn_direct_band(:,jcol,:) = flux%sw_dn_direct_clear_band(:,jcol,:)
327            end if
328          end if
329
330        end if ! Cloud is present in profile
331
332      else
333        ! Set fluxes to zero if sun is below the horizon
334        flux%sw_up(jcol,:) = 0.0_jprb
335        flux%sw_dn(jcol,:) = 0.0_jprb
336        if (allocated(flux%sw_dn_direct)) then
337          flux%sw_dn_direct(jcol,:) = 0.0_jprb
338        end if
339        flux%sw_dn_diffuse_surf_g(:,jcol) = 0.0_jprb
340        flux%sw_dn_direct_surf_g(:,jcol)  = 0.0_jprb
341
342        if (config%do_clear) then
343          flux%sw_up_clear(jcol,:) = 0.0_jprb
344          flux%sw_dn_clear(jcol,:) = 0.0_jprb
345          if (allocated(flux%sw_dn_direct_clear)) then
346            flux%sw_dn_direct_clear(jcol,:) = 0.0_jprb
347          end if
348          flux%sw_dn_diffuse_surf_clear_g(:,jcol) = 0.0_jprb
349          flux%sw_dn_direct_surf_clear_g(:,jcol)  = 0.0_jprb
350        end if
351
352        if (config%do_save_spectral_flux) then
353          flux%sw_dn_band(:,jcol,:) = 0.0_jprb
354          flux%sw_up_band(:,jcol,:) = 0.0_jprb
355          if (allocated(flux%sw_dn_direct_band)) then
356            flux%sw_dn_direct_band(:,jcol,:) = 0.0_jprb
357          end if
358          if (config%do_clear) then
359            flux%sw_dn_clear_band(:,jcol,:) = 0.0_jprb
360            flux%sw_up_clear_band(:,jcol,:) = 0.0_jprb
361            if (allocated(flux%sw_dn_direct_clear_band)) then
362              flux%sw_dn_direct_clear_band(:,jcol,:) = 0.0_jprb
363            end if
364          end if
365        end if
366
367      end if ! sun above horizon
368    end do
369
370    if (lhook) call dr_hook('radiation_homogeneous_sw:solver_homogeneous_sw',1,hook_handle)
371
372  end subroutine solver_homogeneous_sw
373
374end module radiation_homogeneous_sw
Note: See TracBrowser for help on using the repository browser.