source: LMDZ6/trunk/libf/phylmd/ecrad/radiation_monochromatic.F90 @ 4190

Last change on this file since 4190 was 3908, checked in by idelkadi, 3 years ago

Online implementation of the radiative transfer code ECRAD in the LMDZ model.

  • Inclusion of the ecrad directory containing the sources of the ECRAD code
    • interface routine : radiation_scheme.F90
  • Adaptation of compilation scripts :
    • compilation under CPP key CPP_ECRAD
    • compilation with option "-rad ecard" or "-ecard true"
    • The "-rad old/rtm/ecran" build option will need to replace the "-rrtm true" and "-ecrad true" options in the future.
  • Runing LMDZ simulations with ecrad, you need :
    • logical key iflag_rrtm = 2 in physiq.def
    • namelist_ecrad (DefLists?)
    • the directory "data" containing the configuration files is temporarily placed in ../libfphylmd/ecrad/
  • Compilation and execution are tested in the 1D case. The repository under svn would allow to continue the implementation work: tests, verification of the results, ...
File size: 15.2 KB
Line 
1! radiation_interface.F90 - Monochromatic gas/cloud optics for testing
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-04-11  R. Hogan  Receive "surface" dummy argument
17!   2017-09-13  R. Hogan  Revert
18!   2018-08-29  R. Hogan  Particulate single-scattering albedo / asymmetry from namelist
19
20module radiation_monochromatic
21
22  implicit none
23
24  public  :: setup_gas_optics, gas_optics, set_gas_units, &
25       &     setup_cloud_optics, cloud_optics,            &
26       &     setup_aerosol_optics, add_aerosol_optics
27
28contains
29
30  ! Provides elemental function "delta_eddington"
31#include "radiation_delta_eddington.h"
32
33  !---------------------------------------------------------------------
34  ! Setup the arrays in the config object corresponding to the
35  ! monochromatic gas optics model.  The directory argument is not
36  ! used, since no look-up tables need to be loaded.
37  subroutine setup_gas_optics(config, directory)
38
39    use radiation_config, only : config_type
40   
41    type(config_type), intent(inout) :: config
42    character(len=*),  intent(in)    :: directory
43
44    ! In the monochromatic model we have simply one band and g-point
45    ! in both the longwave and shortwave parts of the spectrum
46    config%n_g_sw     = 1
47    config%n_g_lw     = 1
48    config%n_bands_sw = 1
49    config%n_bands_lw = 1
50
51    ! Allocate arrays
52    allocate(config%i_band_from_g_sw          (config%n_g_sw))
53    allocate(config%i_band_from_g_lw          (config%n_g_lw))
54    allocate(config%i_band_from_reordered_g_sw(config%n_g_sw))
55    allocate(config%i_band_from_reordered_g_lw(config%n_g_lw))
56    allocate(config%i_g_from_reordered_g_sw(config%n_g_sw))
57    allocate(config%i_g_from_reordered_g_lw(config%n_g_lw))
58
59    ! Indices are trivial...
60    config%i_band_from_g_sw           = 1
61    config%i_band_from_g_lw           = 1
62    config%i_g_from_reordered_g_sw    = 1
63    config%i_g_from_reordered_g_lw    = 1
64    config%i_band_from_reordered_g_sw = 1
65    config%i_band_from_reordered_g_lw = 1
66
67  end subroutine setup_gas_optics
68
69
70  !---------------------------------------------------------------------
71  ! Dummy routine for scaling gas mixing ratios
72  subroutine set_gas_units(gas)
73
74    use radiation_gas,           only : gas_type
75    type(gas_type),    intent(inout) :: gas
76
77  end subroutine set_gas_units
78
79
80  !---------------------------------------------------------------------
81  ! Dummy setup routine for cloud optics: in fact, no setup is
82  ! required for monochromatic case
83  subroutine setup_cloud_optics(config)
84
85    use radiation_config, only : config_type
86    type(config_type), intent(inout) :: config
87
88  end subroutine setup_cloud_optics
89
90
91  !---------------------------------------------------------------------
92  ! Dummy subroutine since no aerosols are represented in
93  ! monochromatic case
94  subroutine setup_aerosol_optics(config)
95
96    use radiation_config,              only : config_type
97    type(config_type), intent(inout) :: config
98
99  end subroutine setup_aerosol_optics
100
101
102  !---------------------------------------------------------------------
103  ! Compute gas optical depths, shortwave scattering, Planck function
104  ! and incoming shortwave radiation at top-of-atmosphere
105  subroutine gas_optics(ncol,nlev,istartcol,iendcol, &
106       config, single_level, thermodynamics, gas, lw_albedo, &
107       od_lw, od_sw, ssa_sw, planck_hl, lw_emission, &
108       incoming_sw)
109
110    use parkind1,                 only : jprb
111    use radiation_config,         only : config_type
112    use radiation_thermodynamics, only : thermodynamics_type
113    use radiation_single_level,   only : single_level_type
114    use radiation_gas,            only : gas_type
115    use radiation_constants,      only : Pi, StefanBoltzmann
116
117    ! Inputs
118    integer, intent(in) :: ncol               ! number of columns
119    integer, intent(in) :: nlev               ! number of model levels
120    integer, intent(in) :: istartcol, iendcol ! range of columns to process
121    type(config_type),        intent(in) :: config
122    type(single_level_type),  intent(in) :: single_level
123    type(thermodynamics_type),intent(in) :: thermodynamics
124    type(gas_type),           intent(in) :: gas
125
126    ! Longwave albedo of the surface
127    real(jprb), dimension(config%n_g_lw,istartcol:iendcol), &
128         &  intent(in) :: lw_albedo
129
130    ! Outputs
131
132    ! Gaseous layer optical depth in longwave and shortwave, and
133    ! shortwave single scattering albedo (i.e. fraction of extinction
134    ! due to Rayleigh scattering) at each g-point
135    real(jprb), dimension(config%n_g_lw,nlev,istartcol:iendcol), intent(out) :: &
136         &   od_lw
137    real(jprb), dimension(config%n_g_sw,nlev,istartcol:iendcol), intent(out) :: &
138         &   od_sw, ssa_sw
139
140    ! The Planck function (emitted flux from a black body) at half
141    ! levels and at the surface at each longwave g-point
142    real(jprb), dimension(config%n_g_lw,nlev+1,istartcol:iendcol), intent(out) :: &
143         &   planck_hl
144    real(jprb), dimension(config%n_g_lw,istartcol:iendcol), intent(out) :: &
145         &   lw_emission
146
147    ! The incoming shortwave flux into a plane perpendicular to the
148    ! incoming radiation at top-of-atmosphere in each of the shortwave
149    ! g-points
150    real(jprb), dimension(config%n_g_sw,istartcol:iendcol), intent(out) :: &
151         &   incoming_sw
152   
153    ! Ratio of the optical depth of the entire atmospheric column that
154    ! is in the current layer
155    real(jprb), dimension(istartcol:iendcol) :: extinction_fraction
156
157    ! In the monochromatic model, the absorption by the atmosphere is
158    ! assumed proportional to the mass in each layer, so is defined in
159    ! terms of a total zenith optical depth and then distributed with
160    ! height according to the pressure.
161    !real(jprb), parameter :: total_od_sw = 0.10536_jprb ! Transmittance 0.9
162    !real(jprb), parameter :: total_od_lw = 1.6094_jprb  ! Transmittance 0.2
163
164    integer :: jlev
165
166    do jlev = 1,nlev
167      ! The fraction of the total optical depth in the current layer
168      ! is proportional to the fraction of the mass of the atmosphere
169      ! in the current layer, computed from pressure assuming
170      ! hydrostatic balance
171      extinction_fraction = &
172           &   (thermodynamics%pressure_hl(istartcol:iendcol,jlev+1) &
173           &   -thermodynamics%pressure_hl(istartcol:iendcol,jlev)) &
174           &   /thermodynamics%pressure_hl(istartcol:iendcol,nlev)
175     
176      ! Assign longwave and shortwave optical depths
177      od_lw(1,jlev,:) = config%mono_lw_total_od*extinction_fraction
178      od_sw(1,jlev,:) = config%mono_sw_total_od*extinction_fraction
179    end do
180
181    ! Shortwave single-scattering albedo is almost entirely Rayleigh
182    ! scattering
183    ssa_sw = 0.999999_jprb
184
185    ! Entire shortwave spectrum represented in one band
186    incoming_sw(1,:) = single_level%solar_irradiance
187
188    if (single_level%is_simple_surface) then
189      if (config%mono_lw_wavelength <= 0.0_jprb) then
190        ! Entire longwave spectrum represented in one band
191        lw_emission(1,istartcol:iendcol) &
192             &  = StefanBoltzmann * single_level%skin_temperature(istartcol:iendcol)**4 &
193             &  * single_level%lw_emissivity(istartcol:iendcol,1)
194        do jlev = 1,nlev+1
195          planck_hl(1,jlev,istartcol:iendcol) = StefanBoltzmann * thermodynamics%temperature_hl(istartcol:iendcol,jlev)**4
196        end do
197      else
198        ! Single wavelength: multiply by pi to convert W sr-1 m-3 to W m-3
199        lw_emission(1,istartcol:iendcol) = Pi*planck_function(config%mono_lw_wavelength, &
200             &             single_level%skin_temperature(istartcol:iendcol)) &
201             &  * single_level%lw_emissivity(istartcol:iendcol,1)
202        do jlev = 1,nlev+1
203          planck_hl(1,jlev,istartcol:iendcol) = Pi*planck_function(config%mono_lw_wavelength, &
204               &             thermodynamics%temperature_hl(istartcol:iendcol,jlev))
205        end do
206      end if
207    else
208      lw_emission = transpose(single_level%lw_emission)
209    end if
210
211  end subroutine gas_optics
212
213
214  !---------------------------------------------------------------------
215  ! Compute cloud optical depth, single-scattering albedo and
216  ! g factor in the longwave and shortwave
217  subroutine cloud_optics(nlev,istartcol,iendcol, &
218       &   config, thermodynamics, cloud, &
219       &   od_lw_cloud, ssa_lw_cloud, g_lw_cloud, &
220       &   od_sw_cloud, ssa_sw_cloud, g_sw_cloud)
221
222    use parkind1,                 only : jprb
223    use radiation_config,         only : config_type
224    use radiation_thermodynamics, only : thermodynamics_type
225    use radiation_cloud,          only : cloud_type
226    use radiation_constants,      only : AccelDueToGravity, &
227         &   DensityLiquidWater, DensitySolidIce
228
229    ! Inputs
230    integer, intent(in) :: nlev               ! number of model levels
231    integer, intent(in) :: istartcol, iendcol ! range of columns to process
232    type(config_type), intent(in) :: config
233    type(thermodynamics_type),intent(in) :: thermodynamics
234    type(cloud_type),   intent(in) :: cloud
235
236    ! Outputs
237
238    ! Layer optical depth, single scattering albedo and g factor of
239    ! clouds in each longwave band, where the latter two
240    ! variables are only defined if cloud longwave scattering is
241    ! enabled (otherwise both are treated as zero).
242    real(jprb), dimension(config%n_bands_lw,nlev,istartcol:iendcol), intent(out) :: &
243         &   od_lw_cloud
244    real(jprb), dimension(config%n_bands_lw_if_scattering,nlev,istartcol:iendcol), &
245         &   intent(out) :: ssa_lw_cloud, g_lw_cloud
246
247    ! Layer optical depth, single scattering albedo and g factor of
248    ! clouds in each shortwave band
249    real(jprb), dimension(config%n_g_sw,nlev,istartcol:iendcol), intent(out) :: &
250         &   od_sw_cloud, ssa_sw_cloud, g_sw_cloud
251
252    ! In-cloud liquid and ice water path in a layer, in kg m-2
253    real(jprb), dimension(nlev,istartcol:iendcol) :: lwp_kg_m2, iwp_kg_m2
254
255    integer  :: jlev, jcol
256
257    ! Factor to convert from gridbox-mean mass mixing ratio to
258    ! in-cloud water path
259    real(jprb) :: factor
260
261    ! Convert cloud mixing ratio into liquid and ice water path
262    ! in each layer
263    do jlev = 1, nlev
264      do jcol = istartcol, iendcol
265        ! Factor to convert from gridbox-mean mass mixing ratio to
266        ! in-cloud water path involves the pressure difference in
267        ! Pa, acceleration due to gravity and cloud fraction
268        ! adjusted to avoid division by zero.
269        factor = ( thermodynamics%pressure_hl(jcol,jlev+1)    &
270             &    -thermodynamics%pressure_hl(jcol,jlev  )  ) &
271             &   / (AccelDueToGravity &
272             &   * max(epsilon(1.0_jprb), cloud%fraction(jcol,jlev)))
273        lwp_kg_m2(jlev,jcol) = factor * cloud%q_liq(jcol,jlev)
274        iwp_kg_m2(jlev,jcol) = factor * cloud%q_ice(jcol,jlev)
275      end do
276    end do
277
278    ! Geometric optics approximation: particles treated as much larger
279    ! than the wavelength in both shortwave and longwave
280    od_sw_cloud(1,:,:) &
281         &   = (3.0_jprb/(2.0_jprb*DensityLiquidWater)) &
282         &   * lwp_kg_m2 / transpose(cloud%re_liq(istartcol:iendcol,:)) &
283         &   + (3.0_jprb / (2.0_jprb * DensitySolidIce)) &
284         &   * iwp_kg_m2 / transpose(cloud%re_ice(istartcol:iendcol,:))
285    od_lw_cloud(1,:,:) = lwp_kg_m2 * 137.22_jprb &
286         &   + (3.0_jprb / (2.0_jprb * DensitySolidIce)) &
287         &   * iwp_kg_m2 / transpose(cloud%re_ice(istartcol:iendcol,:))
288
289    if (config%iverbose >= 4) then
290      do jcol = istartcol,iendcol
291        write(*,'(a,i0,a,f7.3,a,f7.3)') 'Profile ', jcol, ': shortwave optical depth = ', &
292             &  sum(od_sw_cloud(1,:,jcol)*cloud%fraction(jcol,:)), &
293             &  ', longwave optical depth = ', &
294             &  sum(od_lw_cloud(1,:,jcol)*cloud%fraction(jcol,:))
295        !    print *, 'LWP = ', sum(lwp_kg_m2(:,istartcol)*cloud%fraction(istartcol,:))
296      end do
297    end if
298
299    ssa_sw_cloud = config%mono_sw_single_scattering_albedo
300    g_sw_cloud   = config%mono_sw_asymmetry_factor
301
302    ! In-place delta-Eddington scaling
303    call delta_eddington(od_sw_cloud, ssa_sw_cloud, g_sw_cloud)
304
305    if (config%do_lw_cloud_scattering) then
306      ssa_lw_cloud = config%mono_lw_single_scattering_albedo
307      g_lw_cloud   = config%mono_lw_asymmetry_factor
308      ! In-place delta-Eddington scaling
309      call delta_eddington(od_lw_cloud, ssa_lw_cloud, g_lw_cloud)
310    end if
311
312  end subroutine cloud_optics
313
314
315  !---------------------------------------------------------------------
316  ! Dummy subroutine since no aerosols are represented in
317  ! monochromatic case
318  subroutine add_aerosol_optics(nlev,istartcol,iendcol, &
319       &  config, thermodynamics, gas, aerosol, &
320       &  od_lw, ssa_lw, g_lw, od_sw, ssa_sw, g_sw)
321
322    use parkind1,                      only : jprb
323
324    use radiation_config,              only : config_type
325    use radiation_thermodynamics,      only : thermodynamics_type
326    use radiation_gas,                 only : gas_type
327    use radiation_aerosol,             only : aerosol_type
328
329    integer, intent(in) :: nlev               ! number of model levels
330    integer, intent(in) :: istartcol, iendcol ! range of columns to process
331    type(config_type), intent(in), target :: config
332    type(thermodynamics_type),intent(in)  :: thermodynamics
333    type(gas_type),           intent(in)  :: gas
334    type(aerosol_type),       intent(in)  :: aerosol
335    ! Optical depth, single scattering albedo and asymmetry factor of
336    ! the atmosphere (gases on input, gases and aerosols on output)
337    ! for each g point. Note that longwave ssa and asymmetry and
338    ! shortwave asymmetry are all zero for gases, so are not yet
339    ! defined on input and are therefore intent(out).
340    real(jprb), dimension(config%n_g_lw,nlev,istartcol:iendcol), intent(inout) :: od_lw
341    real(jprb), dimension(config%n_g_lw_if_scattering,nlev,istartcol:iendcol), &
342         &  intent(out) :: ssa_lw, g_lw
343    real(jprb), dimension(config%n_g_sw,nlev,istartcol:iendcol), intent(inout) &
344         &  :: od_sw, ssa_sw
345    real(jprb), dimension(config%n_g_sw,nlev,istartcol:iendcol), intent(out) :: g_sw
346
347    g_sw = 0.0_jprb
348
349    if (config%do_lw_aerosol_scattering) then
350      ssa_lw = 0.0_jprb
351      g_lw   = 0.0_jprb
352    end if
353
354  end subroutine add_aerosol_optics
355
356  !---------------------------------------------------------------------
357  ! Planck function in terms of wavelength
358  elemental function planck_function(wavelength, temperature)
359
360    use parkind1,            only : jprb
361
362    use radiation_constants, only : BoltzmannConstant, PlanckConstant, &
363         &                          SpeedOfLight
364
365    real(jprb), intent(in) :: wavelength  ! metres
366    real(jprb), intent(in) :: temperature ! Kelvin
367
368    ! Output in W sr-1 m-3
369    real(jprb)             :: planck_function
370
371    if (temperature > 0.0_jprb) then
372      planck_function = 2.0_jprb * PlanckConstant * SpeedOfLight**2 &
373           &   / (wavelength**5 &
374           &   * (exp(PlanckConstant*SpeedOfLight &
375           &         /(wavelength*BoltzmannConstant*temperature)) - 1.0_jprb))
376    else
377      planck_function = 0.0_jprb
378    end if
379
380  end function planck_function
381
382end module radiation_monochromatic
Note: See TracBrowser for help on using the repository browser.