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

Last change on this file since 5418 was 4853, checked in by idelkadi, 9 months ago

Ecrad update in LMDZ / Implementation of Ecrad double call in LMDZ

  • version 1.6.1 (commit 210d7911380f53a788c3cad73b3cf9b4e022ef87)
  • interface routines between lmdz and ecrad grouped in a new "lmdz" directory
  • double call method redesigned so as to go through the Ecrad initialization and configuration part only once for the entire simulation
  • clean-up in the read.F routine: delete unitules arguments
  • modification of default gas model in namelist (default: ECCKD)
File size: 13.2 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 .and. config%i_gas_model_sw == IGasModelECCKD) 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      allocate(config%i_band_from_g_sw          (config%n_g_sw))
60      allocate(config%i_band_from_reordered_g_sw(config%n_g_sw))
61      allocate(config%i_g_from_reordered_g_sw   (config%n_g_sw))
62       
63      if (config%do_cloud_aerosol_per_sw_g_point) then
64        config%i_band_from_g_sw           = [ (jj, jj = 1,config%n_g_sw) ]
65        config%i_band_from_reordered_g_sw = [ (jj, jj = 1,config%n_g_sw) ]
66      else
67        config%i_band_from_g_sw &
68             &  = config%gas_optics_sw%spectral_def%i_band_number
69        config%i_band_from_reordered_g_sw &
70             &  = config%gas_optics_sw%spectral_def%i_band_number
71      end if
72      config%i_g_from_reordered_g_sw      = [ (jj, jj = 1,config%n_g_sw) ]
73
74      if (config%iverbosesetup >= 2) then
75        call config%gas_optics_sw%print()
76      end if
77
78      ! Override solar spectral irradiance, if filename provided
79      if (config%use_spectral_solar_cycle) then
80        call config%gas_optics_sw%read_spectral_solar_cycle(config%ssi_file_name, &
81             &           config%iverbosesetup, config%use_updated_solar_spectrum)
82      end if
83
84    end if
85
86    if (config%do_lw .and. config%i_gas_model_lw == IGasModelECCKD) then
87
88      ! Read longwave ecCKD gas optics NetCDF file
89      call config%gas_optics_lw%read(trim(config%gas_optics_lw_file_name), &
90           &                         config%iverbosesetup)
91
92      ! Copy over relevant properties
93      config%n_g_lw     = config%gas_optics_lw%ng
94
95      if (config%do_cloud_aerosol_per_lw_g_point) then
96        ! Bands and g points are the same
97        config%n_bands_lw = config%n_g_lw
98      else
99        ! Bands are groups of g points and span a continuous region of
100        ! wavenumber space
101        config%n_bands_lw = config%gas_optics_lw%spectral_def%nband
102      end if
103
104      allocate(config%i_band_from_g_lw          (config%n_g_lw))
105      allocate(config%i_band_from_reordered_g_lw(config%n_g_lw))
106      allocate(config%i_g_from_reordered_g_lw   (config%n_g_lw))
107
108      if (config%do_cloud_aerosol_per_lw_g_point) then
109        config%i_band_from_g_lw           = [ (jj, jj = 1,config%n_g_lw) ]
110        config%i_band_from_reordered_g_lw = [ (jj, jj = 1,config%n_g_lw) ]
111      else
112        config%i_band_from_g_lw &
113             &  = config%gas_optics_lw%spectral_def%i_band_number
114        config%i_band_from_reordered_g_lw &
115             &  = config%gas_optics_lw%spectral_def%i_band_number
116      end if
117      config%i_g_from_reordered_g_lw      = [ (jj, jj = 1,config%n_g_lw) ]
118
119      if (config%iverbosesetup >= 2) then
120        call config%gas_optics_lw%print()
121      end if
122
123    end if
124
125    ! The i_spec_* variables are used solely for storing spectral
126    ! data, and this can either be by band or by g-point
127    if (config%do_save_spectral_flux) then
128      if (config%do_save_gpoint_flux) then
129        config%n_spec_sw = config%n_g_sw
130        config%n_spec_lw = config%n_g_lw
131        config%i_spec_from_reordered_g_sw => config%i_g_from_reordered_g_sw
132        config%i_spec_from_reordered_g_lw => config%i_g_from_reordered_g_lw
133      else
134        config%n_spec_sw = config%n_bands_sw
135        config%n_spec_lw = config%n_bands_lw
136        config%i_spec_from_reordered_g_sw => config%i_band_from_reordered_g_sw
137        config%i_spec_from_reordered_g_lw => config%i_band_from_reordered_g_lw
138      end if
139    else
140      config%n_spec_sw = 0
141      config%n_spec_lw = 0
142      nullify(config%i_spec_from_reordered_g_sw)
143      nullify(config%i_spec_from_reordered_g_lw)
144    end if
145
146    if (lhook) call dr_hook('radiation_ecckd_interface:setup_gas_optics',1,hook_handle)
147   
148  end subroutine setup_gas_optics
149
150
151  !---------------------------------------------------------------------
152  ! Scale gas mixing ratios according to required units
153  subroutine set_gas_units(gas)
154
155    use radiation_gas, only : gas_type, IVolumeMixingRatio
156    use yomhook,       only : lhook, dr_hook, jphook
157   
158    type(gas_type),    intent(inout) :: gas
159
160    real(jphook) :: hook_handle
161
162    if (lhook) call dr_hook('radiation_ecckd_interface:set_gas_units',0,hook_handle)
163
164    call gas%set_units(IVolumeMixingRatio)
165
166    if (lhook) call dr_hook('radiation_ecckd_interface:set_gas_units',1,hook_handle)
167
168  end subroutine set_gas_units
169
170
171  !---------------------------------------------------------------------
172  ! Compute gas optical depths, shortwave scattering, Planck function
173  ! and incoming shortwave radiation at top-of-atmosphere
174  subroutine gas_optics(ncol,nlev,istartcol,iendcol, &
175       &  config, single_level, thermodynamics, gas, &
176       &  od_lw, od_sw, ssa_sw, lw_albedo, planck_hl, lw_emission, &
177       &  incoming_sw)
178
179    use parkind1, only : jprb
180    use yomhook,  only : lhook, dr_hook, jphook
181
182    use radiation_config,         only : config_type, IGasModelECCKD
183    use radiation_thermodynamics, only : thermodynamics_type
184    use radiation_single_level,   only : single_level_type
185    use radiation_gas_constants,  only : NMaxGases
186    use radiation_gas
187
188    integer,                  intent(in) :: ncol               ! number of columns
189    integer,                  intent(in) :: nlev               ! number of levels
190    integer,                  intent(in) :: istartcol, iendcol ! range of cols to process
191    type(config_type),        intent(in) :: config
192    type(single_level_type),  intent(in) :: single_level
193    type(thermodynamics_type),intent(in) :: thermodynamics
194    type(gas_type),           intent(in) :: gas
195
196    ! Longwave albedo of the surface
197    real(jprb), dimension(config%n_g_lw,istartcol:iendcol), &
198         &  intent(in), optional :: lw_albedo
199
200    ! Gaseous layer optical depth in longwave and shortwave, and
201    ! shortwave single scattering albedo (i.e. fraction of extinction
202    ! due to Rayleigh scattering) at each g-point
203    real(jprb), dimension(config%n_g_lw,nlev,istartcol:iendcol), intent(out) :: &
204         &   od_lw
205    real(jprb), dimension(config%n_g_sw,nlev,istartcol:iendcol), intent(out) :: &
206         &   od_sw, ssa_sw
207
208    ! The Planck function (emitted flux from a black body) at half
209    ! levels at each longwave g-point
210    real(jprb), dimension(config%n_g_lw,nlev+1,istartcol:iendcol), &
211         &   intent(out), optional :: planck_hl
212    ! Planck function for the surface (W m-2)
213    real(jprb), dimension(config%n_g_lw,istartcol:iendcol), &
214         &   intent(out), optional :: lw_emission
215
216    ! The incoming shortwave flux into a plane perpendicular to the
217    ! incoming radiation at top-of-atmosphere in each of the shortwave
218    ! g-points
219    real(jprb), dimension(config%n_g_sw,istartcol:iendcol), &
220         &   intent(out), optional :: incoming_sw
221
222    ! Temperature at full levels (K)
223    real(jprb) :: temperature_fl(istartcol:iendcol,nlev)
224
225    real(jprb) :: concentration_scaling(NMaxGases)
226   
227    logical :: is_volume_mixing_ratio
228   
229    integer :: jcol, jlev, jg
230
231    real(jphook) :: hook_handle
232
233    if (lhook) call dr_hook('radiation_ecckd_interface:gas_optics',0,hook_handle)
234
235    !temperature_fl(istartcol:iendcol,:) &
236    !     &  = 0.5_jprb * (thermodynamics%temperature_hl(istartcol:iendcol,1:nlev) &
237    !     &               +thermodynamics%temperature_hl(istartcol:iendcol,2:nlev+1))
238 
239    temperature_fl(istartcol:iendcol,:) &
240         &  = (thermodynamics%temperature_hl(istartcol:iendcol,1:nlev) &
241         &     *thermodynamics%pressure_hl(istartcol:iendcol,1:nlev) &
242         &    +thermodynamics%temperature_hl(istartcol:iendcol,2:nlev+1) &
243         &     *thermodynamics%pressure_hl(istartcol:iendcol,2:nlev+1)) &
244         &  / (thermodynamics%pressure_hl(istartcol:iendcol,1:nlev) &
245         &    +thermodynamics%pressure_hl(istartcol:iendcol,2:nlev+1))
246
247    ! Check that the gas concentrations are stored in volume mixing
248    ! ratio with no scaling; if not, return a vector of scalings
249    call gas%assert_units(IVolumeMixingRatio, scale_factor=1.0_jprb, &
250         &                istatus=is_volume_mixing_ratio)
251    if (.not. is_volume_mixing_ratio) then
252      call gas%get_scaling(IVolumeMixingRatio, concentration_scaling)
253    else
254      concentration_scaling = 1.0_jprb
255    end if
256   
257    if (config%do_sw .and. config%i_gas_model_sw == IGasModelECCKD) then
258
259      if (is_volume_mixing_ratio) then
260        call config%gas_optics_sw%calc_optical_depth(ncol,nlev,istartcol,iendcol, &
261             &  NMaxGases, thermodynamics%pressure_hl, &
262             &  temperature_fl, gas%mixing_ratio, &
263             &  od_sw, rayleigh_od_fl=ssa_sw)
264      else
265        call config%gas_optics_sw%calc_optical_depth(ncol,nlev,istartcol,iendcol, &
266             &  NMaxGases, thermodynamics%pressure_hl, &
267             &  temperature_fl, gas%mixing_ratio, &
268             &  od_sw, rayleigh_od_fl=ssa_sw, concentration_scaling=concentration_scaling)
269      end if
270     
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      do jcol = istartcol,iendcol
275        do jlev = 1, nlev
276          do jg = 1, config%n_g_sw
277            od_sw(jg,jlev,jcol)  = od_sw(jg,jlev,jcol) + ssa_sw(jg,jlev,jcol)
278            ssa_sw(jg,jlev,jcol) = ssa_sw(jg,jlev,jcol) / od_sw(jg,jlev,jcol)
279          end do
280        end do
281      end do
282
283      if (present(incoming_sw)) then
284        if (single_level%spectral_solar_cycle_multiplier == 0.0_jprb) then
285          call config%gas_optics_sw%calc_incoming_sw(single_level%solar_irradiance, &
286               &        incoming_sw)
287        else
288          call config%gas_optics_sw%calc_incoming_sw(single_level%solar_irradiance, &
289               &        incoming_sw, single_level%spectral_solar_cycle_multiplier)
290        end if
291      end if
292
293    end if
294
295    if (config%do_lw .and. config%i_gas_model_lw == IGasModelECCKD) then
296
297      if (is_volume_mixing_ratio) then
298        call config%gas_optics_lw%calc_optical_depth(ncol,nlev,istartcol,iendcol, &
299             &  NMaxGases, thermodynamics%pressure_hl, &
300             &  temperature_fl, gas%mixing_ratio, &
301             &  od_lw)
302      else
303        call config%gas_optics_lw%calc_optical_depth(ncol,nlev,istartcol,iendcol, &
304             &  NMaxGases, thermodynamics%pressure_hl, &
305             &  temperature_fl, gas%mixing_ratio, &
306             &  od_lw, concentration_scaling=concentration_scaling)
307      end if
308
309      ! Calculate the Planck function for each g point
310      do jcol = istartcol,iendcol
311        call config%gas_optics_lw%calc_planck_function(nlev+1, &
312             &  thermodynamics%temperature_hl(jcol,:), planck_hl(:,:,jcol))
313      end do
314      call config%gas_optics_lw%calc_planck_function(iendcol+1-istartcol, &
315           &  single_level%skin_temperature(istartcol:iendcol), &
316           &  lw_emission(:,:))
317!NEC$ forced_collapse
318      lw_emission = lw_emission * (1.0_jprb - lw_albedo)
319
320    end if
321
322    if (lhook) call dr_hook('radiation_ecckd_interface:gas_optics',1,hook_handle)
323   
324  end subroutine gas_optics
325
326  ! !---------------------------------------------------------------------
327  ! ! Externally facing function for computing the Planck function
328  ! ! without reference to any gas profile; typically this would be used
329  ! ! for computing the emission by a surface.
330  ! subroutine planck_function(config, temperature, planck_surf)
331
332  !   use parkind1,                 only : jprb
333  !   use radiation_config,         only : config_type
334
335  !   type(config_type), intent(in) :: config
336  !   real(jprb),        intent(in) :: temperature
337
338  !   ! Planck function of the surface (W m-2)
339  !   real(jprb), dimension(config%n_g_lw), intent(out) :: planck_surf
340
341  ! end subroutine planck_function
342
343end module radiation_ecckd_interface
Note: See TracBrowser for help on using the repository browser.