source: LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_cloud_optics.F90 @ 5423

Last change on this file since 5423 was 4773, checked in by idelkadi, 13 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: 23.5 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
22  public
23
24contains
25
26  ! Provides elemental function "delta_eddington_scat_od"
27#include "radiation_delta_eddington.h"
28
29  !---------------------------------------------------------------------
30  ! Load cloud scattering data; this subroutine delegates to one
31  ! in radiation_cloud_optics_data.F90, but checks the size of
32  ! what is returned
33  subroutine setup_cloud_optics(config)
34
35    use parkind1,         only : jprb
36    use yomhook,          only : lhook, dr_hook, jphook
37
38    use radiation_io,     only : nulerr, radiation_abort
39    use radiation_config, only : config_type, IIceModelFu, IIceModelBaran, &
40         &                       IIceModelBaran2016, IIceModelBaran2017, &
41         &                       IIceModelYi, &
42         &                       ILiquidModelSOCRATES, ILiquidModelSlingo
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(jphook) :: 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, jphook
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_ice_optics_fu, only     : calc_ice_optics_fu_sw, &
217         &                                  calc_ice_optics_fu_lw
218    use radiation_ice_optics_baran, only  : calc_ice_optics_baran, &
219         &                                  calc_ice_optics_baran2016
220    use radiation_ice_optics_baran2017, only  : calc_ice_optics_baran2017
221    use radiation_ice_optics_yi, only     : calc_ice_optics_yi_sw, &
222         &                                  calc_ice_optics_yi_lw
223    use radiation_liquid_optics_socrates, only:calc_liq_optics_socrates
224    use radiation_liquid_optics_slingo, only:calc_liq_optics_slingo, &
225         &                                   calc_liq_optics_lindner_li
226
227    integer, intent(in) :: nlev               ! number of model levels
228    integer, intent(in) :: istartcol, iendcol ! range of columns to process
229    type(config_type), intent(in), target :: config
230    type(thermodynamics_type),intent(in)  :: thermodynamics
231    type(cloud_type),   intent(in)        :: cloud
232
233    ! Layer optical depth, single scattering albedo and g factor of
234    ! clouds in each longwave band, where the latter two
235    ! variables are only defined if cloud longwave scattering is
236    ! enabled (otherwise both are treated as zero).
237    real(jprb), dimension(config%n_bands_lw,nlev,istartcol:iendcol), intent(out) :: &
238         &   od_lw_cloud
239    real(jprb), dimension(config%n_bands_lw_if_scattering,nlev,istartcol:iendcol), &
240         &   intent(out) :: ssa_lw_cloud, g_lw_cloud
241
242    ! Layer optical depth, single scattering albedo and g factor of
243    ! clouds in each shortwave band
244    real(jprb), dimension(config%n_bands_sw,nlev,istartcol:iendcol), intent(out) :: &
245         &   od_sw_cloud, ssa_sw_cloud, g_sw_cloud
246
247    ! Longwave and shortwave optical depth, scattering optical depth
248    ! and asymmetry factor, for liquid and ice in all bands but a
249    ! single cloud layer
250    real(jprb), dimension(config%n_bands_lw) :: &
251         &  od_lw_liq, scat_od_lw_liq, g_lw_liq, &
252         &  od_lw_ice, scat_od_lw_ice, g_lw_ice
253    real(jprb), dimension(config%n_bands_sw) :: &
254         &  od_sw_liq, scat_od_sw_liq, g_sw_liq, &
255         &  od_sw_ice, scat_od_sw_ice, g_sw_ice
256
257    ! In-cloud water path of cloud liquid or ice (i.e. liquid or ice
258    ! gridbox-mean water path divided by cloud fraction); kg m-2
259    real(jprb) :: lwp_in_cloud, iwp_in_cloud
260
261    ! Full-level temperature (K)
262    real(jprb) :: temperature
263
264    ! Factor to convert gridbox-mean mixing ratio to in-cloud water
265    ! path
266    real(jprb) :: factor
267
268    integer    :: jcol, jlev, jb
269
270    real(jphook) :: hook_handle
271
272    if (lhook) call dr_hook('radiation_cloud_optics:cloud_optics',0,hook_handle)
273
274    if (config%iverbose >= 2) then
275      write(nulout,'(a)') 'Computing cloud absorption/scattering properties'
276    end if
277
278    associate(ho => config%cloud_optics)
279
280      ! Array-wise assignment
281      od_lw_cloud  = 0.0_jprb
282      od_sw_cloud  = 0.0_jprb
283      ssa_sw_cloud = 0.0_jprb
284      g_sw_cloud   = 0.0_jprb
285      if (config%do_lw_cloud_scattering) then
286        ssa_lw_cloud = 0.0_jprb
287        g_lw_cloud   = 0.0_jprb
288      end if
289
290      do jlev = 1,nlev
291        do jcol = istartcol,iendcol
292          ! Only do anything if cloud is present (assume that
293          ! cloud%crop_cloud_fraction has already been called)
294          if (cloud%fraction(jcol,jlev) > 0.0_jprb) then
295
296            ! Compute in-cloud liquid and ice water path
297            if (config%is_homogeneous) then
298              ! Homogeneous solvers assume cloud fills the box
299              ! horizontally, so we don't divide by cloud fraction
300              factor = ( thermodynamics%pressure_hl(jcol,jlev+1)    &
301                  &  -thermodynamics%pressure_hl(jcol,jlev  )  ) &
302                  &  / AccelDueToGravity
303            else
304              factor = ( thermodynamics%pressure_hl(jcol,jlev+1)    &
305                  &  -thermodynamics%pressure_hl(jcol,jlev  )  ) &
306                  &  / (AccelDueToGravity * cloud%fraction(jcol,jlev))
307            end if
308            lwp_in_cloud = factor * cloud%q_liq(jcol,jlev)
309            iwp_in_cloud = factor * cloud%q_ice(jcol,jlev)
310
311            ! Only compute liquid properties if liquid cloud is
312            ! present
313            if (lwp_in_cloud > 0.0_jprb) then
314              if (config%i_liq_model == ILiquidModelSOCRATES) then
315                ! Compute longwave properties
316                call calc_liq_optics_socrates(config%n_bands_lw, &
317                    &  config%cloud_optics%liq_coeff_lw, &
318                    &  lwp_in_cloud, cloud%re_liq(jcol,jlev), &
319                    &  od_lw_liq, scat_od_lw_liq, g_lw_liq)
320                ! Compute shortwave properties
321                call calc_liq_optics_socrates(config%n_bands_sw, &
322                    &  config%cloud_optics%liq_coeff_sw, &
323                    &  lwp_in_cloud, cloud%re_liq(jcol,jlev), &
324                    &  od_sw_liq, scat_od_sw_liq, g_sw_liq)
325              else if (config%i_liq_model == ILiquidModelSlingo) then
326                ! Compute longwave properties
327                call calc_liq_optics_lindner_li(config%n_bands_lw, &
328                    &  config%cloud_optics%liq_coeff_lw, &
329                    &  lwp_in_cloud, cloud%re_liq(jcol,jlev), &
330                    &  od_lw_liq, scat_od_lw_liq, g_lw_liq)
331                ! Compute shortwave properties
332                call calc_liq_optics_slingo(config%n_bands_sw, &
333                    &  config%cloud_optics%liq_coeff_sw, &
334                    &  lwp_in_cloud, cloud%re_liq(jcol,jlev), &
335                    &  od_sw_liq, scat_od_sw_liq, g_sw_liq)
336              else
337                write(nulerr,*) '*** Error: Unknown liquid model with code', &
338                    &          config%i_liq_model
339                call radiation_abort()
340              end if
341
342              ! Delta-Eddington scaling in the shortwave only
343              if (.not. config%do_sw_delta_scaling_with_gases) then
344                call delta_eddington_scat_od(od_sw_liq, scat_od_sw_liq, g_sw_liq)
345              end if
346              ! Originally delta-Eddington has been off in ecRad for
347              ! liquid clouds in the longwave, but it should be on
348              !call delta_eddington_scat_od(od_lw_liq, scat_od_lw_liq, g_lw_liq)
349
350            else
351              ! Liquid not present: set properties to zero
352              od_lw_liq = 0.0_jprb
353              scat_od_lw_liq = 0.0_jprb
354              g_lw_liq = 0.0_jprb
355
356              od_sw_liq = 0.0_jprb
357              scat_od_sw_liq = 0.0_jprb
358              g_sw_liq = 0.0_jprb
359            end if ! Liquid present
360
361            ! Only compute ice properties if ice cloud is present
362            if (iwp_in_cloud > 0.0_jprb) then
363              if (config%i_ice_model == IIceModelBaran) then
364                ! Compute longwave properties
365                call calc_ice_optics_baran(config%n_bands_lw, &
366                    &  config%cloud_optics%ice_coeff_lw, &
367                    &  iwp_in_cloud, cloud%q_ice(jcol,jlev), &
368                    &  od_lw_ice, scat_od_lw_ice, g_lw_ice)
369                ! Compute shortwave properties
370                call calc_ice_optics_baran(config%n_bands_sw, &
371                    &  config%cloud_optics%ice_coeff_sw, &
372                    &  iwp_in_cloud, cloud%q_ice(jcol,jlev), &
373                    &  od_sw_ice, scat_od_sw_ice, g_sw_ice)
374              else if (config%i_ice_model == IIceModelBaran2016) then
375                temperature = 0.5_jprb * (thermodynamics%temperature_hl(jcol,jlev) &
376                    &                   +thermodynamics%temperature_hl(jcol,jlev+1))
377                ! Compute longwave properties
378                call calc_ice_optics_baran2016(config%n_bands_lw, &
379                    &  config%cloud_optics%ice_coeff_lw, &
380                    &  iwp_in_cloud, cloud%q_ice(jcol,jlev), &
381                    &  temperature, &
382                    &  od_lw_ice, scat_od_lw_ice, g_lw_ice)
383                ! Compute shortwave properties
384                call calc_ice_optics_baran2016(config%n_bands_sw, &
385                    &  config%cloud_optics%ice_coeff_sw, &
386                    &  iwp_in_cloud, cloud%q_ice(jcol,jlev), &
387                    &  temperature, &
388                    &  od_sw_ice, scat_od_sw_ice, g_sw_ice)
389              else if (config%i_ice_model == IIceModelBaran2017) then
390                temperature = 0.5_jprb * (thermodynamics%temperature_hl(jcol,jlev) &
391                    &                   +thermodynamics%temperature_hl(jcol,jlev+1))
392                ! Compute longwave properties
393                call calc_ice_optics_baran2017(config%n_bands_lw, &
394                    &  config%cloud_optics%ice_coeff_gen, &
395                    &  config%cloud_optics%ice_coeff_lw, &
396                    &  iwp_in_cloud, cloud%q_ice(jcol,jlev), &
397                    &  temperature, &
398                    &  od_lw_ice, scat_od_lw_ice, g_lw_ice)
399                ! Compute shortwave properties
400                call calc_ice_optics_baran2017(config%n_bands_sw, &
401                    &  config%cloud_optics%ice_coeff_gen, &
402                    &  config%cloud_optics%ice_coeff_sw, &
403                    &  iwp_in_cloud, cloud%q_ice(jcol,jlev), &
404                    &  temperature, &
405                    &  od_sw_ice, scat_od_sw_ice, g_sw_ice)
406              else if (config%i_ice_model == IIceModelFu) then
407                ! Compute longwave properties
408                call calc_ice_optics_fu_lw(config%n_bands_lw, &
409                    &  config%cloud_optics%ice_coeff_lw, &
410                    &  iwp_in_cloud, cloud%re_ice(jcol,jlev), &
411                    &  od_lw_ice, scat_od_lw_ice, g_lw_ice)
412                if (config%do_fu_lw_ice_optics_bug) then
413                  ! Reproduce bug in old IFS scheme
414                  scat_od_lw_ice = od_lw_ice - scat_od_lw_ice
415                end if
416                ! Compute shortwave properties
417                call calc_ice_optics_fu_sw(config%n_bands_sw, &
418                    &  config%cloud_optics%ice_coeff_sw, &
419                    &  iwp_in_cloud, cloud%re_ice(jcol,jlev), &
420                    &  od_sw_ice, scat_od_sw_ice, g_sw_ice)
421              else if (config%i_ice_model == IIceModelYi) then
422                ! Compute longwave properties
423                call calc_ice_optics_yi_lw(config%n_bands_lw, &
424                    &  config%cloud_optics%ice_coeff_lw, &
425                    &  iwp_in_cloud, cloud%re_ice(jcol,jlev), &
426                    &  od_lw_ice, scat_od_lw_ice, g_lw_ice)
427                ! Compute shortwave properties
428                call calc_ice_optics_yi_sw(config%n_bands_sw, &
429                    &  config%cloud_optics%ice_coeff_sw, &
430                    &  iwp_in_cloud, cloud%re_ice(jcol,jlev), &
431                    &  od_sw_ice, scat_od_sw_ice, g_sw_ice)
432              else
433                write(nulerr,*) '*** Error: Unknown ice model with code', &
434                    &          config%i_ice_model
435                call radiation_abort()
436              end if
437
438              ! Delta-Eddington scaling in both longwave and shortwave
439              ! (assume that particles are larger than wavelength even
440              ! in longwave)
441              if (.not. config%do_sw_delta_scaling_with_gases) then
442                call delta_eddington_scat_od(od_sw_ice, scat_od_sw_ice, g_sw_ice)
443              end if
444              call delta_eddington_scat_od(od_lw_ice, scat_od_lw_ice, g_lw_ice)
445
446            else
447              ! Ice not present: set properties to zero
448              od_lw_ice = 0.0_jprb
449              scat_od_lw_ice = 0.0_jprb
450              g_lw_ice = 0.0_jprb
451
452              od_sw_ice = 0.0_jprb
453              scat_od_sw_ice = 0.0_jprb
454              g_sw_ice = 0.0_jprb
455            end if ! Ice present
456
457            ! Combine liquid and ice
458            if (config%do_lw_cloud_scattering) then
459  ! Added for DWD (2020)
460  !NEC$ shortloop
461              do jb = 1, config%n_bands_lw
462                od_lw_cloud(jb,jlev,jcol) = od_lw_liq(jb) + od_lw_ice(jb)
463                if (scat_od_lw_liq(jb)+scat_od_lw_ice(jb) > 0.0_jprb) then
464                  g_lw_cloud(jb,jlev,jcol) = (g_lw_liq(jb) * scat_od_lw_liq(jb) &
465                    &  + g_lw_ice(jb) * scat_od_lw_ice(jb)) &
466                    &  / (scat_od_lw_liq(jb)+scat_od_lw_ice(jb))
467                else
468                  g_lw_cloud(jb,jlev,jcol) = 0.0_jprb
469                end if
470                ssa_lw_cloud(jb,jlev,jcol) = (scat_od_lw_liq(jb) + scat_od_lw_ice(jb)) &
471                  &                    / (od_lw_liq(jb) + od_lw_ice(jb))
472              end do
473            else
474              ! If longwave scattering is to be neglected then the
475              ! best approximation is to set the optical depth equal
476              ! to the absorption optical depth
477  ! Added for DWD (2020)
478  !NEC$ shortloop
479              do jb = 1, config%n_bands_lw
480                od_lw_cloud(jb,jlev,jcol) = od_lw_liq(jb) - scat_od_lw_liq(jb) &
481                      &                   + od_lw_ice(jb) - scat_od_lw_ice(jb)
482              end do
483            end if
484  ! Added for DWD (2020)
485  !NEC$ shortloop
486            do jb = 1, config%n_bands_sw
487              od_sw_cloud(jb,jlev,jcol) = od_sw_liq(jb) + od_sw_ice(jb)
488              g_sw_cloud(jb,jlev,jcol) = (g_sw_liq(jb) * scat_od_sw_liq(jb) &
489                &  + g_sw_ice(jb) * scat_od_sw_ice(jb)) &
490                &  / (scat_od_sw_liq(jb) + scat_od_sw_ice(jb))
491              ssa_sw_cloud(jb,jlev,jcol) &
492                &  = (scat_od_sw_liq(jb) + scat_od_sw_ice(jb)) / (od_sw_liq(jb) + od_sw_ice(jb))
493            end do
494          end if ! Cloud present
495        end do ! Loop over column
496      end do ! Loop over level
497
498    end associate
499
500    if (lhook) call dr_hook('radiation_cloud_optics:cloud_optics',1,hook_handle)
501
502  end subroutine cloud_optics
503
504end module radiation_cloud_optics
Note: See TracBrowser for help on using the repository browser.