source: LMDZ6/branches/blowing_snow/libf/phylmd/ecrad/radiation_cloud_optics.F90 @ 5411

Last change on this file since 5411 was 3908, checked in by idelkadi, 4 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: 22.7 KB
Line 
1! radiation_cloud_optics.F90 - Computing cloud optical properties
2!
3! (C) Copyright 2014- 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-07-22  R. Hogan  Added Yi et al. ice optics model
17
18module radiation_cloud_optics
19
20  implicit none
21  public
22
23contains
24
25  ! Provides elemental function "delta_eddington_scat_od"
26#include "radiation_delta_eddington.h"
27
28  !---------------------------------------------------------------------
29  ! Load cloud scattering data; this subroutine delegates to one
30  ! in radiation_cloud_optics_data.F90, but checks the size of
31  ! what is returned
32  subroutine setup_cloud_optics(config)
33
34    use parkind1,         only : jprb
35    use yomhook,          only : lhook, dr_hook
36
37    use radiation_io,     only : nulerr, radiation_abort
38    use radiation_config, only : config_type, IIceModelFu, IIceModelBaran, &
39         &                       IIceModelBaran2016, IIceModelBaran2017, &
40         &                       IIceModelYi, &
41         &                       ILiquidModelSOCRATES, ILiquidModelSlingo
42    use radiation_cloud_optics_data, only  : cloud_optics_type
43    use radiation_ice_optics_fu, only    : NIceOpticsCoeffsFuSW, &
44         &                                 NIceOpticsCoeffsFuLW
45    use radiation_ice_optics_baran, only : NIceOpticsCoeffsBaran, &
46         &                                 NIceOpticsCoeffsBaran2016
47    use radiation_ice_optics_baran2017, only : NIceOpticsCoeffsBaran2017, &
48         &                                 NIceOpticsGeneralCoeffsBaran2017
49    use radiation_ice_optics_yi, only    : NIceOpticsCoeffsYiSW, &
50         &                                 NIceOpticsCoeffsYiLW
51    use radiation_liquid_optics_socrates, only : NLiqOpticsCoeffsSOCRATES
52    use radiation_liquid_optics_slingo, only : NLiqOpticsCoeffsSlingoSW, &
53         &                                     NLiqOpticsCoeffsLindnerLiLW
54
55    type(config_type), intent(inout) :: config
56
57    real(jprb) :: hook_handle
58
59    if (lhook) call dr_hook('radiation_cloud_optics:setup_cloud_optics',0,hook_handle)
60
61    call config%cloud_optics%setup(trim(config%liq_optics_file_name), &
62         &     trim(config%ice_optics_file_name), iverbose=config%iverbosesetup)
63
64    ! Check liquid coefficients
65    if (size(config%cloud_optics%liq_coeff_lw, 1) /= config%n_bands_lw) then
66      write(nulerr,'(a,i0,a,i0,a)') &
67           &  '*** Error: number of longwave bands for droplets (', &
68           &  size(config%cloud_optics%liq_coeff_lw, 1), &
69           &  ') does not match number for gases (', config%n_bands_lw, ')'
70      call radiation_abort()
71    end if
72    if (size(config%cloud_optics%liq_coeff_sw, 1) /= config%n_bands_sw) then
73      write(nulerr,'(a,i0,a,i0,a)') &
74           &  '*** Error: number of shortwave bands for droplets (', &
75           &  size(config%cloud_optics%liq_coeff_sw, 1), &
76           &  ') does not match number for gases (', config%n_bands_sw, ')'
77      call radiation_abort()       
78    end if
79
80    if (config%i_liq_model == ILiquidModelSOCRATES) then
81      if (size(config%cloud_optics%liq_coeff_lw, 2) /= NLiqOpticsCoeffsSOCRATES) then
82        write(nulerr,'(a,i0,a,i0,a,i0,a)') &
83             &  '*** Error: number of liquid cloud optical coefficients (', &
84             &  size(config%cloud_optics%liq_coeff_lw, 2), &
85             &  ') does not match number expected (', NLiqOpticsCoeffsSOCRATES,')'
86        call radiation_abort()
87      end if
88    else if (config%i_liq_model == ILiquidModelSlingo) then
89      if (size(config%cloud_optics%liq_coeff_sw, 2) /= NLiqOpticsCoeffsSlingoSW) then
90        write(nulerr,'(a,i0,a,i0,a,i0,a)') &
91             &  '*** Error: number of shortwave liquid cloud optical coefficients (', &
92             &  size(config%cloud_optics%liq_coeff_sw, 2), &
93             &  ') does not match number expected (', NLiqOpticsCoeffsSlingoSW,')'
94        call radiation_abort()
95      end if
96      if (size(config%cloud_optics%liq_coeff_lw, 2) /= NLiqOpticsCoeffsLindnerLiLW) then
97        write(nulerr,'(a,i0,a,i0,a,i0,a)') &
98             &  '*** Error: number of longwave liquid cloud optical coefficients (', &
99             &  size(config%cloud_optics%liq_coeff_lw, 2), &
100             &  ') does not match number expected (', NLiqOpticsCoeffsLindnerLiLw,')'
101        call radiation_abort()
102      end if
103    end if
104
105    ! Check ice coefficients
106    if (size(config%cloud_optics%ice_coeff_lw, 1) /= config%n_bands_lw) then
107      write(nulerr,'(a,i0,a,i0,a)') &
108           &  '*** Error: number of longwave bands for ice particles (', &
109           &  size(config%cloud_optics%ice_coeff_lw, 1), &
110           &  ') does not match number for gases (', config%n_bands_lw, ')'
111      call radiation_abort()
112    end if
113    if (size(config%cloud_optics%ice_coeff_sw, 1) /= config%n_bands_sw) then
114      write(nulerr,'(a,i0,a,i0,a)') &
115           &  '*** Error: number of shortwave bands for ice particles (', &
116           &  size(config%cloud_optics%ice_coeff_sw, 1), &
117           &  ') does not match number for gases (', config%n_bands_sw, ')'
118      call radiation_abort()
119    end if
120
121    if (config%i_ice_model == IIceModelFu) then
122      if (size(config%cloud_optics%ice_coeff_lw, 2) &
123           &  /= NIceOpticsCoeffsFuLW) then
124        write(nulerr,'(a,i0,a,i0,a,i0,a)') &
125             &  '*** Error: number of LW ice-particle optical coefficients (', &
126             &  size(config%cloud_optics%ice_coeff_lw, 2), &
127             &  ') does not match number expected (', NIceOpticsCoeffsFuLW,')'
128        call radiation_abort()
129      end if
130      if (size(config%cloud_optics%ice_coeff_sw, 2) &
131           &  /= NIceOpticsCoeffsFuSW) then
132        write(nulerr,'(a,i0,a,i0,a,i0,a)') &
133             &  '*** Error: number of SW ice-particle optical coefficients (', &
134             &  size(config%cloud_optics%ice_coeff_sw, 2), &
135             &  ') does not match number expected (', NIceOpticsCoeffsFuSW,')'
136        call radiation_abort()
137      end if
138    else if (config%i_ice_model == IIceModelBaran &
139         &  .and. size(config%cloud_optics%ice_coeff_lw, 2) &
140         &  /= NIceOpticsCoeffsBaran) then
141      write(nulerr,'(a,i0,a,i0,a,i0,a)') &
142           &  '*** Error: number of ice-particle optical coefficients (', &
143           &  size(config%cloud_optics%ice_coeff_lw, 2), &
144           &  ') does not match number expected (', NIceOpticsCoeffsBaran,')'
145      call radiation_abort()
146    else if (config%i_ice_model == IIceModelBaran2016 &
147         &  .and. size(config%cloud_optics%ice_coeff_lw, 2) &
148         &  /= NIceOpticsCoeffsBaran2016) then
149      write(nulerr,'(a,i0,a,i0,a,i0,a)') &
150           &  '*** Error: number of ice-particle optical coefficients (', &
151           &  size(config%cloud_optics%ice_coeff_lw, 2), &
152           &  ') does not match number expected (', NIceOpticsCoeffsBaran2016,')'
153      call radiation_abort()
154    else if (config%i_ice_model == IIceModelBaran2017) then
155      if (size(config%cloud_optics%ice_coeff_lw, 2) &
156           &  /= NIceOpticsCoeffsBaran2017) then
157        write(nulerr,'(a,i0,a,i0,a,i0,a)') &
158             &  '*** Error: number of ice-particle optical coefficients (', &
159             &  size(config%cloud_optics%ice_coeff_lw, 2), &
160             &  ') does not match number expected (', NIceOpticsCoeffsBaran2017,')'
161        call radiation_abort()
162      else if (.not. allocated(config%cloud_optics%ice_coeff_gen)) then
163        write(nulerr,'(a)') &
164             &  '*** Error: coeff_gen needed for Baran-2017 ice optics parameterization'
165        call radiation_abort()
166      else if (size(config%cloud_optics%ice_coeff_gen) &
167           &  /= NIceOpticsGeneralCoeffsBaran2017) then
168        write(nulerr,'(a,i0,a,i0,a,i0,a)') &
169             &  '*** Error: number of general ice-particle optical coefficients (', &
170             &  size(config%cloud_optics%ice_coeff_gen), &
171             &  ') does not match number expected (', NIceOpticsGeneralCoeffsBaran2017,')'
172        call radiation_abort()
173      end if
174    else if (config%i_ice_model == IIceModelYi) then
175      if (size(config%cloud_optics%ice_coeff_lw, 2) &
176           &  /= NIceOpticsCoeffsYiLW) then
177        write(nulerr,'(a,i0,a,i0,a,i0,a)') &
178             &  '*** Error: number of LW ice-particle optical coefficients (', &
179             &  size(config%cloud_optics%ice_coeff_lw, 2), &
180             &  ') does not match number expected (', NIceOpticsCoeffsYiLW,')'
181        call radiation_abort()
182      end if
183      if (size(config%cloud_optics%ice_coeff_sw, 2) &
184           &  /= NIceOpticsCoeffsYiSW) then
185        write(nulerr,'(a,i0,a,i0,a,i0,a)') &
186             &  '*** Error: number of SW ice-particle optical coefficients (', &
187             &  size(config%cloud_optics%ice_coeff_sw, 2), &
188             &  ') does not match number expected (', NIceOpticsCoeffsYiSW,')'
189        call radiation_abort()
190      end if
191    end if
192
193    if (lhook) call dr_hook('radiation_cloud_optics:setup_cloud_optics',1,hook_handle)
194
195  end subroutine setup_cloud_optics
196
197
198  !---------------------------------------------------------------------
199  ! Compute cloud optical properties
200  subroutine cloud_optics(nlev,istartcol,iendcol, &
201       &  config, thermodynamics, cloud, &
202       &  od_lw_cloud, ssa_lw_cloud, g_lw_cloud, &
203       &  od_sw_cloud, ssa_sw_cloud, g_sw_cloud)
204
205    use parkind1, only           : jprb
206    use yomhook,  only           : lhook, dr_hook
207
208    use radiation_io,     only : nulout, nulerr, radiation_abort
209    use radiation_config, only : config_type, IIceModelFu, IIceModelBaran, &
210         &                       IIceModelBaran2016, IIceModelBaran2017, &
211         &                       IIceModelYi, &
212         &                       ILiquidModelSOCRATES, ILiquidModelSlingo
213    use radiation_thermodynamics, only    : thermodynamics_type
214    use radiation_cloud, only             : cloud_type
215    use radiation_constants, only         : AccelDueToGravity
216    use radiation_cloud_optics_data, only : cloud_optics_type
217    use radiation_ice_optics_fu, only     : calc_ice_optics_fu_sw, &
218         &                                  calc_ice_optics_fu_lw
219    use radiation_ice_optics_baran, only  : calc_ice_optics_baran, &
220         &                                  calc_ice_optics_baran2016
221    use radiation_ice_optics_baran2017, only  : calc_ice_optics_baran2017
222    use radiation_ice_optics_yi, only     : calc_ice_optics_yi_sw, &
223         &                                  calc_ice_optics_yi_lw
224    use radiation_liquid_optics_socrates, only:calc_liq_optics_socrates
225    use radiation_liquid_optics_slingo, only:calc_liq_optics_slingo, &
226         &                                   calc_liq_optics_lindner_li
227
228    integer, intent(in) :: nlev               ! number of model levels
229    integer, intent(in) :: istartcol, iendcol ! range of columns to process
230    type(config_type), intent(in), target :: config
231    type(thermodynamics_type),intent(in)  :: thermodynamics
232    type(cloud_type),   intent(in)        :: cloud
233
234    ! Layer optical depth, single scattering albedo and g factor of
235    ! clouds in each longwave band, where the latter two
236    ! variables are only defined if cloud longwave scattering is
237    ! enabled (otherwise both are treated as zero).
238    real(jprb), dimension(config%n_bands_lw,nlev,istartcol:iendcol), intent(out) :: &
239         &   od_lw_cloud
240    real(jprb), dimension(config%n_bands_lw_if_scattering,nlev,istartcol:iendcol), &
241         &   intent(out) :: ssa_lw_cloud, g_lw_cloud
242
243    ! Layer optical depth, single scattering albedo and g factor of
244    ! clouds in each shortwave band
245    real(jprb), dimension(config%n_bands_sw,nlev,istartcol:iendcol), intent(out) :: &
246         &   od_sw_cloud, ssa_sw_cloud, g_sw_cloud
247
248    ! Longwave and shortwave optical depth, scattering optical depth
249    ! and asymmetry factor, for liquid and ice in all bands but a
250    ! single cloud layer
251    real(jprb), dimension(config%n_bands_lw) :: &
252         &  od_lw_liq, scat_od_lw_liq, g_lw_liq, &
253         &  od_lw_ice, scat_od_lw_ice, g_lw_ice
254    real(jprb), dimension(config%n_bands_sw) :: &
255         &  od_sw_liq, scat_od_sw_liq, g_sw_liq, &
256         &  od_sw_ice, scat_od_sw_ice, g_sw_ice
257
258    ! In-cloud water path of cloud liquid or ice (i.e. liquid or ice
259    ! gridbox-mean water path divided by cloud fraction); kg m-2
260    real(jprb) :: lwp_in_cloud, iwp_in_cloud
261
262    ! Full-level temperature (K)
263    real(jprb) :: temperature
264
265    ! Factor to convert gridbox-mean mixing ratio to in-cloud water
266    ! path
267    real(jprb) :: factor
268
269    ! Pointer to the cloud optics coefficients for brevity of
270    ! access
271    type(cloud_optics_type), pointer :: ho
272
273    integer    :: jcol, jlev
274
275    real(jprb) :: hook_handle
276
277    if (lhook) call dr_hook('radiation_cloud_optics:cloud_optics',0,hook_handle)
278
279    if (config%iverbose >= 2) then
280      write(nulout,'(a)') 'Computing cloud absorption/scattering properties'
281    end if
282
283    ho => config%cloud_optics
284
285    ! Array-wise assignment
286    od_lw_cloud  = 0.0_jprb
287    od_sw_cloud  = 0.0_jprb
288    ssa_sw_cloud = 0.0_jprb
289    g_sw_cloud   = 0.0_jprb
290    if (config%do_lw_cloud_scattering) then
291      ssa_lw_cloud = 0.0_jprb
292      g_lw_cloud   = 0.0_jprb
293    end if
294
295    do jlev = 1,nlev
296      do jcol = istartcol,iendcol
297        ! Only do anything if cloud is present (assume that
298        ! cloud%crop_cloud_fraction has already been called)
299        if (cloud%fraction(jcol,jlev) > 0.0_jprb) then
300
301          ! Compute in-cloud liquid and ice water path
302          if (config%is_homogeneous) then
303            ! Homogeneous solvers assume cloud fills the box
304            ! horizontally, so we don't divide by cloud fraction
305            factor = ( thermodynamics%pressure_hl(jcol,jlev+1)    &
306                 &  -thermodynamics%pressure_hl(jcol,jlev  )  ) &
307                 &  / AccelDueToGravity
308          else
309            factor = ( thermodynamics%pressure_hl(jcol,jlev+1)    &
310                 &  -thermodynamics%pressure_hl(jcol,jlev  )  ) &
311                 &  / (AccelDueToGravity * cloud%fraction(jcol,jlev))
312          end if
313          lwp_in_cloud = factor * cloud%q_liq(jcol,jlev)
314          iwp_in_cloud = factor * cloud%q_ice(jcol,jlev)
315
316          ! Only compute liquid properties if liquid cloud is
317          ! present
318          if (lwp_in_cloud > 0.0_jprb) then
319            if (config%i_liq_model == ILiquidModelSOCRATES) then
320              ! Compute longwave properties
321              call calc_liq_optics_socrates(config%n_bands_lw, &
322                   &  config%cloud_optics%liq_coeff_lw, &
323                   &  lwp_in_cloud, cloud%re_liq(jcol,jlev), &
324                   &  od_lw_liq, scat_od_lw_liq, g_lw_liq)
325              ! Compute shortwave properties
326              call calc_liq_optics_socrates(config%n_bands_sw, &
327                   &  config%cloud_optics%liq_coeff_sw, &
328                   &  lwp_in_cloud, cloud%re_liq(jcol,jlev), &
329                   &  od_sw_liq, scat_od_sw_liq, g_sw_liq)
330            else if (config%i_liq_model == ILiquidModelSlingo) then
331              ! Compute longwave properties
332              call calc_liq_optics_lindner_li(config%n_bands_lw, &
333                   &  config%cloud_optics%liq_coeff_lw, &
334                   &  lwp_in_cloud, cloud%re_liq(jcol,jlev), &
335                   &  od_lw_liq, scat_od_lw_liq, g_lw_liq)
336              ! Compute shortwave properties
337              call calc_liq_optics_slingo(config%n_bands_sw, &
338                   &  config%cloud_optics%liq_coeff_sw, &
339                   &  lwp_in_cloud, cloud%re_liq(jcol,jlev), &
340                   &  od_sw_liq, scat_od_sw_liq, g_sw_liq)
341            else
342              write(nulerr,*) '*** Error: Unknown liquid model with code', &
343                   &          config%i_liq_model
344              call radiation_abort()
345            end if
346
347            if (.not. config%do_sw_delta_scaling_with_gases) then
348              ! Delta-Eddington scaling in the shortwave only
349              call delta_eddington_scat_od(od_sw_liq, scat_od_sw_liq, g_sw_liq)
350            end if
351          else
352            ! Liquid not present: set properties to zero
353            od_lw_liq = 0.0_jprb
354            scat_od_lw_liq = 0.0_jprb
355            g_lw_liq = 0.0_jprb
356
357            od_sw_liq = 0.0_jprb
358            scat_od_sw_liq = 0.0_jprb
359            g_sw_liq = 0.0_jprb
360          end if ! Liquid present
361
362          ! Only compute ice properties if ice cloud is present
363          if (iwp_in_cloud > 0.0_jprb) then
364            if (config%i_ice_model == IIceModelBaran) then
365              ! Compute longwave properties
366              call calc_ice_optics_baran(config%n_bands_lw, &
367                   &  config%cloud_optics%ice_coeff_lw, &
368                   &  iwp_in_cloud, cloud%q_ice(jcol,jlev), &
369                   &  od_lw_ice, scat_od_lw_ice, g_lw_ice)
370              ! Compute shortwave properties
371              call calc_ice_optics_baran(config%n_bands_sw, &
372                   &  config%cloud_optics%ice_coeff_sw, &
373                   &  iwp_in_cloud, cloud%q_ice(jcol,jlev), &
374                   &  od_sw_ice, scat_od_sw_ice, g_sw_ice)
375            else if (config%i_ice_model == IIceModelBaran2016) then
376              temperature = 0.5_jprb * (thermodynamics%temperature_hl(jcol,jlev) &
377                   &                   +thermodynamics%temperature_hl(jcol,jlev+1))
378              ! Compute longwave properties
379              call calc_ice_optics_baran2016(config%n_bands_lw, &
380                   &  config%cloud_optics%ice_coeff_lw, &
381                   &  iwp_in_cloud, cloud%q_ice(jcol,jlev), &
382                   &  temperature, &
383                   &  od_lw_ice, scat_od_lw_ice, g_lw_ice)
384              ! Compute shortwave properties
385              call calc_ice_optics_baran2016(config%n_bands_sw, &
386                   &  config%cloud_optics%ice_coeff_sw, &
387                   &  iwp_in_cloud, cloud%q_ice(jcol,jlev), &
388                   &  temperature, &
389                   &  od_sw_ice, scat_od_sw_ice, g_sw_ice)
390            else if (config%i_ice_model == IIceModelBaran2017) then
391              temperature = 0.5_jprb * (thermodynamics%temperature_hl(jcol,jlev) &
392                   &                   +thermodynamics%temperature_hl(jcol,jlev+1))
393              ! Compute longwave properties
394              call calc_ice_optics_baran2017(config%n_bands_lw, &
395                   &  config%cloud_optics%ice_coeff_gen, &
396                   &  config%cloud_optics%ice_coeff_lw, &
397                   &  iwp_in_cloud, cloud%q_ice(jcol,jlev), &
398                   &  temperature, &
399                   &  od_lw_ice, scat_od_lw_ice, g_lw_ice)
400              ! Compute shortwave properties
401              call calc_ice_optics_baran2017(config%n_bands_sw, &
402                   &  config%cloud_optics%ice_coeff_gen, &
403                   &  config%cloud_optics%ice_coeff_sw, &
404                   &  iwp_in_cloud, cloud%q_ice(jcol,jlev), &
405                   &  temperature, &
406                   &  od_sw_ice, scat_od_sw_ice, g_sw_ice)
407            else if (config%i_ice_model == IIceModelFu) then
408              ! Compute longwave properties
409              call calc_ice_optics_fu_lw(config%n_bands_lw, &
410                   &  config%cloud_optics%ice_coeff_lw, &
411                   &  iwp_in_cloud, cloud%re_ice(jcol,jlev), &
412                   &  od_lw_ice, scat_od_lw_ice, g_lw_ice)
413              if (config%do_fu_lw_ice_optics_bug) then
414                ! Reproduce bug in old IFS scheme
415                scat_od_lw_ice = od_lw_ice - scat_od_lw_ice
416              end if
417              ! Compute shortwave properties
418              call calc_ice_optics_fu_sw(config%n_bands_sw, &
419                   &  config%cloud_optics%ice_coeff_sw, &
420                   &  iwp_in_cloud, cloud%re_ice(jcol,jlev), &
421                   &  od_sw_ice, scat_od_sw_ice, g_sw_ice)
422            else if (config%i_ice_model == IIceModelYi) then
423              ! Compute longwave properties
424              call calc_ice_optics_yi_lw(config%n_bands_lw, &
425                   &  config%cloud_optics%ice_coeff_lw, &
426                   &  iwp_in_cloud, cloud%re_ice(jcol,jlev), &
427                   &  od_lw_ice, scat_od_lw_ice, g_lw_ice)
428              ! Compute shortwave properties
429              call calc_ice_optics_yi_sw(config%n_bands_sw, &
430                   &  config%cloud_optics%ice_coeff_sw, &
431                   &  iwp_in_cloud, cloud%re_ice(jcol,jlev), &
432                   &  od_sw_ice, scat_od_sw_ice, g_sw_ice)
433            else
434              write(nulerr,*) '*** Error: Unknown ice model with code', &
435                   &          config%i_ice_model
436              call radiation_abort()
437            end if
438
439            if (.not. config%do_sw_delta_scaling_with_gases) then
440              ! Delta-Eddington scaling in both longwave and shortwave
441              ! (assume that particles are larger than wavelength even
442              ! in longwave)
443              call delta_eddington_scat_od(od_sw_ice, scat_od_sw_ice, g_sw_ice)
444            end if
445
446            call delta_eddington_scat_od(od_lw_ice, scat_od_lw_ice, g_lw_ice)
447          else
448            ! Ice not present: set properties to zero
449            od_lw_ice = 0.0_jprb
450            scat_od_lw_ice = 0.0_jprb
451            g_lw_ice = 0.0_jprb
452
453            od_sw_ice = 0.0_jprb
454            scat_od_sw_ice = 0.0_jprb
455            g_sw_ice = 0.0_jprb
456          end if ! Ice present
457
458          ! Combine liquid and ice
459          if (config%do_lw_cloud_scattering) then
460            od_lw_cloud(:,jlev,jcol) = od_lw_liq + od_lw_ice
461            where (scat_od_lw_liq+scat_od_lw_ice > 0.0_jprb)
462              g_lw_cloud(:,jlev,jcol) = (g_lw_liq * scat_od_lw_liq &
463                   &  + g_lw_ice * scat_od_lw_ice) &
464                   &  / (scat_od_lw_liq+scat_od_lw_ice)
465            elsewhere
466              g_lw_cloud(:,jlev,jcol) = 0.0_jprb
467            end where
468            ssa_lw_cloud(:,jlev,jcol) = (scat_od_lw_liq + scat_od_lw_ice) &
469                 &                    / (od_lw_liq + od_lw_ice)
470          else
471            ! If longwave scattering is to be neglected then the
472            ! best approximation is to set the optical depth equal
473            ! to the absorption optical depth
474            od_lw_cloud(:,jlev,jcol) = od_lw_liq - scat_od_lw_liq &
475                 &                   + od_lw_ice - scat_od_lw_ice
476          end if
477          od_sw_cloud(:,jlev,jcol) = od_sw_liq + od_sw_ice
478          g_sw_cloud(:,jlev,jcol) = (g_sw_liq * scat_od_sw_liq &
479               &  + g_sw_ice * scat_od_sw_ice) &
480               &  / (scat_od_sw_liq + scat_od_sw_ice)
481          ssa_sw_cloud(:,jlev,jcol) &
482               &  = (scat_od_sw_liq + scat_od_sw_ice) / (od_sw_liq + od_sw_ice)
483        end if ! Cloud present
484      end do ! Loop over column
485    end do ! Loop over level
486
487    if (lhook) call dr_hook('radiation_cloud_optics:cloud_optics',1,hook_handle)
488
489  end subroutine cloud_optics
490
491end module radiation_cloud_optics
Note: See TracBrowser for help on using the repository browser.