source: LMDZ6/trunk/libf/phylmd/ecrad/radiation_aerosol_optics.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: 26.7 KB
Line 
1! radiation_aerosol_optics.F90 - Computing aerosol optical properties
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!   2018-04-15  R Hogan  Add "direct" option
17
18module radiation_aerosol_optics
19
20  implicit none
21  public
22
23contains
24
25  ! Provides the elemental function "delta_eddington_extensive"
26#include "radiation_delta_eddington.h"
27
28  !---------------------------------------------------------------------
29  ! Load aerosol scattering data; this subroutine delegates to one
30  ! in radiation_aerosol_optics_data.F90
31  subroutine setup_aerosol_optics(config)
32
33    use parkind1,                      only : jprb
34    use yomhook,                       only : lhook, dr_hook
35    use radiation_config,              only : config_type
36    use radiation_aerosol_optics_data, only : aerosol_optics_type
37    use radiation_io,                  only : nulerr, radiation_abort
38
39    type(config_type), intent(inout) :: config
40
41    real(jprb) :: hook_handle
42
43    if (lhook) call dr_hook('radiation_aerosol_optics:setup_aerosol_optics',0,hook_handle)
44
45    if (config%n_aerosol_types > 0) then
46      ! Load data from file and prepare to map config%n_aerosol_types
47      ! aerosol types
48      call config%aerosol_optics%setup(trim(config%aerosol_optics_file_name), &
49           &                        config%n_aerosol_types, iverbose=config%iverbosesetup)
50
51      ! Check agreement in number of bands
52      if (config%n_bands_lw /= config%aerosol_optics%n_bands_lw) then
53        write(nulerr,'(a)') '*** Error: number of longwave bands does not match aerosol optics look-up table'
54        call radiation_abort()
55      end if
56      if (config%n_bands_sw /= config%aerosol_optics%n_bands_sw) then
57        write(nulerr,'(a)') '*** Error: number of shortwave bands does not match aerosol optics look-up table'
58        call radiation_abort()
59      end if
60
61      ! Map aerosol types to those loaded from the data file
62      call config%aerosol_optics%set_types(config%i_aerosol_type_map(1:config%n_aerosol_types))
63    end if
64
65    call config%aerosol_optics%print_description(config%i_aerosol_type_map(1:config%n_aerosol_types))
66
67    if (lhook) call dr_hook('radiation_aerosol_optics:setup_aerosol_optics',1,hook_handle)
68
69  end subroutine setup_aerosol_optics
70
71
72  !---------------------------------------------------------------------
73  ! Compute aerosol optical properties and add to existing gas optical
74  ! depth and scattering properties
75  subroutine add_aerosol_optics(nlev,istartcol,iendcol, &
76       &  config, thermodynamics, gas, aerosol, &
77       &  od_lw, ssa_lw, g_lw, od_sw, ssa_sw, g_sw)
78
79    use parkind1,                      only : jprb
80    use radiation_io,                  only : nulout, nulerr, radiation_abort
81    use yomhook,                       only : lhook, dr_hook
82    use radiation_config,              only : config_type
83    use radiation_thermodynamics,      only : thermodynamics_type
84    use radiation_gas,                 only : gas_type, IH2O, IMassMixingRatio
85    use radiation_aerosol,             only : aerosol_type
86    use radiation_constants,           only : AccelDueToGravity
87    use radiation_aerosol_optics_data, only : aerosol_optics_type, &
88         &  IAerosolClassUndefined,   IAerosolClassIgnored, &
89         &  IAerosolClassHydrophobic, IAerosolClassHydrophilic
90
91    integer, intent(in) :: nlev               ! number of model levels
92    integer, intent(in) :: istartcol, iendcol ! range of columns to process
93    type(config_type), intent(in), target :: config
94    type(thermodynamics_type),intent(in)  :: thermodynamics
95    type(gas_type),           intent(in)  :: gas
96    type(aerosol_type),       intent(in)  :: aerosol
97    ! Optical depth, single scattering albedo and asymmetry factor of
98    ! the atmosphere (gases on input, gases and aerosols on output)
99    ! for each g point. Note that longwave ssa and asymmetry and
100    ! shortwave asymmetry are all zero for gases, so are not yet
101    ! defined on input and are therefore intent(out).
102    real(jprb), dimension(config%n_g_lw,nlev,istartcol:iendcol), &
103         &   intent(inout) :: od_lw
104    real(jprb), dimension(config%n_g_lw_if_scattering,nlev,istartcol:iendcol), &
105         &   intent(out)   :: ssa_lw, g_lw
106    real(jprb), dimension(config%n_g_sw,nlev,istartcol:iendcol), &
107         &   intent(inout) :: od_sw, ssa_sw
108    real(jprb), dimension(config%n_g_sw,nlev,istartcol:iendcol), &
109         &   intent(out)   :: g_sw
110
111    ! Extinction optical depth, scattering optical depth and
112    ! asymmetry-times-scattering-optical-depth for all the aerosols at
113    ! a point in space for each spectral band of the shortwave and
114    ! longwave spectrum
115    real(jprb), dimension(config%n_bands_sw) &
116         & :: od_sw_aerosol, scat_sw_aerosol, scat_g_sw_aerosol, local_od_sw
117    real(jprb), dimension(config%n_bands_lw) :: od_lw_aerosol, local_od_lw
118    real(jprb), dimension(config%n_bands_lw_if_scattering) &
119         & :: scat_lw_aerosol, scat_g_lw_aerosol
120
121    real(jprb) :: h2o_mmr(istartcol:iendcol,nlev)
122
123    real(jprb) :: rh ! Relative humidity with respect to liquid water
124
125    ! Factor (kg m-2) to convert mixing ratio (kg kg-1) to mass in
126    ! path (kg m-2)
127    real(jprb) :: factor
128
129    ! Temporary extinction and scattering optical depths of aerosol
130    ! plus gas
131    real(jprb) :: local_od, local_scat
132
133    ! Loop indices for column, level, g point, band and aerosol type
134    integer :: jcol, jlev, jg, jtype
135
136    ! Range of levels over which aerosols are present
137    integer :: istartlev, iendlev
138
139    ! Indices to spectral band and relative humidity look-up table
140    integer :: iband, irh
141
142    ! Pointer to the aerosol optics coefficients for brevity of access
143    type(aerosol_optics_type), pointer :: ao
144
145    real(jprb) :: hook_handle
146
147    if (lhook) call dr_hook('radiation_aerosol_optics:add_aerosol_optics',0,hook_handle)
148
149    if (aerosol%is_direct) then
150      ! Aerosol optical properties have been provided in each band
151      ! directly by the user
152      call add_aerosol_optics_direct(nlev,istartcol,iendcol, &
153           &  config, aerosol, &
154           &  od_lw, ssa_lw, g_lw, od_sw, ssa_sw, g_sw)
155    else
156      ! Aerosol mixing ratios have been provided
157
158      do jtype = 1,config%n_aerosol_types
159        if (config%aerosol_optics%iclass(jtype) == IAerosolClassUndefined) then
160          write(nulerr,'(a)') '*** Error: not all aerosol types are defined'
161          call radiation_abort()
162        end if
163      end do
164
165      if (config%iverbose >= 2) then
166        write(nulout,'(a)') 'Computing aerosol absorption/scattering properties'
167      end if
168
169      ao => config%aerosol_optics
170
171      istartlev = lbound(aerosol%mixing_ratio,2)
172      iendlev   = ubound(aerosol%mixing_ratio,2)
173
174      if (ubound(aerosol%mixing_ratio,3) /= config%n_aerosol_types) then
175        write(nulerr,'(a,i0,a,i0)') '*** Error: aerosol%mixing_ratio contains ', &
176             &  ubound(aerosol%mixing_ratio,3), ' aerosol types, expected ', &
177             &  config%n_aerosol_types
178        call radiation_abort()
179      end if
180
181      ! Set variables to zero that may not have been previously
182      g_sw = 0.0_jprb
183      if (config%do_lw_aerosol_scattering) then
184        ssa_lw = 0.0_jprb
185        g_lw   = 0.0_jprb
186      end if
187
188      call gas%get(IH2O, IMassMixingRatio, h2o_mmr, istartcol=istartcol)
189
190      ! Loop over position
191      do jlev = istartlev,iendlev
192        do jcol = istartcol,iendcol
193          ! Compute relative humidity with respect to liquid
194          ! saturation and the index to the relative-humidity index of
195          ! hydrophilic-aerosol data
196          rh  = h2o_mmr(jcol,jlev) / thermodynamics%h2o_sat_liq(jcol,jlev)
197          irh = ao%calc_rh_index(rh)
198
199          factor = ( thermodynamics%pressure_hl(jcol,jlev+1) &
200               &    -thermodynamics%pressure_hl(jcol,jlev  )  ) &
201               &   / AccelDueToGravity 
202
203          ! Reset temporary arrays
204          od_sw_aerosol     = 0.0_jprb
205          scat_sw_aerosol   = 0.0_jprb
206          scat_g_sw_aerosol = 0.0_jprb
207          od_lw_aerosol     = 0.0_jprb
208          scat_lw_aerosol   = 0.0_jprb
209          scat_g_lw_aerosol = 0.0_jprb
210
211          do jtype = 1,config%n_aerosol_types
212            ! Add the optical depth, scattering optical depth and
213            ! scattering optical depth-weighted asymmetry factor for
214            ! this aerosol type to the total for all aerosols.  Note
215            ! that the following expressions are array-wise, the
216            ! dimension being spectral band.
217            if (ao%iclass(jtype) == IAerosolClassHydrophobic) then
218              local_od_sw = factor * aerosol%mixing_ratio(jcol,jlev,jtype) &
219                   &  * ao%mass_ext_sw_phobic(:,ao%itype(jtype))
220              od_sw_aerosol = od_sw_aerosol + local_od_sw
221              scat_sw_aerosol = scat_sw_aerosol &
222                   &  + local_od_sw * ao%ssa_sw_phobic(:,ao%itype(jtype))
223              scat_g_sw_aerosol = scat_g_sw_aerosol &
224                   &  + local_od_sw * ao%ssa_sw_phobic(:,ao%itype(jtype)) &
225                   &  * ao%g_sw_phobic(:,ao%itype(jtype))
226              if (config%do_lw_aerosol_scattering) then
227                local_od_lw = factor * aerosol%mixing_ratio(jcol,jlev,jtype) &
228                     &  * ao%mass_ext_lw_phobic(:,ao%itype(jtype))
229                od_lw_aerosol = od_lw_aerosol + local_od_lw
230                scat_lw_aerosol = scat_lw_aerosol &
231                     &  + local_od_lw * ao%ssa_lw_phobic(:,ao%itype(jtype))
232                scat_g_lw_aerosol = scat_g_lw_aerosol &
233                     &  + local_od_lw * ao%ssa_lw_phobic(:,ao%itype(jtype)) &
234                     &  * ao%g_lw_phobic(:,ao%itype(jtype))
235              else
236                ! If aerosol longwave scattering is not included then we
237                ! weight the optical depth by the single scattering
238                ! co-albedo
239                od_lw_aerosol = od_lw_aerosol &
240                     &  + factor * aerosol%mixing_ratio(jcol,jlev,jtype) &
241                     &  * ao%mass_ext_lw_phobic(:,ao%itype(jtype)) &
242                     &  * (1.0_jprb - ao%ssa_lw_phobic(:,ao%itype(jtype)))
243              end if
244            else if (ao%iclass(jtype) == IAerosolClassHydrophilic) then
245              ! Hydrophilic aerosols require the look-up tables to
246              ! be indexed with irh
247              local_od_sw = factor * aerosol%mixing_ratio(jcol,jlev,jtype) &
248                   &  * ao%mass_ext_sw_philic(:,irh,ao%itype(jtype))
249              od_sw_aerosol = od_sw_aerosol + local_od_sw
250              scat_sw_aerosol = scat_sw_aerosol &
251                   &  + local_od_sw * ao%ssa_sw_philic(:,irh,ao%itype(jtype))
252              scat_g_sw_aerosol = scat_g_sw_aerosol &
253                   &  + local_od_sw * ao%ssa_sw_philic(:,irh,ao%itype(jtype)) &
254                   &  * ao%g_sw_philic(:,irh,ao%itype(jtype))
255              if (config%do_lw_aerosol_scattering) then
256                local_od_lw = factor * aerosol%mixing_ratio(jcol,jlev,jtype) &
257                     &  * ao%mass_ext_lw_philic(:,irh,ao%itype(jtype))
258                od_lw_aerosol = od_lw_aerosol + local_od_lw
259                scat_lw_aerosol = scat_lw_aerosol &
260                     &  + local_od_lw * ao%ssa_lw_philic(:,irh,ao%itype(jtype))
261                scat_g_lw_aerosol = scat_g_lw_aerosol &
262                     &  + local_od_lw * ao%ssa_lw_philic(:,irh,ao%itype(jtype)) &
263                     &  * ao%g_lw_philic(:,irh,ao%itype(jtype))
264              else
265                ! If aerosol longwave scattering is not included then we
266                ! weight the optical depth by the single scattering
267                ! co-albedo
268                od_lw_aerosol = od_lw_aerosol &
269                     &  + factor * aerosol%mixing_ratio(jcol,jlev,jtype) &
270                     &  * ao%mass_ext_lw_philic(:,irh,ao%itype(jtype)) &
271                     &  * (1.0_jprb - ao%ssa_lw_philic(:,irh,ao%itype(jtype)))
272              end if
273            end if
274            ! Implicitly, if ao%iclass(jtype) == IAerosolClassNone, then
275            ! no aerosol scattering properties are added
276
277          end do ! Loop over aerosol type
278
279          if (.not. config%do_sw_delta_scaling_with_gases) then
280            ! Delta-Eddington scaling on aerosol only.  Note that if
281            ! do_sw_delta_scaling_with_gases==.true. then the delta
282            ! scaling is done to the cloud-aerosol-gas mixture inside
283            ! the solver
284            call delta_eddington_extensive(od_sw_aerosol, scat_sw_aerosol, &
285                 &                         scat_g_sw_aerosol)
286          end if
287
288          ! Combine aerosol shortwave scattering properties with gas
289          ! properties (noting that any gas scattering will have an
290          ! asymmetry factor of zero)
291          if (od_sw_aerosol(1) > 0.0_jprb) then
292            do jg = 1,config%n_g_sw
293              iband = config%i_band_from_reordered_g_sw(jg)
294              local_od = od_sw(jg,jlev,jcol) + od_sw_aerosol(iband)
295              local_scat = ssa_sw(jg,jlev,jcol) * od_sw(jg,jlev,jcol) &
296                   &  + scat_sw_aerosol(iband)
297              ! Note that asymmetry_sw of gases is zero so the following
298              ! simply weights the aerosol asymmetry by the scattering
299              ! optical depth
300              g_sw(jg,jlev,jcol) = scat_g_sw_aerosol(iband) / local_scat
301              ssa_sw(jg,jlev,jcol) = local_scat / local_od
302              od_sw (jg,jlev,jcol) = local_od
303            end do
304          end if
305
306          ! Combine aerosol longwave scattering properties with gas
307          ! properties, noting that in the longwave, gases do not
308          ! scatter at all
309          if (config%do_lw_aerosol_scattering) then
310
311            call delta_eddington_extensive(od_lw_aerosol, scat_lw_aerosol, &
312                 &                         scat_g_lw_aerosol) 
313
314            do jg = 1,config%n_g_lw
315              iband = config%i_band_from_reordered_g_lw(jg)
316              if (od_lw_aerosol(iband) > 0.0_jprb) then
317                ! All scattering is due to aerosols, therefore the
318                ! asymmetry factor is equal to the value for aerosols
319                if (scat_lw_aerosol(iband) > 0.0_jprb) then
320                  g_lw(jg,jlev,jcol) = scat_g_lw_aerosol(iband) &
321                       &  / scat_lw_aerosol(iband)
322                else
323                  g_lw(jg,jlev,jcol) = 0.0_jprb
324                end if
325                local_od = od_lw(jg,jlev,jcol) + od_lw_aerosol(iband)
326                ssa_lw(jg,jlev,jcol) = scat_lw_aerosol(iband) / local_od
327                od_lw (jg,jlev,jcol) = local_od
328              end if
329            end do
330          else
331            do jg = 1,config%n_g_lw
332              od_lw(jg,jlev,jcol) = od_lw(jg,jlev,jcol) &
333                   &  + od_lw_aerosol(config%i_band_from_reordered_g_lw(jg))
334            end do
335          end if
336
337        end do ! Loop over column
338      end do ! Loop over level
339
340    end if
341
342    if (lhook) call dr_hook('radiation_aerosol_optics:add_aerosol_optics',1,hook_handle)
343
344  end subroutine add_aerosol_optics
345
346
347  !---------------------------------------------------------------------
348  ! Add precomputed optical properties to gas optical depth and
349  ! scattering properties
350  subroutine add_aerosol_optics_direct(nlev,istartcol,iendcol, &
351       &  config, aerosol, &
352       &  od_lw, ssa_lw, g_lw, od_sw, ssa_sw, g_sw)
353
354    use parkind1,                      only : jprb
355    use radiation_io,                  only : nulerr, radiation_abort
356    use yomhook,                       only : lhook, dr_hook
357    use radiation_config,              only : config_type
358    use radiation_aerosol,             only : aerosol_type
359
360    integer, intent(in) :: nlev               ! number of model levels
361    integer, intent(in) :: istartcol, iendcol ! range of columns to process
362    type(config_type), intent(in), target :: config
363    type(aerosol_type),       intent(in)  :: aerosol
364    ! Optical depth, single scattering albedo and asymmetry factor of
365    ! the atmosphere (gases on input, gases and aerosols on output)
366    ! for each g point. Note that longwave ssa and asymmetry and
367    ! shortwave asymmetry are all zero for gases, so are not yet
368    ! defined on input and are therefore intent(out).
369    real(jprb), dimension(config%n_g_lw,nlev,istartcol:iendcol), &
370         &   intent(inout) :: od_lw
371    real(jprb), dimension(config%n_g_lw_if_scattering,nlev,istartcol:iendcol), &
372         &   intent(out)   :: ssa_lw, g_lw
373    real(jprb), dimension(config%n_g_sw,nlev,istartcol:iendcol), &
374         &   intent(inout) :: od_sw, ssa_sw
375    real(jprb), dimension(config%n_g_sw,nlev,istartcol:iendcol), &
376         &   intent(out)   :: g_sw
377
378    ! Temporary extinction and scattering optical depths of aerosol
379    ! plus gas
380    real(jprb) :: local_od, local_scat
381
382    ! Extinction optical depth, scattering optical depth and
383    ! asymmetry-times-scattering-optical-depth for all the aerosols at
384    ! a point in space for each spectral band of the shortwave and
385    ! longwave spectrum
386    real(jprb), dimension(config%n_bands_sw) &
387         & :: od_sw_aerosol, scat_sw_aerosol, scat_g_sw_aerosol
388    real(jprb), dimension(config%n_bands_lw) :: od_lw_aerosol
389    real(jprb), dimension(config%n_bands_lw_if_scattering) &
390         & :: scat_lw_aerosol, scat_g_lw_aerosol
391
392    ! Loop indices for column, level, g point and band
393    integer :: jcol, jlev, jg
394
395    ! Range of levels over which aerosols are present
396    integer :: istartlev, iendlev
397
398    ! Indices to spectral band
399    integer :: iband
400
401    real(jprb) :: hook_handle
402
403    if (lhook) call dr_hook('radiation_aerosol_optics:add_aerosol_optics_direct',0,hook_handle)
404
405    if (config%do_sw) then
406      ! Check array dimensions
407      if (ubound(aerosol%od_sw,1) /= config%n_bands_sw) then
408        write(nulerr,'(a,i0,a,i0)') '*** Error: aerosol%od_sw contains ', &
409           &  ubound(aerosol%od_sw,1), ' band, expected ', &
410           &  config%n_bands_sw
411        call radiation_abort()
412      end if
413
414      istartlev = lbound(aerosol%od_sw,2)
415      iendlev   = ubound(aerosol%od_sw,2)
416
417      ! Set variables to zero that may not have been previously
418      g_sw = 0.0_jprb
419
420      ! Loop over position
421      do jcol = istartcol,iendcol
422        do jlev = istartlev,iendlev
423          od_sw_aerosol = aerosol%od_sw(:,jlev,jcol)
424          scat_sw_aerosol = aerosol%ssa_sw(:,jlev,jcol) * od_sw_aerosol
425          scat_g_sw_aerosol = aerosol%g_sw(:,jlev,jcol) * scat_sw_aerosol
426
427          if (.not. config%do_sw_delta_scaling_with_gases) then
428            ! Delta-Eddington scaling on aerosol only.  Note that if
429            ! do_sw_delta_scaling_with_gases==.true. then the delta
430            ! scaling is done to the cloud-aerosol-gas mixture inside
431            ! the solver
432            call delta_eddington_extensive(od_sw_aerosol, scat_sw_aerosol, &
433                 &                         scat_g_sw_aerosol)
434          end if
435
436          ! Combine aerosol shortwave scattering properties with gas
437          ! properties (noting that any gas scattering will have an
438          ! asymmetry factor of zero)
439          if (od_sw_aerosol(1) > 0.0_jprb) then
440            do jg = 1,config%n_g_sw
441              iband = config%i_band_from_reordered_g_sw(jg)
442              local_od = od_sw(jg,jlev,jcol) + od_sw_aerosol(iband)
443              local_scat = ssa_sw(jg,jlev,jcol) * od_sw(jg,jlev,jcol) &
444                   &  + scat_sw_aerosol(iband)
445              ! Note that asymmetry_sw of gases is zero so the following
446              ! simply weights the aerosol asymmetry by the scattering
447              ! optical depth
448              g_sw(jg,jlev,jcol) = scat_g_sw_aerosol(iband) / local_scat
449              ssa_sw(jg,jlev,jcol) = local_scat / local_od
450              od_sw (jg,jlev,jcol) = local_od
451            end do
452          end if
453        end do
454      end do
455
456    end if
457
458    if (config%do_lw) then
459
460      if (ubound(aerosol%od_lw,1) /= config%n_bands_lw) then
461        write(nulerr,'(a,i0,a,i0)') '*** Error: aerosol%od_lw contains ', &
462           &  ubound(aerosol%od_lw,1), ' band, expected ', &
463           &  config%n_bands_lw
464        call radiation_abort()
465      end if
466
467      istartlev = lbound(aerosol%od_lw,2)
468      iendlev   = ubound(aerosol%od_lw,2)
469
470      if (config%do_lw_aerosol_scattering) then
471        ssa_lw = 0.0_jprb
472        g_lw   = 0.0_jprb
473 
474        ! Loop over position
475        do jcol = istartcol,iendcol
476          do jlev = istartlev,iendlev
477            od_lw_aerosol = aerosol%od_lw(:,jlev,jcol)
478            scat_lw_aerosol = aerosol%ssa_lw(:,jlev,jcol) * od_lw_aerosol
479            scat_g_lw_aerosol = aerosol%g_lw(:,jlev,jcol) * scat_lw_aerosol
480           
481            call delta_eddington_extensive(od_lw_aerosol, scat_lw_aerosol, &
482                 &                         scat_g_lw_aerosol)
483           
484            do jg = 1,config%n_g_lw
485              iband = config%i_band_from_reordered_g_lw(jg)
486              if (od_lw_aerosol(iband) > 0.0_jprb) then
487                ! All scattering is due to aerosols, therefore the
488                ! asymmetry factor is equal to the value for aerosols
489                if (scat_lw_aerosol(iband) > 0.0_jprb) then
490                  g_lw(jg,jlev,jcol) = scat_g_lw_aerosol(iband) &
491                       &  / scat_lw_aerosol(iband)
492                else
493                  g_lw(jg,jlev,jcol) = 0.0_jprb
494                end if
495                local_od = od_lw(jg,jlev,jcol) + od_lw_aerosol(iband)
496                ssa_lw(jg,jlev,jcol) = scat_lw_aerosol(iband) / local_od
497                od_lw (jg,jlev,jcol) = local_od
498              end if
499            end do
500          end do
501        end do
502
503      else ! No longwave scattering
504
505        ! Loop over position
506        do jcol = istartcol,iendcol
507          do jlev = istartlev,iendlev
508            ! If aerosol longwave scattering is not included then we
509            ! weight the optical depth by the single scattering
510            ! co-albedo
511            od_lw_aerosol = aerosol%od_lw(:,jlev,jcol) &
512                 &  * (1.0_jprb - aerosol%ssa_lw(:,jlev,jcol))
513            do jg = 1,config%n_g_lw
514              od_lw(jg,jlev,jcol) = od_lw(jg,jlev,jcol) &
515                   &  + od_lw_aerosol(config%i_band_from_reordered_g_lw(jg))
516            end do
517          end do
518        end do
519
520      end if
521    end if
522
523
524    if (lhook) call dr_hook('radiation_aerosol_optics:add_aerosol_optics_direct',1,hook_handle)
525
526  end subroutine add_aerosol_optics_direct
527 
528
529  !---------------------------------------------------------------------
530  ! Sometimes it is useful to specify aerosol in terms of its optical
531  ! depth at a particular wavelength.  This function returns the dry
532  ! shortwave mass-extinction coefficient, i.e. the extinction cross
533  ! section per unit mass, for aerosol of type "itype" at shortwave
534  ! band "iband". For hydrophilic types, the value at the first
535  ! relative humidity bin is taken.
536  function dry_aerosol_sw_mass_extinction(config, itype, iband)
537
538    use parkind1,                      only : jprb
539    use radiation_config,              only : config_type
540    use radiation_aerosol_optics_data, only : aerosol_optics_type, &
541         &  IAerosolClassUndefined,   IAerosolClassIgnored, &
542         &  IAerosolClassHydrophobic, IAerosolClassHydrophilic
543
544    type(config_type), intent(in), target :: config
545
546    ! Aerosol type and shortwave band as indices to the array
547    integer, intent(in) :: itype, iband
548   
549    real(jprb) dry_aerosol_sw_mass_extinction
550
551    ! Pointer to the aerosol optics coefficients for brevity of access
552    type(aerosol_optics_type), pointer :: ao
553
554    ao => config%aerosol_optics
555
556    if (ao%iclass(itype) == IAerosolClassHydrophobic) then
557      dry_aerosol_sw_mass_extinction = ao%mass_ext_sw_phobic(iband,ao%itype(itype))
558    else if (ao%iclass(itype) == IAerosolClassHydrophilic) then
559      ! Take the value at the first relative-humidity bin for the
560      ! "dry" aerosol value
561      dry_aerosol_sw_mass_extinction = ao%mass_ext_sw_philic(iband,1,ao%itype(itype))
562    else
563      dry_aerosol_sw_mass_extinction = 0.0_jprb
564    end if
565
566  end function dry_aerosol_sw_mass_extinction
567
568
569  !---------------------------------------------------------------------
570  ! Compute aerosol extinction coefficient at a particular shortwave
571  ! band and a single height - this is useful for visibility
572  ! diagnostics
573  subroutine aerosol_sw_extinction(ncol,istartcol,iendcol, &
574       &  config, iband, mixing_ratio, relative_humidity, extinction)
575
576    use parkind1,                      only : jprb
577    use yomhook,                       only : lhook, dr_hook
578    use radiation_io,                  only : nulerr, radiation_abort
579    use radiation_config,              only : config_type
580    use radiation_aerosol_optics_data, only : aerosol_optics_type, &
581         &  IAerosolClassUndefined,   IAerosolClassIgnored, &
582         &  IAerosolClassHydrophobic, IAerosolClassHydrophilic
583
584    integer, intent(in) :: ncol               ! number of columns
585    integer, intent(in) :: istartcol, iendcol ! range of columns to process
586    type(config_type), intent(in), target :: config
587    integer, intent(in)     :: iband ! Index of required spectral band
588    real(jprb), intent(in)  :: mixing_ratio(ncol,config%n_aerosol_types)
589    real(jprb), intent(in)  :: relative_humidity(ncol)
590    real(jprb), intent(out) :: extinction(ncol)
591
592    ! Local aerosol extinction
593    real(jprb) :: ext
594
595    ! Pointer to the aerosol optics coefficients for brevity of access
596    type(aerosol_optics_type), pointer :: ao
597   
598    ! Loop indices for column and aerosol type
599    integer :: jcol, jtype
600
601    ! Relative humidity index
602    integer :: irh
603
604    real(jprb) :: hook_handle
605
606    if (lhook) call dr_hook('radiation_aerosol_optics:aerosol_sw_extinction',0,hook_handle)
607
608    do jtype = 1,config%n_aerosol_types
609      if (config%aerosol_optics%iclass(jtype) == IAerosolClassUndefined) then
610        write(nulerr,'(a)') '*** Error: not all aerosol types are defined'
611        call radiation_abort()
612      end if
613    end do
614
615    ao => config%aerosol_optics
616
617    ! Loop over position
618    do jcol = istartcol,iendcol
619      ext = 0.0_jprb
620      ! Get relative-humidity index
621      irh = ao%calc_rh_index(relative_humidity(jcol))
622      ! Add extinction coefficients from each aerosol type
623      do jtype = 1,config%n_aerosol_types
624        if (ao%iclass(jtype) == IAerosolClassHydrophobic) then
625          ext = ext + mixing_ratio(jcol,jtype) &
626               &    * ao%mass_ext_sw_phobic(iband,ao%itype(jtype))
627        else if (ao%iclass(jtype) == IAerosolClassHydrophilic) then
628          ext = ext + mixing_ratio(jcol,jtype) &
629               &    * ao%mass_ext_sw_philic(iband,irh,ao%itype(jtype))
630        end if
631      end do
632
633      extinction(jcol) = ext
634    end do
635
636    if (lhook) call dr_hook('radiation_aerosol_optics:aerosol_sw_extinction',1,hook_handle)
637
638  end subroutine aerosol_sw_extinction
639
640end module radiation_aerosol_optics
Note: See TracBrowser for help on using the repository browser.