source: LMDZ6/branches/Test_modipsl/libf/phylmd/ecrad/radiation_ecckd_interface.F90 @ 5425

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

Merge LMDZ_ECRad branch back into trunk!

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