source: LMDZ6/trunk/libf/phylmd/ecrad/radiation_mcica_lw.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: 15.7 KB
Line 
1! radiation_mcica_lw.F90 - Monte-Carlo Independent Column Approximation longtwave solver
2!
3! (C) Copyright 2015- 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-07-12  R. Hogan  Call fast adding method if only clouds scatter
19!   2017-10-23  R. Hogan  Renamed single-character variables
20
21module radiation_mcica_lw
22
23  public
24
25contains
26
27  !---------------------------------------------------------------------
28  ! Longwave Monte Carlo Independent Column Approximation
29  ! (McICA). This implementation performs a clear-sky and a cloudy-sky
30  ! calculation, and then weights the two to get the all-sky fluxes
31  ! according to the total cloud cover. This method reduces noise for
32  ! low cloud cover situations, and exploits the clear-sky
33  ! calculations that are usually performed for diagnostic purposes
34  ! simultaneously. The cloud generator has been carefully written
35  ! such that the stochastic cloud field satisfies the prescribed
36  ! overlap parameter accounting for this weighting.
37  subroutine solver_mcica_lw(nlev,istartcol,iendcol, &
38       &  config, single_level, cloud, &
39       &  od, ssa, g, od_cloud, ssa_cloud, g_cloud, planck_hl, &
40       &  emission, albedo, &
41       &  flux)
42
43    use parkind1, only           : jprb
44    use yomhook,  only           : lhook, dr_hook
45
46    use radiation_io,   only           : nulerr, radiation_abort
47    use radiation_config, only         : config_type
48    use radiation_single_level, only   : single_level_type
49    use radiation_cloud, only          : cloud_type
50    use radiation_flux, only           : flux_type
51    use radiation_two_stream, only     : calc_two_stream_gammas_lw, &
52         &                               calc_reflectance_transmittance_lw, &
53         &                               calc_no_scattering_transmittance_lw
54    use radiation_adding_ica_lw, only  : adding_ica_lw, fast_adding_ica_lw, &
55         &                               calc_fluxes_no_scattering_lw
56    use radiation_lw_derivatives, only : calc_lw_derivatives_ica, modify_lw_derivatives_ica
57    use radiation_cloud_generator, only: cloud_generator
58
59    implicit none
60
61    ! Inputs
62    integer, intent(in) :: nlev               ! number of model levels
63    integer, intent(in) :: istartcol, iendcol ! range of columns to process
64    type(config_type),        intent(in) :: config
65    type(single_level_type),  intent(in) :: single_level
66    type(cloud_type),         intent(in) :: cloud
67
68    ! Gas and aerosol optical depth, single-scattering albedo and
69    ! asymmetry factor at each longwave g-point
70    real(jprb), intent(in), dimension(config%n_g_lw, nlev, istartcol:iendcol) :: &
71         &  od
72    real(jprb), intent(in), dimension(config%n_g_lw_if_scattering, nlev, istartcol:iendcol) :: &
73         &  ssa, g
74
75    ! Cloud and precipitation optical depth, single-scattering albedo and
76    ! asymmetry factor in each longwave band
77    real(jprb), intent(in), dimension(config%n_bands_lw,nlev,istartcol:iendcol)   :: &
78         &  od_cloud
79    real(jprb), intent(in), dimension(config%n_bands_lw_if_scattering, &
80         &  nlev,istartcol:iendcol) :: ssa_cloud, g_cloud
81
82    ! Planck function at each half-level and the surface
83    real(jprb), intent(in), dimension(config%n_g_lw,nlev+1,istartcol:iendcol) :: &
84         &  planck_hl
85
86    ! Emission (Planck*emissivity) and albedo (1-emissivity) at the
87    ! surface at each longwave g-point
88    real(jprb), intent(in), dimension(config%n_g_lw, istartcol:iendcol) :: emission, albedo
89
90    ! Output
91    type(flux_type), intent(inout):: flux
92
93    ! Local variables
94
95    ! Diffuse reflectance and transmittance for each layer in clear
96    ! and all skies
97    real(jprb), dimension(config%n_g_lw, nlev) :: ref_clear, trans_clear, reflectance, transmittance
98
99    ! Emission by a layer into the upwelling or downwelling diffuse
100    ! streams, in clear and all skies
101    real(jprb), dimension(config%n_g_lw, nlev) :: source_up_clear, source_dn_clear, source_up, source_dn
102
103    ! Fluxes per g point
104    real(jprb), dimension(config%n_g_lw, nlev+1) :: flux_up, flux_dn
105    real(jprb), dimension(config%n_g_lw, nlev+1) :: flux_up_clear, flux_dn_clear
106
107    ! Combined gas+aerosol+cloud optical depth, single scattering
108    ! albedo and asymmetry factor
109    real(jprb), dimension(config%n_g_lw) :: od_total, ssa_total, g_total
110
111    ! Combined scattering optical depth
112    real(jprb) :: scat_od, scat_od_total(config%n_g_lw)
113
114    ! Two-stream coefficients
115    real(jprb), dimension(config%n_g_lw) :: gamma1, gamma2
116
117    ! Optical depth scaling from the cloud generator, zero indicating
118    ! clear skies
119    real(jprb), dimension(config%n_g_lw,nlev) :: od_scaling
120
121    ! Modified optical depth after McICA scaling to represent cloud
122    ! inhomogeneity
123    real(jprb), dimension(config%n_g_lw) :: od_cloud_new
124
125    ! Total cloud cover output from the cloud generator
126    real(jprb) :: total_cloud_cover
127
128    ! Identify clear-sky layers
129    logical :: is_clear_sky_layer(nlev)
130
131    ! Index of the highest cloudy layer
132    integer :: i_cloud_top
133
134    ! Number of g points
135    integer :: ng
136
137    ! Loop indices for level, column and g point
138    integer :: jlev, jcol, jg
139
140    real(jprb) :: hook_handle
141
142    if (lhook) call dr_hook('radiation_mcica_lw:solver_mcica_lw',0,hook_handle)
143
144    if (.not. config%do_clear) then
145      write(nulerr,'(a)') '*** Error: longwave McICA requires clear-sky calculation to be performed'
146      call radiation_abort()     
147    end if
148
149    ng = config%n_g_lw
150
151    ! Loop through columns
152    do jcol = istartcol,iendcol
153
154      ! Clear-sky calculation
155      if (config%do_lw_aerosol_scattering) then
156        ! Scattering case: first compute clear-sky reflectance,
157        ! transmittance etc at each model level
158        do jlev = 1,nlev
159          ssa_total = ssa(:,jlev,jcol)
160          g_total   = g(:,jlev,jcol)
161          call calc_two_stream_gammas_lw(ng, ssa_total, g_total, &
162               &  gamma1, gamma2)
163          call calc_reflectance_transmittance_lw(ng, &
164               &  od(:,jlev,jcol), gamma1, gamma2, &
165               &  planck_hl(:,jlev,jcol), planck_hl(:,jlev+1,jcol), &
166               &  ref_clear(:,jlev), trans_clear(:,jlev), &
167               &  source_up_clear(:,jlev), source_dn_clear(:,jlev))
168        end do
169        ! Then use adding method to compute fluxes
170        call adding_ica_lw(ng, nlev, &
171             &  ref_clear, trans_clear, source_up_clear, source_dn_clear, &
172             &  emission(:,jcol), albedo(:,jcol), &
173             &  flux_up_clear, flux_dn_clear)
174       
175      else
176        ! Non-scattering case: use simpler functions for
177        ! transmission and emission
178        do jlev = 1,nlev
179          call calc_no_scattering_transmittance_lw(ng, od(:,jlev,jcol), &
180               &  planck_hl(:,jlev,jcol), planck_hl(:,jlev+1, jcol), &
181               &  trans_clear(:,jlev), source_up_clear(:,jlev), source_dn_clear(:,jlev))
182        end do
183        ! Simpler down-then-up method to compute fluxes
184        call calc_fluxes_no_scattering_lw(ng, nlev, &
185             &  trans_clear, source_up_clear, source_dn_clear, &
186             &  emission(:,jcol), albedo(:,jcol), &
187             &  flux_up_clear, flux_dn_clear)
188       
189        ! Ensure that clear-sky reflectance is zero since it may be
190        ! used in cloudy-sky case
191        ref_clear = 0.0_jprb
192      end if
193
194      ! Sum over g-points to compute broadband fluxes
195      flux%lw_up_clear(jcol,:) = sum(flux_up_clear,1)
196      flux%lw_dn_clear(jcol,:) = sum(flux_dn_clear,1)
197      ! Store surface spectral downwelling fluxes
198      flux%lw_dn_surf_clear_g(:,jcol) = flux_dn_clear(:,nlev+1)
199
200      ! Do cloudy-sky calculation; add a prime number to the seed in
201      ! the longwave
202      call cloud_generator(ng, nlev, config%i_overlap_scheme, &
203           &  single_level%iseed(jcol) + 997, &
204           &  config%cloud_fraction_threshold, &
205           &  cloud%fraction(jcol,:), cloud%overlap_param(jcol,:), &
206           &  config%cloud_inhom_decorr_scaling, cloud%fractional_std(jcol,:), &
207           &  config%pdf_sampler, od_scaling, total_cloud_cover, &
208           &  is_beta_overlap=config%use_beta_overlap)
209     
210      ! Store total cloud cover
211      flux%cloud_cover_lw(jcol) = total_cloud_cover
212     
213      if (total_cloud_cover >= config%cloud_fraction_threshold) then
214        ! Total-sky calculation
215
216        is_clear_sky_layer = .true.
217        i_cloud_top = nlev+1
218        do jlev = 1,nlev
219          ! Compute combined gas+aerosol+cloud optical properties
220          if (cloud%fraction(jcol,jlev) >= config%cloud_fraction_threshold) then
221            is_clear_sky_layer(jlev) = .false.
222            ! Get index to the first cloudy layer from the top
223            if (i_cloud_top > jlev) then
224              i_cloud_top = jlev
225            end if
226
227            od_cloud_new = od_scaling(:,jlev) &
228                 &  * od_cloud(config%i_band_from_reordered_g_lw,jlev,jcol)
229            od_total = od(:,jlev,jcol) + od_cloud_new
230            ssa_total = 0.0_jprb
231            g_total   = 0.0_jprb
232
233            if (config%do_lw_cloud_scattering) then
234              ! Scattering case: calculate reflectance and
235              ! transmittance at each model level
236
237              if (config%do_lw_aerosol_scattering) then
238                ! In single precision we need to protect against the
239                ! case that od_total > 0.0 and ssa_total > 0.0 but
240                ! od_total*ssa_total == 0 due to underflow
241                scat_od_total = ssa(:,jlev,jcol)*od(:,jlev,jcol) &
242                     &     + ssa_cloud(config%i_band_from_reordered_g_lw,jlev,jcol) &
243                     &     *  od_cloud_new
244                where (scat_od_total > 0.0_jprb)
245                  g_total = (g(:,jlev,jcol)*ssa(:,jlev,jcol)*od(:,jlev,jcol) &
246                       &     +   g_cloud(config%i_band_from_reordered_g_lw,jlev,jcol) &
247                       &     * ssa_cloud(config%i_band_from_reordered_g_lw,jlev,jcol) &
248                       &     *  od_cloud_new) &
249                       &     / scat_od_total
250                end where               
251                where (od_total > 0.0_jprb)
252                  ssa_total = scat_od_total / od_total
253                end where
254              else
255                do jg = 1,ng
256                  if (od_total(jg) > 0.0_jprb) then
257                    scat_od = ssa_cloud(config%i_band_from_reordered_g_lw(jg),jlev,jcol) &
258                         &     * od_cloud_new(jg)
259                    ssa_total(jg) = scat_od / od_total(jg)
260                    if (scat_od > 0.0_jprb) then
261                      g_total(jg) = g_cloud(config%i_band_from_reordered_g_lw(jg),jlev,jcol) &
262                           &     * ssa_cloud(config%i_band_from_reordered_g_lw(jg),jlev,jcol) &
263                           &     *  od_cloud_new(jg) / scat_od
264                    end if
265                  end if
266                end do
267              end if
268           
269              ! Compute cloudy-sky reflectance, transmittance etc at
270              ! each model level
271              call calc_two_stream_gammas_lw(ng, ssa_total, g_total, &
272                   &  gamma1, gamma2)
273              call calc_reflectance_transmittance_lw(ng, &
274                   &  od_total, gamma1, gamma2, &
275                   &  planck_hl(:,jlev,jcol), planck_hl(:,jlev+1,jcol), &
276                   &  reflectance(:,jlev), transmittance(:,jlev), source_up(:,jlev), source_dn(:,jlev))
277            else
278              ! No-scattering case: use simpler functions for
279              ! transmission and emission
280              call calc_no_scattering_transmittance_lw(ng, od_total, &
281                   &  planck_hl(:,jlev,jcol), planck_hl(:,jlev+1, jcol), &
282                   &  transmittance(:,jlev), source_up(:,jlev), source_dn(:,jlev))
283            end if
284
285          else
286            ! Clear-sky layer: copy over clear-sky values
287            reflectance(:,jlev) = ref_clear(:,jlev)
288            transmittance(:,jlev) = trans_clear(:,jlev)
289            source_up(:,jlev) = source_up_clear(:,jlev)
290            source_dn(:,jlev) = source_dn_clear(:,jlev)
291          end if
292        end do
293       
294        if (config%do_lw_aerosol_scattering) then
295          ! Use adding method to compute fluxes for an overcast sky,
296          ! allowing for scattering in all layers
297          call adding_ica_lw(ng, nlev, reflectance, transmittance, source_up, source_dn, &
298               &  emission(:,jcol), albedo(:,jcol), &
299               &  flux_up, flux_dn)
300        else if (config%do_lw_cloud_scattering) then
301          ! Use adding method to compute fluxes but optimize for the
302          ! presence of clear-sky layers
303!          call adding_ica_lw(ng, nlev, reflectance, transmittance, source_up, source_dn, &
304!               &  emission(:,jcol), albedo(:,jcol), &
305!               &  flux_up, flux_dn)
306          call fast_adding_ica_lw(ng, nlev, reflectance, transmittance, source_up, source_dn, &
307               &  emission(:,jcol), albedo(:,jcol), &
308               &  is_clear_sky_layer, i_cloud_top, flux_dn_clear, &
309               &  flux_up, flux_dn)
310        else
311          ! Simpler down-then-up method to compute fluxes
312          call calc_fluxes_no_scattering_lw(ng, nlev, &
313               &  transmittance, source_up, source_dn, emission(:,jcol), albedo(:,jcol), &
314               &  flux_up, flux_dn)
315        end if
316       
317        ! Store overcast broadband fluxes
318        flux%lw_up(jcol,:) = sum(flux_up,1)
319        flux%lw_dn(jcol,:) = sum(flux_dn,1)
320
321        ! Cloudy flux profiles currently assume completely overcast
322        ! skies; perform weighted average with clear-sky profile
323        flux%lw_up(jcol,:) =  total_cloud_cover *flux%lw_up(jcol,:) &
324             &  + (1.0_jprb - total_cloud_cover)*flux%lw_up_clear(jcol,:)
325        flux%lw_dn(jcol,:) =  total_cloud_cover *flux%lw_dn(jcol,:) &
326             &  + (1.0_jprb - total_cloud_cover)*flux%lw_dn_clear(jcol,:)
327        ! Store surface spectral downwelling fluxes
328        flux%lw_dn_surf_g(:,jcol) = total_cloud_cover*flux_dn(:,nlev+1) &
329             &  + (1.0_jprb - total_cloud_cover)*flux%lw_dn_surf_clear_g(:,jcol)
330
331        ! Compute the longwave derivatives needed by Hogan and Bozzo
332        ! (2015) approximate radiation update scheme
333        if (config%do_lw_derivatives) then
334          call calc_lw_derivatives_ica(ng, nlev, jcol, transmittance, flux_up(:,nlev+1), &
335               &                       flux%lw_derivatives)
336          if (total_cloud_cover < 1.0_jprb - config%cloud_fraction_threshold) then
337            ! Modify the existing derivative with the contribution from the clear sky
338            call modify_lw_derivatives_ica(ng, nlev, jcol, trans_clear, flux_up_clear(:,nlev+1), &
339                 &                         1.0_jprb-total_cloud_cover, flux%lw_derivatives)
340          end if
341        end if
342
343      else
344        ! No cloud in profile and clear-sky fluxes already
345        ! calculated: copy them over
346        flux%lw_up(jcol,:) = flux%lw_up_clear(jcol,:)
347        flux%lw_dn(jcol,:) = flux%lw_dn_clear(jcol,:)
348        flux%lw_dn_surf_g(:,jcol) = flux%lw_dn_surf_clear_g(:,jcol)
349        if (config%do_lw_derivatives) then
350          call calc_lw_derivatives_ica(ng, nlev, jcol, trans_clear, flux_up_clear(:,nlev+1), &
351               &                       flux%lw_derivatives)
352 
353        end if
354      end if ! Cloud is present in profile
355    end do
356
357    if (lhook) call dr_hook('radiation_mcica_lw:solver_mcica_lw',1,hook_handle)
358   
359  end subroutine solver_mcica_lw
360
361end module radiation_mcica_lw
Note: See TracBrowser for help on using the repository browser.