source: LMDZ6/trunk/libf/phylmd/ecrad.v1.5.1/radiation_mcica_lw.F90 @ 5353

Last change on this file since 5353 was 4489, checked in by lguez, 21 months ago

Merge LMDZ_ECRad branch back into trunk!

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          call calc_two_stream_gammas_lw(ng, ssa(:,jlev,jcol), g(:,jlev,jcol), &
160               &  gamma1, gamma2)
161          call calc_reflectance_transmittance_lw(ng, &
162               &  od(:,jlev,jcol), gamma1, gamma2, &
163               &  planck_hl(:,jlev,jcol), planck_hl(:,jlev+1,jcol), &
164               &  ref_clear(:,jlev), trans_clear(:,jlev), &
165               &  source_up_clear(:,jlev), source_dn_clear(:,jlev))
166        end do
167        ! Then use adding method to compute fluxes
168        call adding_ica_lw(ng, nlev, &
169             &  ref_clear, trans_clear, source_up_clear, source_dn_clear, &
170             &  emission(:,jcol), albedo(:,jcol), &
171             &  flux_up_clear, flux_dn_clear)
172       
173      else
174        ! Non-scattering case: use simpler functions for
175        ! transmission and emission
176        do jlev = 1,nlev
177          call calc_no_scattering_transmittance_lw(ng, od(:,jlev,jcol), &
178               &  planck_hl(:,jlev,jcol), planck_hl(:,jlev+1, jcol), &
179               &  trans_clear(:,jlev), source_up_clear(:,jlev), source_dn_clear(:,jlev))
180        end do
181        ! Simpler down-then-up method to compute fluxes
182        call calc_fluxes_no_scattering_lw(ng, nlev, &
183             &  trans_clear, source_up_clear, source_dn_clear, &
184             &  emission(:,jcol), albedo(:,jcol), &
185             &  flux_up_clear, flux_dn_clear)
186       
187        ! Ensure that clear-sky reflectance is zero since it may be
188        ! used in cloudy-sky case
189        ref_clear = 0.0_jprb
190      end if
191
192      ! Sum over g-points to compute broadband fluxes
193      flux%lw_up_clear(jcol,:) = sum(flux_up_clear,1)
194      flux%lw_dn_clear(jcol,:) = sum(flux_dn_clear,1)
195      ! Store surface spectral downwelling fluxes
196      flux%lw_dn_surf_clear_g(:,jcol) = flux_dn_clear(:,nlev+1)
197
198      ! Do cloudy-sky calculation; add a prime number to the seed in
199      ! the longwave
200      call cloud_generator(ng, nlev, config%i_overlap_scheme, &
201           &  single_level%iseed(jcol) + 997, &
202           &  config%cloud_fraction_threshold, &
203           &  cloud%fraction(jcol,:), cloud%overlap_param(jcol,:), &
204           &  config%cloud_inhom_decorr_scaling, cloud%fractional_std(jcol,:), &
205           &  config%pdf_sampler, od_scaling, total_cloud_cover, &
206           &  use_beta_overlap=config%use_beta_overlap, &
207           &  use_vectorizable_generator=config%use_vectorizable_generator)
208     
209      ! Store total cloud cover
210      flux%cloud_cover_lw(jcol) = total_cloud_cover
211     
212      if (total_cloud_cover >= config%cloud_fraction_threshold) then
213        ! Total-sky calculation
214
215        is_clear_sky_layer = .true.
216        i_cloud_top = nlev+1
217        do jlev = 1,nlev
218          ! Compute combined gas+aerosol+cloud optical properties
219          if (cloud%fraction(jcol,jlev) >= config%cloud_fraction_threshold) then
220            is_clear_sky_layer(jlev) = .false.
221            ! Get index to the first cloudy layer from the top
222            if (i_cloud_top > jlev) then
223              i_cloud_top = jlev
224            end if
225
226            do jg = 1,ng
227              od_cloud_new(jg) = od_scaling(jg,jlev) &
228                 &  * od_cloud(config%i_band_from_reordered_g_lw(jg),jlev,jcol)
229              od_total(jg)  = od(jg,jlev,jcol) + od_cloud_new(jg)
230              ssa_total(jg) = 0.0_jprb
231              g_total(jg)   = 0.0_jprb
232            end do
233
234            if (config%do_lw_cloud_scattering) then
235              ! Scattering case: calculate reflectance and
236              ! transmittance at each model level
237
238              if (config%do_lw_aerosol_scattering) then
239                ! In single precision we need to protect against the
240                ! case that od_total > 0.0 and ssa_total > 0.0 but
241                ! od_total*ssa_total == 0 due to underflow
242                do jg = 1,ng
243                  if (od_total(jg) > 0.0_jprb) then
244                    scat_od_total(jg) = ssa(jg,jlev,jcol)*od(jg,jlev,jcol) &
245                     &     + ssa_cloud(config%i_band_from_reordered_g_lw(jg),jlev,jcol) &
246                     &     *  od_cloud_new(jg)
247                    ssa_total(jg) = scat_od_total(jg) / od_total(jg)
248
249                    if (scat_od_total(jg) > 0.0_jprb) then
250                      g_total(jg) = (g(jg,jlev,jcol)*ssa(jg,jlev,jcol)*od(jg,jlev,jcol) &
251                         &     +   g_cloud(config%i_band_from_reordered_g_lw(jg),jlev,jcol) &
252                         &     * ssa_cloud(config%i_band_from_reordered_g_lw(jg),jlev,jcol) &
253                         &     *  od_cloud_new(jg)) &
254                         &     / scat_od_total(jg)
255                    end if
256                  end if
257                end do
258
259              else
260
261                do jg = 1,ng
262                  if (od_total(jg) > 0.0_jprb) then
263                    scat_od = ssa_cloud(config%i_band_from_reordered_g_lw(jg),jlev,jcol) &
264                         &     * od_cloud_new(jg)
265                    ssa_total(jg) = scat_od / od_total(jg)
266                    if (scat_od > 0.0_jprb) then
267                      g_total(jg) = g_cloud(config%i_band_from_reordered_g_lw(jg),jlev,jcol) &
268                           &     * ssa_cloud(config%i_band_from_reordered_g_lw(jg),jlev,jcol) &
269                           &     *  od_cloud_new(jg) / scat_od
270                    end if
271                  end if
272                end do
273
274              end if
275           
276              ! Compute cloudy-sky reflectance, transmittance etc at
277              ! each model level
278              call calc_two_stream_gammas_lw(ng, ssa_total, g_total, &
279                   &  gamma1, gamma2)
280              call calc_reflectance_transmittance_lw(ng, &
281                   &  od_total, gamma1, gamma2, &
282                   &  planck_hl(:,jlev,jcol), planck_hl(:,jlev+1,jcol), &
283                   &  reflectance(:,jlev), transmittance(:,jlev), source_up(:,jlev), source_dn(:,jlev))
284            else
285              ! No-scattering case: use simpler functions for
286              ! transmission and emission
287              call calc_no_scattering_transmittance_lw(ng, od_total, &
288                   &  planck_hl(:,jlev,jcol), planck_hl(:,jlev+1, jcol), &
289                   &  transmittance(:,jlev), source_up(:,jlev), source_dn(:,jlev))
290            end if
291
292          else
293            ! Clear-sky layer: copy over clear-sky values
294            reflectance(:,jlev) = ref_clear(:,jlev)
295            transmittance(:,jlev) = trans_clear(:,jlev)
296            source_up(:,jlev) = source_up_clear(:,jlev)
297            source_dn(:,jlev) = source_dn_clear(:,jlev)
298          end if
299        end do
300       
301        if (config%do_lw_aerosol_scattering) then
302          ! Use adding method to compute fluxes for an overcast sky,
303          ! allowing for scattering in all layers
304          call adding_ica_lw(ng, nlev, reflectance, transmittance, source_up, source_dn, &
305               &  emission(:,jcol), albedo(:,jcol), &
306               &  flux_up, flux_dn)
307        else if (config%do_lw_cloud_scattering) then
308          ! Use adding method to compute fluxes but optimize for the
309          ! presence of clear-sky layers
310          call fast_adding_ica_lw(ng, nlev, reflectance, transmittance, source_up, source_dn, &
311               &  emission(:,jcol), albedo(:,jcol), &
312               &  is_clear_sky_layer, i_cloud_top, flux_dn_clear, &
313               &  flux_up, flux_dn)
314        else
315          ! Simpler down-then-up method to compute fluxes
316          call calc_fluxes_no_scattering_lw(ng, nlev, &
317               &  transmittance, source_up, source_dn, emission(:,jcol), albedo(:,jcol), &
318               &  flux_up, flux_dn)
319        end if
320       
321        ! Store overcast broadband fluxes
322        flux%lw_up(jcol,:) = sum(flux_up,1)
323        flux%lw_dn(jcol,:) = sum(flux_dn,1)
324
325        ! Cloudy flux profiles currently assume completely overcast
326        ! skies; perform weighted average with clear-sky profile
327        flux%lw_up(jcol,:) =  total_cloud_cover *flux%lw_up(jcol,:) &
328             &  + (1.0_jprb - total_cloud_cover)*flux%lw_up_clear(jcol,:)
329        flux%lw_dn(jcol,:) =  total_cloud_cover *flux%lw_dn(jcol,:) &
330             &  + (1.0_jprb - total_cloud_cover)*flux%lw_dn_clear(jcol,:)
331        ! Store surface spectral downwelling fluxes
332        flux%lw_dn_surf_g(:,jcol) = total_cloud_cover*flux_dn(:,nlev+1) &
333             &  + (1.0_jprb - total_cloud_cover)*flux%lw_dn_surf_clear_g(:,jcol)
334
335        ! Compute the longwave derivatives needed by Hogan and Bozzo
336        ! (2015) approximate radiation update scheme
337        if (config%do_lw_derivatives) then
338          call calc_lw_derivatives_ica(ng, nlev, jcol, transmittance, flux_up(:,nlev+1), &
339               &                       flux%lw_derivatives)
340          if (total_cloud_cover < 1.0_jprb - config%cloud_fraction_threshold) then
341            ! Modify the existing derivative with the contribution from the clear sky
342            call modify_lw_derivatives_ica(ng, nlev, jcol, trans_clear, flux_up_clear(:,nlev+1), &
343                 &                         1.0_jprb-total_cloud_cover, flux%lw_derivatives)
344          end if
345        end if
346
347      else
348        ! No cloud in profile and clear-sky fluxes already
349        ! calculated: copy them over
350        flux%lw_up(jcol,:) = flux%lw_up_clear(jcol,:)
351        flux%lw_dn(jcol,:) = flux%lw_dn_clear(jcol,:)
352        flux%lw_dn_surf_g(:,jcol) = flux%lw_dn_surf_clear_g(:,jcol)
353        if (config%do_lw_derivatives) then
354          call calc_lw_derivatives_ica(ng, nlev, jcol, trans_clear, flux_up_clear(:,nlev+1), &
355               &                       flux%lw_derivatives)
356 
357        end if
358      end if ! Cloud is present in profile
359    end do
360
361    if (lhook) call dr_hook('radiation_mcica_lw:solver_mcica_lw',1,hook_handle)
362   
363  end subroutine solver_mcica_lw
364
365end module radiation_mcica_lw
Note: See TracBrowser for help on using the repository browser.