source: LMDZ6/trunk/libf/phylmd/ecrad/radiation_single_level.F90 @ 4013

Last change on this file since 4013 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: 14.2 KB
RevLine 
[3908]1! radiation_single_level.F90 - Derived type for single-level fields
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!   2019-01-14  R. Hogan  Added out_of_physical_bounds routine
17!   2019-01-18  R. Hogan  Added weighted albedo mapping
18
19module radiation_single_level
20
21  use parkind1, only : jprb
22
23  implicit none
24  public
25
26  !---------------------------------------------------------------------
27  ! Derived type to contain variables that don't vary with height;
28  ! mostly they are surface properties
29  type single_level_type
30    ! Note that skin temperature is only defined if
31    ! is_simple_surface==.true.
32    real(jprb), allocatable, dimension(:) :: &
33         &   cos_sza, &       ! (ncol) Cosine of solar zenith angle
34         &   skin_temperature ! (ncol) Skin temperature (K)
35
36    ! Shortwave albedo: if sw_albedo_direct is not allocated then
37    ! sw_albedo will be used for both direct and diffuse solar
38    ! radiation; otherwise, sw_albedo will be used for diffuse
39    ! radiation and sw_albedo_direct for direct radiation.
40    real(jprb), allocatable, dimension(:,:) :: &
41         &   sw_albedo, &        ! (ncol,nalbedobands)
42         &   sw_albedo_direct    ! (ncol,nalbedobands)
43   
44    ! If is_simple_surface==.true., we provide longwave emissivity
45    ! coarser spectral intervals, while if
46    ! is_simple_surface==.false. it will be generated by the surface
47    ! radiation scheme and may optionally be on the full spectral grid
48    ! of the atmospheric model.
49    real(jprb), allocatable, dimension(:,:) :: &
50         &   lw_emissivity       ! (ncol,nemissbands) If
51    ! is_simple_surface==.false. then we specify the emission instead
52    ! of the skin temperature, and again this may be on the emissivity
53    ! bands or the full spectral grid
54    real(jprb), allocatable, dimension(:,:) :: &
55         &   lw_emission         ! (ncol,nemissbands)
56
57    ! Incoming solar irradiance at the Earth is specified as the total
58    ! across the shortwave spectrum. Note that this needs to be
59    ! adjusted for Earth-Sun distance, solar cycle etc. To get
60    ! normalized fluxes out, simply set this to 1.0.
61    real(jprb) :: solar_irradiance = 1366.0_jprb ! W m-2
62
63    ! If config%use_spectral_solar_irradiance==true then this will be
64    ! scaled by spectral_solar_scaling
65    real(jprb), allocatable, dimension(:) :: &
66         &   spectral_solar_scaling ! (n_bands_sw)
67 
68    ! Seed for random number generator in McICA; it is expected that
69    ! the host model generates this from the model time, longitude
70    ! and latitude, in order that the model is deterministic
71    integer, allocatable, dimension(:) :: iseed ! (ncol)
72
73    ! Is the underlying surface a simple single "flat" tile? If so we
74    ! describe it with skin_temperature and lw_emissivity. Otherwise
75    ! we describe it with lw_emission and lw_emissivity coming out of
76    ! the surface radiation scheme.
77    logical :: is_simple_surface = .true.
78
79    ! If is_simple_surface==.false., do we use the full number of g
80    ! points for dimensioning sw_albedo, sw_albedo_direct and
81    ! lw_emission?
82    !logical :: use_full_spectrum_sw = .false.
83    !logical :: use_full_spectrum_lw = .true.
84
85  contains
86    procedure :: allocate   => allocate_single_level
87    procedure :: deallocate => deallocate_single_level
88    procedure :: init_seed_simple
89    procedure :: get_albedos
90    procedure :: out_of_physical_bounds
91
92  end type single_level_type
93
94contains
95 
96  !---------------------------------------------------------------------
97  subroutine allocate_single_level(this, ncol, nalbedobands, nemisbands, &
98       &                           use_sw_albedo_direct, is_simple_surface)
99
100    use yomhook, only : lhook, dr_hook
101
102    class(single_level_type), intent(inout) :: this
103    integer,                  intent(in)    :: ncol, nalbedobands, nemisbands
104    logical,        optional, intent(in)    :: use_sw_albedo_direct
105    logical,        optional, intent(in)    :: is_simple_surface
106
107    real(jprb) :: hook_handle
108
109    if (lhook) call dr_hook('radiation_single_level:allocate',0,hook_handle)
110
111    call this%deallocate()
112
113    if (present(is_simple_surface)) then
114      this%is_simple_surface = is_simple_surface
115    else
116      this%is_simple_surface = .true.
117    end if
118
119    allocate(this%cos_sza(ncol))
120
121    if (this%is_simple_surface) then
122      allocate(this%skin_temperature(ncol))
123    else
124      allocate(this%lw_emission(ncol, nemisbands))
125    end if
126    allocate(this%lw_emissivity(ncol, nemisbands))
127
128    allocate(this%sw_albedo(ncol, nalbedobands))
129
130    if (present(use_sw_albedo_direct)) then
131      if (use_sw_albedo_direct) then
132        allocate(this%sw_albedo_direct(ncol, nalbedobands))
133      end if
134    end if
135
136    allocate(this%iseed(ncol))
137
138    if (lhook) call dr_hook('radiation_single_level:allocate',1,hook_handle)
139
140  end subroutine allocate_single_level
141
142
143  !---------------------------------------------------------------------
144  subroutine deallocate_single_level(this)
145
146    use yomhook, only : lhook, dr_hook
147
148    class(single_level_type), intent(inout) :: this
149
150    real(jprb) :: hook_handle
151
152    if (lhook) call dr_hook('radiation_single_level:deallocate',0,hook_handle)
153
154    if (allocated(this%cos_sza)) then
155      deallocate(this%cos_sza)
156    end if
157    if (allocated(this%skin_temperature)) then
158      deallocate(this%skin_temperature)
159    end if
160    if (allocated(this%sw_albedo)) then
161      deallocate(this%sw_albedo)
162    end if
163    if (allocated(this%sw_albedo_direct)) then
164      deallocate(this%sw_albedo_direct)
165    end if
166    if (allocated(this%lw_emissivity)) then
167      deallocate(this%lw_emissivity)
168    end if
169    if (allocated(this%lw_emission)) then
170      deallocate(this%lw_emission)
171    end if
172    if (allocated(this%spectral_solar_scaling)) then
173      deallocate(this%spectral_solar_scaling)
174    end if
175    if (allocated(this%iseed)) then
176      deallocate(this%iseed)
177    end if
178
179    if (lhook) call dr_hook('radiation_single_level:deallocate',1,hook_handle)
180
181  end subroutine deallocate_single_level
182
183
184  !---------------------------------------------------------------------
185  ! Unimaginative initialization of random-number seeds
186  subroutine init_seed_simple(this, istartcol, iendcol)
187    class(single_level_type), intent(inout) :: this
188    integer, intent(in)                     :: istartcol, iendcol
189
190    integer :: jcol
191
192    if (.not. allocated(this%iseed)) then
193      allocate(this%iseed(istartcol:iendcol))
194    end if
195
196    do jcol = istartcol,iendcol
197      this%iseed(jcol) = jcol
198    end do
199  end subroutine init_seed_simple
200
201
202  !---------------------------------------------------------------------
203  ! Extract the shortwave and longwave surface albedos in each g-point
204  subroutine get_albedos(this, istartcol, iendcol, config, &
205       &                 sw_albedo_direct, sw_albedo_diffuse, lw_albedo)
206
207    use radiation_config, only : config_type
208    use radiation_io,     only : nulerr, radiation_abort
209    use yomhook,          only : lhook, dr_hook
210
211    class(single_level_type), intent(in) :: this
212    type(config_type),        intent(in) :: config
213    integer,                  intent(in) :: istartcol, iendcol
214
215    ! The longwave albedo of the surface in each longwave g-point;
216    ! note that these are weighted averages of the values from
217    ! individual tiles
218    real(jprb), intent(out), optional &
219         &  :: lw_albedo(config%n_g_lw, istartcol:iendcol)
220
221    ! Direct and diffuse shortwave surface albedo in each shortwave
222    ! g-point; note that these are weighted averages of the values
223    ! from individual tiles
224    real(jprb), intent(out), dimension(config%n_g_sw, istartcol:iendcol) &
225         &  :: sw_albedo_direct, sw_albedo_diffuse
226
227    ! Temporary storage of albedo in ecRad bands
228    real(jprb) :: sw_albedo_band(istartcol:iendcol, config%n_bands_sw)
229    real(jprb) :: lw_albedo_band (istartcol:iendcol, config%n_bands_lw)
230
231    ! Number of albedo bands
232    integer :: nalbedoband
233
234    ! Loop indices for ecRad bands and albedo bands
235    integer :: jband, jalbedoband
236
237    real(jprb) :: hook_handle
238
239    if (lhook) call dr_hook('radiation_single_level:get_albedos',0,hook_handle)
240
241    ! Albedos/emissivities are stored in single_level in their own
242    ! spectral intervals and with column as the first dimension
243    if (config%use_canopy_full_spectrum_sw) then
244      ! Albedos provided in each g point
245      sw_albedo_diffuse = transpose(this%sw_albedo(istartcol:iendcol,:))
246      if (allocated(this%sw_albedo_direct)) then
247        sw_albedo_direct = transpose(this%sw_albedo_direct(istartcol:iendcol,:))
248      end if
249    elseif (.not. config%do_nearest_spectral_sw_albedo) then
250      ! Albedos averaged accurately to ecRad spectral bands
251      nalbedoband = size(config%sw_albedo_weights,1)
252      sw_albedo_band = 0.0_jprb
253      do jband = 1,config%n_bands_sw
254        do jalbedoband = 1,nalbedoband
255          if (config%sw_albedo_weights(jalbedoband,jband) /= 0.0_jprb) then
256            sw_albedo_band(istartcol:iendcol,jband) &
257                 &  = sw_albedo_band(istartcol:iendcol,jband) &
258                 &  + config%sw_albedo_weights(jalbedoband,jband) &
259                 &    * this%sw_albedo(istartcol:iendcol, jalbedoband)
260          end if
261        end do
262      end do
263
264      sw_albedo_diffuse = transpose(sw_albedo_band(istartcol:iendcol, &
265           &                              config%i_band_from_reordered_g_sw))
266      if (allocated(this%sw_albedo_direct)) then
267        sw_albedo_band = 0.0_jprb
268        do jband = 1,config%n_bands_sw
269          do jalbedoband = 1,nalbedoband
270            if (config%sw_albedo_weights(jalbedoband,jband) /= 0.0_jprb) then
271              sw_albedo_band(istartcol:iendcol,jband) &
272                   &  = sw_albedo_band(istartcol:iendcol,jband) &
273                   &  + config%sw_albedo_weights(jalbedoband,jband) &
274                   &    * this%sw_albedo_direct(istartcol:iendcol, jalbedoband)
275            end if
276          end do
277        end do
278        sw_albedo_direct = transpose(sw_albedo_band(istartcol:iendcol, &
279             &                             config%i_band_from_reordered_g_sw))
280      else
281        sw_albedo_direct = sw_albedo_diffuse
282      end if
283    else
284      ! Albedos mapped less accurately to ecRad spectral bands
285      sw_albedo_diffuse = transpose(this%sw_albedo(istartcol:iendcol, &
286           &  config%i_albedo_from_band_sw(config%i_band_from_reordered_g_sw)))
287      if (allocated(this%sw_albedo_direct)) then
288        sw_albedo_direct = transpose(this%sw_albedo_direct(istartcol:iendcol, &
289             &  config%i_albedo_from_band_sw(config%i_band_from_reordered_g_sw)))
290      else
291        sw_albedo_direct = sw_albedo_diffuse
292      end if
293    end if
294
295    if (present(lw_albedo)) then
296      if (config%use_canopy_full_spectrum_lw) then
297        if (config%n_g_lw /= size(this%lw_emissivity,2)) then
298          write(nulerr,'(a)') '*** Error: single_level%lw_emissivity has the wrong number of spectral intervals'
299          call radiation_abort()   
300        end if
301        lw_albedo = 1.0_jprb - transpose(this%lw_emissivity(istartcol:iendcol,:))
302      else if (.not. config%do_nearest_spectral_lw_emiss) then
303        ! Albedos averaged accurately to ecRad spectral bands
304        nalbedoband = size(config%lw_emiss_weights,1)
305        lw_albedo_band = 0.0_jprb
306        do jband = 1,config%n_bands_lw
307          do jalbedoband = 1,nalbedoband
308            if (config%lw_emiss_weights(jalbedoband,jband) /= 0.0_jprb) then
309              lw_albedo_band(istartcol:iendcol,jband) &
310                   &  = lw_albedo_band(istartcol:iendcol,jband) &
311                   &  + config%lw_emiss_weights(jalbedoband,jband) &
312                   &    * (1.0_jprb-this%lw_emissivity(istartcol:iendcol, jalbedoband))
313            end if
314          end do
315        end do
316
317        lw_albedo = transpose(lw_albedo_band(istartcol:iendcol, &
318             &                config%i_band_from_reordered_g_lw))
319      else
320        lw_albedo = 1.0_jprb - transpose(this%lw_emissivity(istartcol:iendcol, &
321             &  config%i_emiss_from_band_lw(config%i_band_from_reordered_g_lw)))
322      end if
323    end if
324
325    if (lhook) call dr_hook('radiation_single_level:get_albedos',1,hook_handle)
326
327  end subroutine get_albedos
328
329
330  !---------------------------------------------------------------------
331  ! Return .true. if the contents of a single_level structure are out
332  ! of a physically sensible range, optionally only considering only
333  ! columns between istartcol and iendcol
334  function out_of_physical_bounds(this, istartcol, iendcol, do_fix) result(is_bad)
335
336    use yomhook,          only : lhook, dr_hook
337    use radiation_config, only : out_of_bounds_1d, out_of_bounds_2d
338
339    class(single_level_type), intent(inout) :: this
340    integer,         optional,intent(in) :: istartcol, iendcol
341    logical,         optional,intent(in) :: do_fix
342    logical                              :: is_bad
343
344    logical    :: do_fix_local
345
346    real(jprb) :: hook_handle
347
348    if (lhook) call dr_hook('radiation_single_level:out_of_physical_bounds',0,hook_handle)
349
350    if (present(do_fix)) then
351      do_fix_local = do_fix
352    else
353      do_fix_local = .false.
354    end if
355
356    is_bad =    out_of_bounds_1d(this%cos_sza, "cos_sza", -1.0_jprb, 1.0_jprb, &
357         &                       do_fix_local, i1=istartcol, i2=iendcol) &
358         & .or. out_of_bounds_1d(this%skin_temperature, "skin_temperature", 173.0_jprb, 373.0_jprb, &
359         &                       do_fix_local, i1=istartcol, i2=iendcol) &
360         & .or. out_of_bounds_2d(this%sw_albedo, "sw_albedo", 0.0_jprb, 1.0_jprb, &
361         &                       do_fix_local, i1=istartcol, i2=iendcol) &
362         & .or. out_of_bounds_2d(this%sw_albedo_direct, "sw_albedo", 0.0_jprb, 1.0_jprb, &
363         &                       do_fix_local, i1=istartcol, i2=iendcol) &
364         & .or. out_of_bounds_2d(this%lw_emissivity, "lw_emissivity", 0.0_jprb, 1.0_jprb, &
365         &                       do_fix_local, i1=istartcol, i2=iendcol)
366
367    if (lhook) call dr_hook('radiation_single_level:out_of_physical_bounds',1,hook_handle)
368
369  end function out_of_physical_bounds
370
371end module radiation_single_level
Note: See TracBrowser for help on using the repository browser.