source: LMDZ6/branches/cirrus/libf/phylmd/ecrad/radiation/radiation_homogeneous_lw.F90

Last change on this file was 4773, checked in by idelkadi, 11 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: 13.1 KB
Line 
1! radiation_homogeneous_lw.F90 - Longwave 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 emission/albedo rather than planck/emissivity
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_lw
22
23  public
24
25contains
26
27  !---------------------------------------------------------------------
28  ! Longwave homogeneous solver, in which clouds are assumed to fill
29  ! the gridbox horizontally
30  subroutine solver_homogeneous_lw(nlev,istartcol,iendcol, &
31       &  config, cloud, &
32       &  od, ssa, g, od_cloud, ssa_cloud, g_cloud, planck_hl, &
33       &  emission, albedo, &
34       &  flux)
35
36    use parkind1, only           : jprb
37    use yomhook,  only           : lhook, dr_hook, jphook
38
39    use radiation_config, only         : config_type
40    use radiation_cloud, only          : cloud_type
41    use radiation_flux, only           : flux_type, indexed_sum_profile
42    use radiation_two_stream, only     : calc_two_stream_gammas_lw, &
43         &                               calc_reflectance_transmittance_lw, &
44         &                               calc_no_scattering_transmittance_lw
45    use radiation_adding_ica_lw, only  : adding_ica_lw, calc_fluxes_no_scattering_lw
46    use radiation_lw_derivatives, only : calc_lw_derivatives_ica
47
48    implicit none
49
50    ! Inputs
51    integer, intent(in) :: nlev               ! number of model levels
52    integer, intent(in) :: istartcol, iendcol ! range of columns to process
53    type(config_type),        intent(in) :: config
54    type(cloud_type),         intent(in) :: cloud
55
56    ! Gas and aerosol optical depth, single-scattering albedo and
57    ! asymmetry factor at each longwave g-point
58    real(jprb), intent(in), dimension(config%n_g_lw, nlev, istartcol:iendcol) :: &
59         &  od
60    real(jprb), intent(in), dimension(config%n_g_lw_if_scattering, nlev, istartcol:iendcol) :: &
61         &  ssa, g
62
63    ! Cloud and precipitation optical depth, single-scattering albedo and
64    ! asymmetry factor in each longwave band
65    real(jprb), intent(in), dimension(config%n_bands_lw,nlev,istartcol:iendcol)   :: &
66         &  od_cloud
67    real(jprb), intent(in), dimension(config%n_bands_lw_if_scattering, &
68         &  nlev,istartcol:iendcol) :: ssa_cloud, g_cloud
69
70    ! Planck function at each half-level and the surface
71    real(jprb), intent(in), dimension(config%n_g_lw,nlev+1,istartcol:iendcol) :: &
72         &  planck_hl
73 
74    ! Emission (Planck*emissivity) and albedo (1-emissivity) at the
75    ! surface at each longwave g-point
76    real(jprb), intent(in), dimension(config%n_g_lw, istartcol:iendcol) &
77         &  :: emission, albedo
78
79    ! Output
80    type(flux_type), intent(inout):: flux
81
82    ! Local variables
83
84    ! Diffuse reflectance and transmittance for each layer in clear
85    ! and all skies
86    real(jprb), dimension(config%n_g_lw, nlev) :: reflectance, transmittance
87
88    ! Emission by a layer into the upwelling or downwelling diffuse
89    ! streams, in clear and all skies
90    real(jprb), dimension(config%n_g_lw, nlev) :: source_up, source_dn
91
92    ! Fluxes per g point
93    real(jprb), dimension(config%n_g_lw, nlev+1) :: flux_up, flux_dn
94
95    ! Combined gas+aerosol+cloud optical depth, single scattering
96    ! albedo and asymmetry factor
97    real(jprb), dimension(config%n_g_lw) :: od_total, ssa_total, g_total
98
99    ! Two-stream coefficients
100    real(jprb), dimension(config%n_g_lw) :: gamma1, gamma2
101
102    ! Optical depth of cloud in g-point space
103    real(jprb), dimension(config%n_g_lw) :: od_cloud_g
104
105    ! Is there any cloud in the profile?
106    logical :: is_cloudy_profile
107
108    ! Number of g points
109    integer :: ng
110
111    ! Loop indices for level and column
112    integer :: jlev, jcol
113
114    real(jphook) :: hook_handle
115
116    if (lhook) call dr_hook('radiation_homogeneous_lw:solver_homogeneous_lw',0,hook_handle)
117
118    ng = config%n_g_lw
119
120    ! Loop through columns
121    do jcol = istartcol,iendcol
122
123      ! Is there any cloud in the profile?
124      is_cloudy_profile = .false.
125      do jlev = 1,nlev
126        if (cloud%fraction(jcol,jlev) >= config%cloud_fraction_threshold) then
127          is_cloudy_profile = .true.
128          exit
129        end if
130      end do
131
132      ! If clear-sky fluxes need to be computed then we first compute
133      ! the reflectance and transmittance of all layers, neglecting
134      ! clouds. If clear-sky fluxes are not required then we only do
135      ! the clear-sky layers since these will be needed when we come
136      ! to do the total-sky fluxes.
137      do jlev = 1,nlev
138        if (config%do_clear .or. cloud%fraction(jcol,jlev) &
139             &                 < config%cloud_fraction_threshold) then
140          if (config%do_lw_aerosol_scattering) then
141            ! Scattering case: first compute clear-sky reflectance,
142            ! transmittance etc at each model level
143            ssa_total = ssa(:,jlev,jcol)
144            g_total   = g(:,jlev,jcol)
145            call calc_two_stream_gammas_lw(ng, ssa_total, g_total, &
146                 &  gamma1, gamma2)
147            call calc_reflectance_transmittance_lw(ng, &
148                 &  od(:,jlev,jcol), gamma1, gamma2, &
149                 &  planck_hl(:,jlev,jcol), planck_hl(:,jlev+1,jcol), &
150                 &  reflectance(:,jlev), transmittance(:,jlev), &
151                 &  source_up(:,jlev), source_dn(:,jlev))
152          else
153            ! Non-scattering case: use simpler functions for
154            ! transmission and emission
155            call calc_no_scattering_transmittance_lw(ng, od(:,jlev,jcol), &
156                 &  planck_hl(:,jlev,jcol), planck_hl(:,jlev+1, jcol), &
157                 &  transmittance(:,jlev), source_up(:,jlev), source_dn(:,jlev))         
158            ! Ensure that clear-sky reflectance is zero since it may be
159            ! used in cloudy-sky case
160            reflectance(:,jlev) = 0.0_jprb
161          end if
162         
163        end if
164      end do
165
166      if (config%do_clear) then
167        if (config%do_lw_aerosol_scattering) then
168          ! Then use adding method to compute fluxes
169          call adding_ica_lw(ng, nlev, &
170               &  reflectance, transmittance, source_up, source_dn, &
171               &  emission(:,jcol), albedo(:,jcol), &
172               &  flux_up, flux_dn)
173        else
174          ! Simpler down-then-up method to compute fluxes
175          call calc_fluxes_no_scattering_lw(ng, nlev, &
176               &  transmittance, source_up, source_dn, &
177               &  emission(:,jcol), albedo(:,jcol), &
178               &  flux_up, flux_dn)
179         
180        end if
181
182        ! Sum over g-points to compute broadband fluxes
183        flux%lw_up_clear(jcol,:) = sum(flux_up,1)
184        flux%lw_dn_clear(jcol,:) = sum(flux_dn,1)
185        ! Store surface spectral downwelling fluxes
186        flux%lw_dn_surf_clear_g(:,jcol) = flux_dn(:,nlev+1)
187
188        ! Save the spectral fluxes if required
189        if (config%do_save_spectral_flux) then
190          call indexed_sum_profile(flux_up, config%i_spec_from_reordered_g_lw, &
191               &                   flux%lw_up_clear_band(:,jcol,:))
192          call indexed_sum_profile(flux_dn, config%i_spec_from_reordered_g_lw, &
193               &                   flux%lw_dn_clear_band(:,jcol,:))
194        end if
195
196      end if ! Do clear-sky calculations
197
198      ! Now the total-sky calculation.  If this is a clear profile and
199      ! clear-sky fluxes have been calculated then we can simply copy
200      ! over the clear-sky fluxes, otherwise we need to compute fluxes
201      ! now.
202      if (is_cloudy_profile .or. .not. config%do_clear) then
203        do jlev = 1,nlev
204          ! Compute combined gas+aerosol+cloud optical properties;
205          ! note that for clear layers, the reflectance and
206          ! transmittance have already been calculated
207          if (cloud%fraction(jcol,jlev) >= config%cloud_fraction_threshold) then
208            od_cloud_g = od_cloud(config%i_band_from_reordered_g_lw,jlev,jcol)
209            od_total = od(:,jlev,jcol) + od_cloud_g
210            ssa_total = 0.0_jprb
211            g_total   = 0.0_jprb
212
213            if (config%do_lw_cloud_scattering) then
214              ! Scattering case: calculate reflectance and
215              ! transmittance at each model level
216              if (config%do_lw_aerosol_scattering) then
217                where (od_total > 0.0_jprb)
218                  ssa_total = (ssa(:,jlev,jcol)*od(:,jlev,jcol) &
219                       &     + ssa_cloud(config%i_band_from_reordered_g_lw,jlev,jcol) &
220                       &     *  od_cloud_g) &
221                       &     / od_total
222                end where
223                where (ssa_total > 0.0_jprb .and. od_total > 0.0_jprb)
224                  g_total = (g(:,jlev,jcol)*ssa(:,jlev,jcol)*od(:,jlev,jcol) &
225                       &     +   g_cloud(config%i_band_from_reordered_g_lw,jlev,jcol) &
226                       &     * ssa_cloud(config%i_band_from_reordered_g_lw,jlev,jcol) &
227                       &     *  od_cloud_g) &
228                       &     / (ssa_total*od_total)
229                end where
230              else
231                where (od_total > 0.0_jprb)
232                  ssa_total = ssa_cloud(config%i_band_from_reordered_g_lw,jlev,jcol) &
233                       &     * od_cloud_g / od_total
234                end where
235                where (ssa_total > 0.0_jprb .and. od_total > 0.0_jprb)
236                  g_total = g_cloud(config%i_band_from_reordered_g_lw,jlev,jcol) &
237                       &     * ssa_cloud(config%i_band_from_reordered_g_lw,jlev,jcol) &
238                       &     *  od_cloud_g / (ssa_total*od_total)
239                end where
240              end if
241           
242              ! Compute cloudy-sky reflectance, transmittance etc at
243              ! each model level
244              call calc_two_stream_gammas_lw(ng, ssa_total, g_total, &
245                   &  gamma1, gamma2)
246              call calc_reflectance_transmittance_lw(ng, &
247                   &  od_total, gamma1, gamma2, &
248                   &  planck_hl(:,jlev,jcol), planck_hl(:,jlev+1,jcol), &
249                   &  reflectance(:,jlev), transmittance(:,jlev), &
250                   &  source_up(:,jlev), source_dn(:,jlev))
251            else
252              ! No-scattering case: use simpler functions for
253              ! transmission and emission
254              call calc_no_scattering_transmittance_lw(ng, od_total, &
255                   &  planck_hl(:,jlev,jcol), planck_hl(:,jlev+1, jcol), &
256                   &  transmittance(:,jlev), source_up(:,jlev), source_dn(:,jlev))
257            end if
258          end if ! is cloudy layer
259        end do
260       
261        if (config%do_lw_cloud_scattering) then
262          ! Use adding method to compute fluxes for an overcast sky
263          call adding_ica_lw(ng, nlev, reflectance, transmittance, source_up, source_dn, &
264               &  emission(:,jcol), albedo(:,jcol), &
265               &  flux_up, flux_dn)
266        else
267          ! Simpler down-then-up method to compute fluxes
268          call calc_fluxes_no_scattering_lw(ng, nlev, &
269               &  transmittance, source_up, source_dn, emission(:,jcol), albedo(:,jcol), &
270               &  flux_up, flux_dn)
271        end if
272       
273        ! Store overcast broadband fluxes
274        flux%lw_up(jcol,:) = sum(flux_up,1)
275        flux%lw_dn(jcol,:) = sum(flux_dn,1)
276        ! Store surface spectral downwelling fluxes
277        flux%lw_dn_surf_g(:,jcol) = flux_dn(:,nlev+1)
278
279        ! Save the spectral fluxes if required
280        if (config%do_save_spectral_flux) then
281          call indexed_sum_profile(flux_up, config%i_spec_from_reordered_g_lw, &
282               &                   flux%lw_up_band(:,jcol,:))
283          call indexed_sum_profile(flux_dn, config%i_spec_from_reordered_g_lw, &
284               &                   flux%lw_dn_band(:,jcol,:))
285        end if
286
287      else
288        ! No cloud in profile and clear-sky fluxes already
289        ! calculated: copy them over
290        flux%lw_up(jcol,:) = flux%lw_up_clear(jcol,:)
291        flux%lw_dn(jcol,:) = flux%lw_dn_clear(jcol,:)
292        flux%lw_dn_surf_g(:,jcol) = flux%lw_dn_surf_clear_g(:,jcol)
293
294        if (config%do_save_spectral_flux) then
295          flux%lw_up_band(:,jcol,:) = flux%lw_up_clear_band(:,jcol,:)
296          flux%lw_dn_band(:,jcol,:) = flux%lw_dn_clear_band(:,jcol,:)
297        end if
298
299     end if
300
301      ! Compute the longwave derivatives needed by Hogan and Bozzo
302      ! (2015) approximate radiation update scheme, using clear-sky
303      ! transmittance if no clouds were present in the profile,
304      ! all-sky transmittance otherwise
305      if (config%do_lw_derivatives) then
306        call calc_lw_derivatives_ica(ng, nlev, jcol, transmittance, flux_up(:,nlev+1), &
307             &                       flux%lw_derivatives)
308 
309      end if
310
311    end do
312
313    if (lhook) call dr_hook('radiation_homogeneous_lw:solver_homogeneous_lw',1,hook_handle)
314   
315  end subroutine solver_homogeneous_lw
316
317end module radiation_homogeneous_lw
Note: See TracBrowser for help on using the repository browser.