source: LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_ecckd_interface.F90 @ 4787

Last change on this file since 4787 was 4773, checked in by idelkadi, 12 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: 13.1 KB
Line 
1! radiation_ecckd_interface.F90 - Interface to ecCKD gas optics model
2!
3! (C) Copyright 2020- 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! License: see the COPYING file for details
15!
16
17module radiation_ecckd_interface
18
19  implicit none
20
21  public  :: setup_gas_optics, set_gas_units, gas_optics !, planck_function
22
23contains
24
25  !---------------------------------------------------------------------
26  ! Setup the ecCKD generalized gas optics model
27  subroutine setup_gas_optics(config)
28
29    use parkind1, only : jprb
30    use radiation_config
31    use yomhook,  only : lhook, dr_hook, jphook
32
33    type(config_type), intent(inout), target :: config
34
35    integer :: jj
36   
37    real(jphook) :: hook_handle
38
39    if (lhook) call dr_hook('radiation_ecckd_interface:setup_gas_optics',0,hook_handle)
40
41    if (config%do_sw) then
42
43      ! Read shortwave ecCKD gas optics NetCDF file
44      call config%gas_optics_sw%read(trim(config%gas_optics_sw_file_name), &
45           &                         config%iverbosesetup)
46
47      ! Copy over relevant properties
48      config%n_g_sw     = config%gas_optics_sw%ng
49
50      if (config%do_cloud_aerosol_per_sw_g_point) then
51        ! Bands and g points are the same
52        config%n_bands_sw = config%n_g_sw
53      else
54        ! Bands are groups of g points and span a continuous region of
55        ! wavenumber space
56        config%n_bands_sw = config%gas_optics_sw%spectral_def%nband
57      end if
58
59!      if (allocated(config%i_band_from_g_sw)) deallocate(config%i_band_from_g_sw)
60!      allocate(config%i_band_from_g_sw(config%n_g_sw))
61!      if (allocated(config%i_band_from_reordered_g_sw)) deallocate(config%i_band_from_reordered_g_sw)
62!      allocate(config%i_band_from_reordered_g_sw(config%n_g_sw))
63!      if (allocated(config%i_g_from_reordered_g_sw)) deallocate(config%i_g_from_reordered_g_sw)
64!      allocate(config%i_g_from_reordered_g_sw(config%n_g_sw))
65      if (.not.allocated(config%i_band_from_g_sw)) &
66              allocate(config%i_band_from_g_sw(config%n_g_sw))
67      if (.not.allocated(config%i_band_from_reordered_g_sw)) &
68              allocate(config%i_band_from_reordered_g_sw(config%n_g_sw))
69      if (.not.allocated(config%i_g_from_reordered_g_sw)) &
70              allocate(config%i_g_from_reordered_g_sw(config%n_g_sw))
71       
72      if (config%do_cloud_aerosol_per_sw_g_point) then
73        config%i_band_from_g_sw           = [ (jj, jj = 1,config%n_g_sw) ]
74        config%i_band_from_reordered_g_sw = [ (jj, jj = 1,config%n_g_sw) ]
75      else
76        config%i_band_from_g_sw &
77             &  = config%gas_optics_sw%spectral_def%i_band_number
78        config%i_band_from_reordered_g_sw &
79             &  = config%gas_optics_sw%spectral_def%i_band_number
80      end if
81      config%i_g_from_reordered_g_sw      = [ (jj, jj = 1,config%n_g_sw) ]
82
83      if (config%iverbosesetup >= 2) then
84        call config%gas_optics_sw%print()
85      end if
86
87      ! Override solar spectral irradiance, if filename provided
88      if (config%use_spectral_solar_cycle) then
89        call config%gas_optics_sw%read_spectral_solar_cycle(config%ssi_file_name, &
90             &           config%iverbosesetup, config%use_updated_solar_spectrum)
91      end if
92
93    end if
94
95    if (config%do_lw) then
96
97      ! Read longwave ecCKD gas optics NetCDF file
98      call config%gas_optics_lw%read(trim(config%gas_optics_lw_file_name), &
99           &                         config%iverbosesetup)
100
101      ! Copy over relevant properties
102      config%n_g_lw     = config%gas_optics_lw%ng
103
104      if (config%do_cloud_aerosol_per_lw_g_point) then
105        ! Bands and g points are the same
106        config%n_bands_lw = config%n_g_lw
107      else
108        ! Bands are groups of g points and span a continuous region of
109        ! wavenumber space
110        config%n_bands_lw = config%gas_optics_lw%spectral_def%nband
111      end if
112
113!      if (allocated(config%i_band_from_g_lw)) deallocate(config%i_band_from_g_lw)
114!      allocate(config%i_band_from_g_lw          (config%n_g_lw))
115!      if (allocated(config%i_band_from_reordered_g_lw)) deallocate(config%i_band_from_reordered_g_lw)
116!      allocate(config%i_band_from_reordered_g_lw(config%n_g_lw))
117!      if (allocated(config%i_g_from_reordered_g_lw)) deallocate(config%i_g_from_reordered_g_lw)
118!      allocate(config%i_g_from_reordered_g_lw   (config%n_g_lw))
119      if (.not.allocated(config%i_band_from_g_lw)) &
120         allocate(config%i_band_from_g_lw          (config%n_g_lw))
121      if (.not.allocated(config%i_band_from_reordered_g_lw)) &
122         allocate(config%i_band_from_reordered_g_lw(config%n_g_lw))
123      if (.not.allocated(config%i_g_from_reordered_g_lw)) &
124         allocate(config%i_g_from_reordered_g_lw   (config%n_g_lw))
125
126
127      if (config%do_cloud_aerosol_per_lw_g_point) then
128        config%i_band_from_g_lw           = [ (jj, jj = 1,config%n_g_lw) ]
129        config%i_band_from_reordered_g_lw = [ (jj, jj = 1,config%n_g_lw) ]
130      else
131        config%i_band_from_g_lw &
132             &  = config%gas_optics_lw%spectral_def%i_band_number
133        config%i_band_from_reordered_g_lw &
134             &  = config%gas_optics_lw%spectral_def%i_band_number
135      end if
136      config%i_g_from_reordered_g_lw      = [ (jj, jj = 1,config%n_g_lw) ]
137
138      if (config%iverbosesetup >= 2) then
139        call config%gas_optics_lw%print()
140      end if
141
142    end if
143
144    ! The i_spec_* variables are used solely for storing spectral
145    ! data, and this can either be by band or by g-point
146    if (config%do_save_spectral_flux) then
147      if (config%do_save_gpoint_flux) then
148        config%n_spec_sw = config%n_g_sw
149        config%n_spec_lw = config%n_g_lw
150        config%i_spec_from_reordered_g_sw => config%i_g_from_reordered_g_sw
151        config%i_spec_from_reordered_g_lw => config%i_g_from_reordered_g_lw
152      else
153        config%n_spec_sw = config%n_bands_sw
154        config%n_spec_lw = config%n_bands_lw
155        config%i_spec_from_reordered_g_sw => config%i_band_from_reordered_g_sw
156        config%i_spec_from_reordered_g_lw => config%i_band_from_reordered_g_lw
157      end if
158    else
159      config%n_spec_sw = 0
160      config%n_spec_lw = 0
161      nullify(config%i_spec_from_reordered_g_sw)
162      nullify(config%i_spec_from_reordered_g_lw)
163    end if
164
165    if (lhook) call dr_hook('radiation_ecckd_interface:setup_gas_optics',1,hook_handle)
166   
167  end subroutine setup_gas_optics
168
169
170  !---------------------------------------------------------------------
171  ! Scale gas mixing ratios according to required units
172  subroutine set_gas_units(gas)
173
174    use radiation_gas, only : gas_type, IVolumeMixingRatio
175    use yomhook,       only : lhook, dr_hook, jphook
176   
177    type(gas_type),    intent(inout) :: gas
178
179    real(jphook) :: hook_handle
180
181    if (lhook) call dr_hook('radiation_ecckd_interface:set_gas_units',0,hook_handle)
182
183    call gas%set_units(IVolumeMixingRatio)
184
185    if (lhook) call dr_hook('radiation_ecckd_interface:set_gas_units',1,hook_handle)
186
187  end subroutine set_gas_units
188
189
190  !---------------------------------------------------------------------
191  ! Compute gas optical depths, shortwave scattering, Planck function
192  ! and incoming shortwave radiation at top-of-atmosphere
193  subroutine gas_optics(ncol,nlev,istartcol,iendcol, &
194       &  config, single_level, thermodynamics, gas, &
195       &  od_lw, od_sw, ssa_sw, lw_albedo, planck_hl, lw_emission, &
196       &  incoming_sw)
197
198    use parkind1, only : jprb
199    use yomhook,  only : lhook, dr_hook, jphook
200
201    use radiation_config,         only : config_type
202    use radiation_thermodynamics, only : thermodynamics_type
203    use radiation_single_level,   only : single_level_type
204    use radiation_gas_constants,  only : NMaxGases
205    use radiation_gas
206
207    integer,                  intent(in) :: ncol               ! number of columns
208    integer,                  intent(in) :: nlev               ! number of levels
209    integer,                  intent(in) :: istartcol, iendcol ! range of cols to process
210    type(config_type),        intent(in) :: config
211    type(single_level_type),  intent(in) :: single_level
212    type(thermodynamics_type),intent(in) :: thermodynamics
213    type(gas_type),           intent(in) :: gas
214
215    ! Longwave albedo of the surface
216    real(jprb), dimension(config%n_g_lw,istartcol:iendcol), &
217         &  intent(in), optional :: lw_albedo
218
219    ! Gaseous layer optical depth in longwave and shortwave, and
220    ! shortwave single scattering albedo (i.e. fraction of extinction
221    ! due to Rayleigh scattering) at each g-point
222    real(jprb), dimension(config%n_g_lw,nlev,istartcol:iendcol), intent(out) :: &
223         &   od_lw
224    real(jprb), dimension(config%n_g_sw,nlev,istartcol:iendcol), intent(out) :: &
225         &   od_sw, ssa_sw
226
227    ! The Planck function (emitted flux from a black body) at half
228    ! levels at each longwave g-point
229    real(jprb), dimension(config%n_g_lw,nlev+1,istartcol:iendcol), &
230         &   intent(out), optional :: planck_hl
231    ! Planck function for the surface (W m-2)
232    real(jprb), dimension(config%n_g_lw,istartcol:iendcol), &
233         &   intent(out), optional :: lw_emission
234
235    ! The incoming shortwave flux into a plane perpendicular to the
236    ! incoming radiation at top-of-atmosphere in each of the shortwave
237    ! g-points
238    real(jprb), dimension(config%n_g_sw,istartcol:iendcol), &
239         &   intent(out), optional :: incoming_sw
240
241    ! Temperature at full levels (K)
242    real(jprb) :: temperature_fl(istartcol:iendcol,nlev)
243
244    integer :: jcol
245
246    real(jphook) :: hook_handle
247
248    if (lhook) call dr_hook('radiation_ecckd_interface:gas_optics',0,hook_handle)
249
250    !temperature_fl(istartcol:iendcol,:) &
251    !     &  = 0.5_jprb * (thermodynamics%temperature_hl(istartcol:iendcol,1:nlev) &
252    !     &               +thermodynamics%temperature_hl(istartcol:iendcol,2:nlev+1))
253 
254    temperature_fl(istartcol:iendcol,:) &
255         &  = (thermodynamics%temperature_hl(istartcol:iendcol,1:nlev) &
256         &     *thermodynamics%pressure_hl(istartcol:iendcol,1:nlev) &
257         &    +thermodynamics%temperature_hl(istartcol:iendcol,2:nlev+1) &
258         &     *thermodynamics%pressure_hl(istartcol:iendcol,2:nlev+1)) &
259         &  / (thermodynamics%pressure_hl(istartcol:iendcol,1:nlev) &
260         &    +thermodynamics%pressure_hl(istartcol:iendcol,2:nlev+1))
261 
262    if (config%do_sw) then
263
264      call config%gas_optics_sw%calc_optical_depth(ncol,nlev,istartcol,iendcol, &
265           &  NMaxGases, thermodynamics%pressure_hl, &
266           &  temperature_fl, &
267           &  gas%mixing_ratio, &
268!           &  reshape(gas%mixing_ratio(istartcol:iendcol,:,:), &
269!           &          [nlev,iendcol-istartcol+1,NMaxGases],order=[2,1,3]), &
270           &  od_sw, rayleigh_od_fl=ssa_sw)
271      ! At this point od_sw = absorption optical depth and ssa_sw =
272      ! rayleigh optical depth: convert to total optical depth and
273      ! single-scattering albedo
274      od_sw = od_sw + ssa_sw
275      ssa_sw = ssa_sw / od_sw
276
277      if (present(incoming_sw)) then
278        if (single_level%spectral_solar_cycle_multiplier == 0.0_jprb) then
279          call config%gas_optics_sw%calc_incoming_sw(single_level%solar_irradiance, &
280               &        incoming_sw)
281        else
282          call config%gas_optics_sw%calc_incoming_sw(single_level%solar_irradiance, &
283               &        incoming_sw, single_level%spectral_solar_cycle_multiplier)
284        end if
285      end if
286
287    end if
288
289    if (config%do_lw) then
290
291      call config%gas_optics_lw%calc_optical_depth(ncol,nlev,istartcol,iendcol, &
292           &  NMaxGases, thermodynamics%pressure_hl, &
293           &  temperature_fl, &
294           &  gas%mixing_ratio, &
295!           &  reshape(gas%mixing_ratio(istartcol:iendcol,:,:), &
296!           &          [nlev,iendcol-istartcol+1,NMaxGases],order=[2,1,3]), &
297           &  od_lw)
298
299      ! Calculate the Planck function for each g point
300      do jcol = istartcol,iendcol
301        call config%gas_optics_lw%calc_planck_function(nlev+1, &
302             &  thermodynamics%temperature_hl(jcol,:), planck_hl(:,:,jcol))
303      end do
304      call config%gas_optics_lw%calc_planck_function(iendcol+1-istartcol, &
305           &  single_level%skin_temperature(istartcol:iendcol), &
306           &  lw_emission(:,:))
307      lw_emission = lw_emission * (1.0_jprb - lw_albedo)
308
309    end if
310
311    if (lhook) call dr_hook('radiation_ecckd_interface:gas_optics',1,hook_handle)
312   
313  end subroutine gas_optics
314
315  ! !---------------------------------------------------------------------
316  ! ! Externally facing function for computing the Planck function
317  ! ! without reference to any gas profile; typically this would be used
318  ! ! for computing the emission by a surface.
319  ! subroutine planck_function(config, temperature, planck_surf)
320
321  !   use parkind1,                 only : jprb
322  !   use radiation_config,         only : config_type
323
324  !   type(config_type), intent(in) :: config
325  !   real(jprb),        intent(in) :: temperature
326
327  !   ! Planck function of the surface (W m-2)
328  !   real(jprb), dimension(config%n_g_lw), intent(out) :: planck_surf
329
330  ! end subroutine planck_function
331
332end module radiation_ecckd_interface
Note: See TracBrowser for help on using the repository browser.