source: LMDZ6/trunk/libf/phylmd/ecrad/radiation_cloudless_sw.F90 @ 4190

Last change on this file since 4190 was 3908, checked in by idelkadi, 3 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: 9.6 KB
Line 
1! radiation_cloudless_sw.F90 - Shortwave homogeneous cloudless solver
2!
3! (C) Copyright 2019- 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
16module radiation_cloudless_sw
17
18public :: solver_cloudless_sw
19
20contains
21
22  ! Provides elemental function "delta_eddington"
23#include "radiation_delta_eddington.h"
24
25  !---------------------------------------------------------------------
26  ! Shortwave homogeneous solver containing no clouds
27  subroutine solver_cloudless_sw(nlev,istartcol,iendcol, &
28       &  config, single_level, &
29       &  od, ssa, g, albedo_direct, albedo_diffuse, incoming_sw, &
30       &  flux)
31
32    use parkind1, only           : jprb
33    use yomhook,  only           : lhook, dr_hook
34
35    use radiation_config, only         : config_type
36    use radiation_single_level, only   : single_level_type
37    use radiation_flux, only           : flux_type, indexed_sum_profile, &
38         &                               add_indexed_sum_profile
39    use radiation_two_stream, only     : calc_two_stream_gammas_sw, &
40         &                       calc_reflectance_transmittance_sw
41    use radiation_constants, only      : Pi, GasConstantDryAir, &
42         &                               AccelDueToGravity
43    use radiation_adding_ica_sw, only  : adding_ica_sw
44
45    implicit none
46
47    ! Inputs
48    integer, intent(in) :: nlev               ! number of model levels
49    integer, intent(in) :: istartcol, iendcol ! range of columns to process
50    type(config_type),        intent(in) :: config
51    type(single_level_type),  intent(in) :: single_level
52
53    ! Gas and aerosol optical depth, single-scattering albedo and
54    ! asymmetry factor at each shortwave g-point
55    real(jprb), intent(in), dimension(config%n_g_sw, nlev, istartcol:iendcol) :: &
56         &  od, ssa, g
57
58    ! Direct and diffuse surface albedos, and the incoming shortwave
59    ! flux into a plane perpendicular to the incoming radiation at
60    ! top-of-atmosphere in each of the shortwave g points
61    real(jprb), intent(in), dimension(config%n_g_sw,istartcol:iendcol) :: &
62         &  albedo_direct, albedo_diffuse, incoming_sw
63
64    ! Output
65    type(flux_type), intent(inout):: flux
66
67    ! Local variables
68
69    ! Cosine of solar zenith angle
70    real(jprb)                                 :: cos_sza
71
72    ! Diffuse reflectance and transmittance for each layer
73    real(jprb), dimension(config%n_g_sw, nlev) :: reflectance, transmittance
74
75    ! Fraction of direct beam scattered by a layer into the upwelling
76    ! or downwelling diffuse streams
77    real(jprb), dimension(config%n_g_sw, nlev) :: ref_dir, trans_dir_diff
78
79    ! Transmittance for the direct beam in clear and all skies
80    real(jprb), dimension(config%n_g_sw, nlev) :: trans_dir_dir
81
82    ! Fluxes per g point
83    real(jprb), dimension(config%n_g_sw, nlev+1) :: flux_up, flux_dn_diffuse, flux_dn_direct
84
85    ! Combined optical depth, single scattering albedo and asymmetry
86    ! factor
87    real(jprb), dimension(config%n_g_sw) :: od_total, ssa_total, g_total
88
89    ! Two-stream coefficients
90    real(jprb), dimension(config%n_g_sw) :: gamma1, gamma2, gamma3
91
92    ! Number of g points
93    integer :: ng
94
95    ! Loop indices for level and column
96    integer :: jlev, jcol
97
98    real(jprb) :: hook_handle
99
100    if (lhook) call dr_hook('radiation_cloudless_sw:solver_cloudless_sw',0,hook_handle)
101
102    ng = config%n_g_sw
103
104    ! Loop through columns
105    do jcol = istartcol,iendcol
106      ! Only perform calculation if sun above the horizon
107      if (single_level%cos_sza(jcol) > 0.0_jprb) then
108
109        cos_sza = single_level%cos_sza(jcol)
110       
111        ! The following is the same as the clear-sky part of
112        ! solver_homogeneous_sw
113        if (.not. config%do_sw_delta_scaling_with_gases) then
114          ! Delta-Eddington scaling has already been performed to the
115          ! aerosol part of od, ssa and g
116          do jlev = 1,nlev
117            call calc_two_stream_gammas_sw(ng, cos_sza, &
118                 &  ssa(:,jlev,jcol), g(:,jlev,jcol), &
119                 &  gamma1, gamma2, gamma3)
120            call calc_reflectance_transmittance_sw(ng, &
121                 &  cos_sza, &
122                 &  od(:,jlev,jcol), ssa(:,jlev,jcol), &
123                 &  gamma1, gamma2, gamma3, &
124                 &  reflectance(:,jlev), transmittance(:,jlev), &
125                 &  ref_dir(:,jlev), trans_dir_diff(:,jlev), &
126                 &  trans_dir_dir(:,jlev) )
127          end do
128        else
129          ! Apply delta-Eddington scaling to the aerosol-gas mixture
130          do jlev = 1,nlev
131            od_total  =  od(:,jlev,jcol)
132            ssa_total = ssa(:,jlev,jcol)
133            g_total   =   g(:,jlev,jcol)
134            call delta_eddington(od_total, ssa_total, g_total)
135            call calc_two_stream_gammas_sw(ng, &
136                 &  cos_sza, ssa_total, g_total, &
137                 &  gamma1, gamma2, gamma3)
138            call calc_reflectance_transmittance_sw(ng, &
139                 &  cos_sza, od_total, ssa_total, &
140                 &  gamma1, gamma2, gamma3, &
141                 &  reflectance(:,jlev), transmittance(:,jlev), &
142                 &  ref_dir(:,jlev), trans_dir_diff(:,jlev), &
143                 &  trans_dir_dir(:,jlev) )
144          end do
145        end if
146         
147        ! Use adding method to compute fluxes
148        call adding_ica_sw(ng, nlev, incoming_sw(:,jcol), &
149             &  albedo_diffuse(:,jcol), albedo_direct(:,jcol), &
150             &  spread(cos_sza,1,ng), reflectance, transmittance, ref_dir, trans_dir_diff, &
151             &  trans_dir_dir, flux_up, flux_dn_diffuse, flux_dn_direct)
152       
153        ! Sum over g-points to compute and save clear-sky broadband
154        ! fluxes
155        flux%sw_up(jcol,:) = sum(flux_up,1)
156        if (allocated(flux%sw_dn_direct)) then
157          flux%sw_dn_direct(jcol,:) = sum(flux_dn_direct,1)
158          flux%sw_dn(jcol,:) = sum(flux_dn_diffuse,1) &
159               &  + flux%sw_dn_direct(jcol,:)
160        else
161          flux%sw_dn(jcol,:) = sum(flux_dn_diffuse,1) + sum(flux_dn_direct,1)
162        end if
163        ! Store spectral downwelling fluxes at surface
164        flux%sw_dn_diffuse_surf_g(:,jcol) = flux_dn_diffuse(:,nlev+1)
165        flux%sw_dn_direct_surf_g(:,jcol)  = flux_dn_direct(:,nlev+1)
166
167        ! Save the spectral fluxes if required
168        if (config%do_save_spectral_flux) then
169          call indexed_sum_profile(flux_up, config%i_spec_from_reordered_g_sw, &
170               &                   flux%sw_up_band(:,jcol,:))
171          call indexed_sum_profile(flux_dn_direct, config%i_spec_from_reordered_g_sw, &
172               &                   flux%sw_dn_band(:,jcol,:))
173          if (allocated(flux%sw_dn_direct_band)) then
174            flux%sw_dn_direct_band(:,jcol,:) &
175                 &  = flux%sw_dn_band(:,jcol,:)
176          end if
177          call add_indexed_sum_profile(flux_dn_diffuse, &
178               &                       config%i_spec_from_reordered_g_sw, &
179               &                       flux%sw_dn_band(:,jcol,:))
180        end if
181
182        if (config%do_clear) then
183          ! Clear-sky calculations are equal to all-sky for this
184          ! solver: copy fluxes over
185          flux%sw_up_clear(jcol,:) = flux%sw_up(jcol,:)
186          flux%sw_dn_clear(jcol,:) = flux%sw_dn(jcol,:)
187          if (allocated(flux%sw_dn_direct_clear)) then
188            flux%sw_dn_direct_clear(jcol,:) = flux%sw_dn_direct(jcol,:)
189          end if
190          flux%sw_dn_diffuse_surf_clear_g(:,jcol) = flux%sw_dn_diffuse_surf_g(:,jcol)
191          flux%sw_dn_direct_surf_clear_g(:,jcol)  = flux%sw_dn_direct_surf_g(:,jcol)
192
193          if (config%do_save_spectral_flux) then
194            flux%sw_up_clear_band(:,jcol,:) = flux%sw_up_band(:,jcol,:)
195            flux%sw_dn_clear_band(:,jcol,:) = flux%sw_dn_band(:,jcol,:)
196            if (allocated(flux%sw_dn_direct_clear_band)) then
197              flux%sw_dn_direct_clear_band(:,jcol,:) = flux%sw_dn_direct_band(:,jcol,:)
198            end if
199          end if
200
201        end if ! do_clear
202
203      else
204        ! Set fluxes to zero if sun is below the horizon
205        flux%sw_up(jcol,:) = 0.0_jprb
206        flux%sw_dn(jcol,:) = 0.0_jprb
207        if (allocated(flux%sw_dn_direct)) then
208          flux%sw_dn_direct(jcol,:) = 0.0_jprb
209        end if
210        flux%sw_dn_diffuse_surf_g(:,jcol) = 0.0_jprb
211        flux%sw_dn_direct_surf_g(:,jcol)  = 0.0_jprb
212
213        if (config%do_clear) then
214          flux%sw_up_clear(jcol,:) = 0.0_jprb
215          flux%sw_dn_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          flux%sw_dn_diffuse_surf_clear_g(:,jcol) = 0.0_jprb
220          flux%sw_dn_direct_surf_clear_g(:,jcol)  = 0.0_jprb
221        end if
222
223        if (config%do_save_spectral_flux) then
224          flux%sw_dn_band(:,jcol,:) = 0.0_jprb
225          flux%sw_up_band(:,jcol,:) = 0.0_jprb
226          if (allocated(flux%sw_dn_direct_band)) then
227            flux%sw_dn_direct_band(:,jcol,:) = 0.0_jprb
228          end if
229          if (config%do_clear) then
230            flux%sw_dn_clear_band(:,jcol,:) = 0.0_jprb
231            flux%sw_up_clear_band(:,jcol,:) = 0.0_jprb
232            if (allocated(flux%sw_dn_direct_clear_band)) then
233              flux%sw_dn_direct_clear_band(:,jcol,:) = 0.0_jprb
234            end if
235          end if
236        end if
237
238      end if ! sun above horizon
239    end do
240
241    if (lhook) call dr_hook('radiation_cloudless_sw:solver_cloudless_sw',1,hook_handle)
242
243  end subroutine solver_cloudless_sw
244
245end module radiation_cloudless_sw
Note: See TracBrowser for help on using the repository browser.