source: LMDZ6/trunk/libf/phylmd/ecrad.v1.5.1/radiation_single_level.F90 @ 5407

Last change on this file since 5407 was 4489, checked in by lguez, 21 months ago

Merge LMDZ_ECRad branch back into trunk!

File size: 15.3 KB
Line 
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  ! Allocate the arrays of a single-level type
98  subroutine allocate_single_level(this, ncol, nalbedobands, nemisbands, &
99       &                           use_sw_albedo_direct, is_simple_surface)
100
101    use yomhook, only : lhook, dr_hook
102
103    class(single_level_type), intent(inout) :: this
104    integer,                  intent(in)    :: ncol, nalbedobands, nemisbands
105    logical,        optional, intent(in)    :: use_sw_albedo_direct
106    logical,        optional, intent(in)    :: is_simple_surface
107
108    real(jprb) :: hook_handle
109
110    if (lhook) call dr_hook('radiation_single_level:allocate',0,hook_handle)
111
112    call this%deallocate()
113
114    if (present(is_simple_surface)) then
115      this%is_simple_surface = is_simple_surface
116    else
117      this%is_simple_surface = .true.
118    end if
119
120    allocate(this%cos_sza(ncol))
121
122    if (this%is_simple_surface) then
123      allocate(this%skin_temperature(ncol))
124    else
125      allocate(this%lw_emission(ncol, nemisbands))
126    end if
127    allocate(this%lw_emissivity(ncol, nemisbands))
128
129    allocate(this%sw_albedo(ncol, nalbedobands))
130
131    if (present(use_sw_albedo_direct)) then
132      if (use_sw_albedo_direct) then
133        allocate(this%sw_albedo_direct(ncol, nalbedobands))
134      end if
135    end if
136
137    allocate(this%iseed(ncol))
138
139    if (lhook) call dr_hook('radiation_single_level:allocate',1,hook_handle)
140
141  end subroutine allocate_single_level
142
143
144  !---------------------------------------------------------------------
145  ! Deallocate the arrays of a single-level type
146  subroutine deallocate_single_level(this)
147
148    use yomhook, only : lhook, dr_hook
149
150    class(single_level_type), intent(inout) :: this
151
152    real(jprb) :: hook_handle
153
154    if (lhook) call dr_hook('radiation_single_level:deallocate',0,hook_handle)
155
156    if (allocated(this%cos_sza)) then
157      deallocate(this%cos_sza)
158    end if
159    if (allocated(this%skin_temperature)) then
160      deallocate(this%skin_temperature)
161    end if
162    if (allocated(this%sw_albedo)) then
163      deallocate(this%sw_albedo)
164    end if
165    if (allocated(this%sw_albedo_direct)) then
166      deallocate(this%sw_albedo_direct)
167    end if
168    if (allocated(this%lw_emissivity)) then
169      deallocate(this%lw_emissivity)
170    end if
171    if (allocated(this%lw_emission)) then
172      deallocate(this%lw_emission)
173    end if
174    if (allocated(this%spectral_solar_scaling)) then
175      deallocate(this%spectral_solar_scaling)
176    end if
177    if (allocated(this%iseed)) then
178      deallocate(this%iseed)
179    end if
180
181    if (lhook) call dr_hook('radiation_single_level:deallocate',1,hook_handle)
182
183  end subroutine deallocate_single_level
184
185
186  !---------------------------------------------------------------------
187  ! Unimaginative initialization of random-number seeds
188  subroutine init_seed_simple(this, istartcol, iendcol)
189    class(single_level_type), intent(inout) :: this
190    integer, intent(in)                     :: istartcol, iendcol
191
192    integer :: jcol
193
194    if (.not. allocated(this%iseed)) then
195      allocate(this%iseed(istartcol:iendcol))
196    end if
197
198    do jcol = istartcol,iendcol
199      this%iseed(jcol) = jcol
200    end do
201  end subroutine init_seed_simple
202
203
204  !---------------------------------------------------------------------
205  ! Extract the shortwave and longwave surface albedos in each g-point
206  subroutine get_albedos(this, istartcol, iendcol, config, &
207       &                 sw_albedo_direct, sw_albedo_diffuse, lw_albedo)
208
209    use radiation_config, only : config_type
210    use radiation_io,     only : nulerr, radiation_abort
211    use yomhook,          only : lhook, dr_hook
212
213    class(single_level_type), intent(in) :: this
214    type(config_type),        intent(in) :: config
215    integer,                  intent(in) :: istartcol, iendcol
216
217    ! The longwave albedo of the surface in each longwave g-point;
218    ! note that these are weighted averages of the values from
219    ! individual tiles
220    real(jprb), intent(out), optional &
221         &  :: lw_albedo(config%n_g_lw, istartcol:iendcol)
222
223    ! Direct and diffuse shortwave surface albedo in each shortwave
224    ! g-point; note that these are weighted averages of the values
225    ! from individual tiles
226    real(jprb), intent(out), dimension(config%n_g_sw, istartcol:iendcol) &
227         &  :: sw_albedo_direct, sw_albedo_diffuse
228
229    ! Temporary storage of albedo in ecRad bands
230    real(jprb) :: sw_albedo_band(istartcol:iendcol, config%n_bands_sw)
231    real(jprb) :: lw_albedo_band(istartcol:iendcol, config%n_bands_lw)
232
233    ! Number of albedo bands
234    integer :: nalbedoband
235
236    ! Loop indices for ecRad bands and albedo bands
237    integer :: jband, jalbedoband, jcol
238
239    real(jprb) :: hook_handle
240
241    if (lhook) call dr_hook('radiation_single_level:get_albedos',0,hook_handle)
242
243    if (config%do_sw) then
244      ! Albedos/emissivities are stored in single_level in their own
245      ! spectral intervals and with column as the first dimension
246      if (config%use_canopy_full_spectrum_sw) then
247        ! Albedos provided in each g point
248        if (size(this%sw_albedo,2) /= config%n_g_sw) then
249          write(nulerr,'(a,i0,a)') '*** Error: single_level%sw_albedo does not have the expected ', &
250               &  config%n_g_sw, ' spectral intervals'
251          call radiation_abort()
252        end if
253        sw_albedo_diffuse = transpose(this%sw_albedo(istartcol:iendcol,:))
254        if (allocated(this%sw_albedo_direct)) then
255          sw_albedo_direct = transpose(this%sw_albedo_direct(istartcol:iendcol,:))
256        end if
257      else if (.not. config%do_nearest_spectral_sw_albedo) then
258        ! Albedos averaged accurately to ecRad spectral bands
259        nalbedoband = size(config%sw_albedo_weights,1)
260        if (size(this%sw_albedo,2) /= nalbedoband) then
261          write(nulerr,'(a,i0,a)') '*** Error: single_level%sw_albedo does not have the expected ', &
262               &  nalbedoband, ' bands'
263          call radiation_abort()
264        end if
265
266        sw_albedo_band = 0.0_jprb
267        do jband = 1,config%n_bands_sw
268          do jalbedoband = 1,nalbedoband
269            if (config%sw_albedo_weights(jalbedoband,jband) /= 0.0_jprb) then
270              do jcol = istartcol,iendcol
271                sw_albedo_band(jcol,jband) &
272                    &  = sw_albedo_band(jcol,jband) &
273                    &  + config%sw_albedo_weights(jalbedoband,jband) &
274                    &    * this%sw_albedo(jcol, jalbedoband)
275              end do
276            end if
277          end do
278        end do
279
280        sw_albedo_diffuse = transpose(sw_albedo_band(istartcol:iendcol, &
281             &                              config%i_band_from_reordered_g_sw))
282        if (allocated(this%sw_albedo_direct)) then
283          sw_albedo_band = 0.0_jprb
284          do jband = 1,config%n_bands_sw
285            do jalbedoband = 1,nalbedoband
286              if (config%sw_albedo_weights(jalbedoband,jband) /= 0.0_jprb) then
287                sw_albedo_band(istartcol:iendcol,jband) &
288                     &  = sw_albedo_band(istartcol:iendcol,jband) &
289                     &  + config%sw_albedo_weights(jalbedoband,jband) &
290                     &    * this%sw_albedo_direct(istartcol:iendcol, jalbedoband)
291              end if
292            end do
293          end do
294          sw_albedo_direct = transpose(sw_albedo_band(istartcol:iendcol, &
295               &                             config%i_band_from_reordered_g_sw))
296        else
297          sw_albedo_direct = sw_albedo_diffuse
298        end if
299      else
300        ! Albedos mapped less accurately to ecRad spectral bands
301        sw_albedo_diffuse = transpose(this%sw_albedo(istartcol:iendcol, &
302             &  config%i_albedo_from_band_sw(config%i_band_from_reordered_g_sw)))
303        if (allocated(this%sw_albedo_direct)) then
304          sw_albedo_direct = transpose(this%sw_albedo_direct(istartcol:iendcol, &
305               &  config%i_albedo_from_band_sw(config%i_band_from_reordered_g_sw)))
306        else
307          sw_albedo_direct = sw_albedo_diffuse
308        end if
309      end if
310    end if
311
312    if (config%do_lw .and. present(lw_albedo)) then
313      if (config%use_canopy_full_spectrum_lw) then
314        if (config%n_g_lw /= size(this%lw_emissivity,2)) then
315          write(nulerr,'(a,i0,a)') '*** Error: single_level%lw_emissivity does not have the expected ', &
316               &  config%n_g_lw, ' spectral intervals'
317          call radiation_abort()
318        end if
319        lw_albedo = 1.0_jprb - transpose(this%lw_emissivity(istartcol:iendcol,:))
320      else if (.not. config%do_nearest_spectral_lw_emiss) then
321        ! Albedos averaged accurately to ecRad spectral bands
322        nalbedoband = size(config%lw_emiss_weights,1)
323        if (nalbedoband /= size(this%lw_emissivity,2)) then
324          write(nulerr,'(a,i0,a)') '*** Error: single_level%lw_emissivity does not have the expected ', &
325               &  nalbedoband, ' bands'
326          call radiation_abort()
327        end if
328        lw_albedo_band = 0.0_jprb
329        do jband = 1,config%n_bands_lw
330          do jalbedoband = 1,nalbedoband
331            if (config%lw_emiss_weights(jalbedoband,jband) /= 0.0_jprb) then
332              do jcol = istartcol,iendcol
333                lw_albedo_band(jcol,jband) &
334                    &  = lw_albedo_band(jcol,jband) &
335                    &  + config%lw_emiss_weights(jalbedoband,jband) &
336                    &    * (1.0_jprb-this%lw_emissivity(jcol, jalbedoband))
337              end do
338            end if
339          end do
340        end do
341
342        lw_albedo = transpose(lw_albedo_band(istartcol:iendcol, &
343             &                config%i_band_from_reordered_g_lw))
344      else
345        lw_albedo = 1.0_jprb - transpose(this%lw_emissivity(istartcol:iendcol, &
346             &  config%i_emiss_from_band_lw(config%i_band_from_reordered_g_lw)))
347      end if
348    end if
349
350    if (lhook) call dr_hook('radiation_single_level:get_albedos',1,hook_handle)
351
352  end subroutine get_albedos
353
354
355  !---------------------------------------------------------------------
356  ! Return .true. if the contents of a single_level structure are out
357  ! of a physically sensible range, optionally only considering only
358  ! columns between istartcol and iendcol
359  function out_of_physical_bounds(this, istartcol, iendcol, do_fix) result(is_bad)
360
361    use yomhook,          only : lhook, dr_hook
362    use radiation_check,  only : out_of_bounds_1d, out_of_bounds_2d
363
364    class(single_level_type), intent(inout) :: this
365    integer,         optional,intent(in) :: istartcol, iendcol
366    logical,         optional,intent(in) :: do_fix
367    logical                              :: is_bad
368
369    logical    :: do_fix_local
370
371    real(jprb) :: hook_handle
372
373    if (lhook) call dr_hook('radiation_single_level:out_of_physical_bounds',0,hook_handle)
374
375    if (present(do_fix)) then
376      do_fix_local = do_fix
377    else
378      do_fix_local = .false.
379    end if
380
381    is_bad =    out_of_bounds_1d(this%cos_sza, "cos_sza", -1.0_jprb, 1.0_jprb, &
382         &                       do_fix_local, i1=istartcol, i2=iendcol) &
383         & .or. out_of_bounds_1d(this%skin_temperature, "skin_temperature", 173.0_jprb, 373.0_jprb, &
384         &                       do_fix_local, i1=istartcol, i2=iendcol) &
385         & .or. out_of_bounds_2d(this%sw_albedo, "sw_albedo", 0.0_jprb, 1.0_jprb, &
386         &                       do_fix_local, i1=istartcol, i2=iendcol) &
387         & .or. out_of_bounds_2d(this%sw_albedo_direct, "sw_albedo", 0.0_jprb, 1.0_jprb, &
388         &                       do_fix_local, i1=istartcol, i2=iendcol) &
389         & .or. out_of_bounds_2d(this%lw_emissivity, "lw_emissivity", 0.0_jprb, 1.0_jprb, &
390         &                       do_fix_local, i1=istartcol, i2=iendcol)
391
392    if (lhook) call dr_hook('radiation_single_level:out_of_physical_bounds',1,hook_handle)
393
394  end function out_of_physical_bounds
395
396end module radiation_single_level
Note: See TracBrowser for help on using the repository browser.