Line data Source code
1 : module radconstants
2 :
3 : ! This module contains constants that are specific to the radiative transfer
4 : ! code used in the RRTMGP model.
5 :
6 : use shr_kind_mod, only: r8 => shr_kind_r8
7 : use mo_gas_optics_rrtmgp, only: ty_gas_optics_rrtmgp
8 : use cam_abortutils, only: endrun
9 :
10 : implicit none
11 : private
12 : save
13 :
14 : ! Number of bands in SW and LW. These values must match data in the RRTMGP coefficients datasets.
15 : ! But they are needed to allocate space in the physics buffer and need to be available before the
16 : ! RRTMGP datasets are read. So they are set as parameters here and checked in the
17 : ! set_wavenumber_bands subroutine after the datasets are read.
18 : integer, parameter, public :: nswbands = 14
19 : integer, parameter, public :: nlwbands = 16
20 :
21 : ! Band limits (set from data in RRTMGP coefficient datasets)
22 : real(r8), target :: wavenumber_low_shortwave(nswbands)
23 : real(r8), target :: wavenumber_high_shortwave(nswbands)
24 : real(r8), target :: wavenumber_low_longwave(nlwbands)
25 : real(r8), target :: wavenumber_high_longwave(nlwbands)
26 :
27 : logical :: wavenumber_boundaries_set = .false.
28 :
29 : ! First and last g-point for each band.
30 : integer, public, protected :: band2gpt_sw(2,nswbands)
31 :
32 : integer, public, protected :: nswgpts ! number of SW g-points
33 : integer, public, protected :: nlwgpts ! number of LW g-points
34 :
35 : ! These are indices to specific bands for diagnostic output and COSP input.
36 : integer, public, protected :: idx_sw_diag = -1 ! band contains 500-nm wave
37 : integer, public, protected :: idx_nir_diag = -1 ! band contains 1000-nm wave
38 : integer, public, protected :: idx_uv_diag = -1 ! band contains 400-nm wave
39 : integer, public, protected :: idx_lw_diag = -1 ! band contains 1000 cm-1 wave (H20 window)
40 : integer, public, protected :: idx_sw_cloudsim = -1 ! band contains 670-nm wave (for COSP)
41 : integer, public, protected :: idx_lw_cloudsim = -1 ! band contains 10.5 micron wave (for COSP)
42 :
43 : ! GASES TREATED BY RADIATION (line spectra)
44 : ! These names are recognized by RRTMGP. They are in the coefficients files as
45 : ! lower case strings. These upper case names are used by CAM's namelist and
46 : ! rad_constituents module.
47 : integer, public, parameter :: gasnamelength = 5
48 : integer, public, parameter :: nradgas = 8
49 : character(len=gasnamelength), public, parameter :: gaslist(nradgas) &
50 : = (/'H2O ','O3 ', 'O2 ', 'CO2 ', 'N2O ', 'CH4 ', 'CFC11', 'CFC12'/)
51 :
52 : ! what is the minimum mass mixing ratio that can be supported by radiation implementation?
53 : real(r8), public, parameter :: minmmr(nradgas) = epsilon(1._r8)
54 :
55 : public :: &
56 : set_wavenumber_bands, &
57 : get_sw_spectral_boundaries, &
58 : get_lw_spectral_boundaries, &
59 : get_band_index_by_value, &
60 : rad_gas_index
61 :
62 : !=========================================================================================
63 : contains
64 : !=========================================================================================
65 :
66 2304 : subroutine set_wavenumber_bands(kdist_sw, kdist_lw)
67 :
68 : ! Set the low and high limits of the wavenumber grid for sw and lw.
69 : ! Values come from RRTMGP coefficients datasets, and are stored in the
70 : ! kdist objects.
71 : !
72 : ! Set band indices for bands containing specific wavelengths.
73 :
74 : ! Arguments
75 : type(ty_gas_optics_rrtmgp), intent(in) :: kdist_sw
76 : type(ty_gas_optics_rrtmgp), intent(in) :: kdist_lw
77 :
78 : ! Local variables
79 : integer :: istat
80 2304 : real(r8), allocatable :: values(:,:)
81 :
82 : character(len=128) :: errmsg
83 : character(len=*), parameter :: sub = 'set_wavenumber_bands'
84 : !----------------------------------------------------------------------------
85 :
86 : ! Check that number of sw/lw bands in gas optics files matches the parameters.
87 2304 : if (kdist_sw%get_nband() /= nswbands) then
88 0 : write(errmsg,'(a,i4,a,i4)') 'number of sw bands in file, ', kdist_sw%get_nband(), &
89 0 : ", doesn't match parameter nswbands= ", nswbands
90 0 : call endrun(sub//': ERROR: '//trim(errmsg))
91 : end if
92 2304 : if (kdist_lw%get_nband() /= nlwbands) then
93 0 : write(errmsg,'(a,i4,a,i4)') 'number of lw bands in file, ', kdist_lw%get_nband(), &
94 0 : ", doesn't match parameter nlwbands= ", nlwbands
95 0 : call endrun(sub//': ERROR: '//trim(errmsg))
96 : end if
97 :
98 2304 : nswgpts = kdist_sw%get_ngpt()
99 2304 : nlwgpts = kdist_lw%get_ngpt()
100 :
101 : ! SW band bounds in cm^-1
102 2304 : allocate( values(2,nswbands), stat=istat )
103 2304 : if (istat/=0) then
104 0 : call endrun(sub//': ERROR allocating array: values(2,nswbands)')
105 : end if
106 2304 : values = kdist_sw%get_band_lims_wavenumber()
107 34560 : wavenumber_low_shortwave = values(1,:)
108 34560 : wavenumber_high_shortwave = values(2,:)
109 :
110 : ! First and last g-point for each SW band:
111 99072 : band2gpt_sw = kdist_sw%get_band_lims_gpoint()
112 :
113 : ! Indices into specific bands
114 2304 : idx_sw_diag = get_band_index_by_value('sw', 500.0_r8, 'nm')
115 2304 : idx_nir_diag = get_band_index_by_value('sw', 1000.0_r8, 'nm')
116 2304 : idx_uv_diag = get_band_index_by_value('sw', 400._r8, 'nm')
117 2304 : idx_sw_cloudsim = get_band_index_by_value('sw', 0.67_r8, 'micron')
118 :
119 2304 : deallocate(values)
120 :
121 : ! LW band bounds in cm^-1
122 2304 : allocate( values(2,nlwbands), stat=istat )
123 2304 : if (istat/=0) then
124 0 : call endrun(sub//': ERROR allocating array: values(2,nlwbands)')
125 : end if
126 2304 : values = kdist_lw%get_band_lims_wavenumber()
127 39168 : wavenumber_low_longwave = values(1,:)
128 39168 : wavenumber_high_longwave = values(2,:)
129 :
130 : ! Indices into specific bands
131 2304 : idx_lw_diag = get_band_index_by_value('lw', 1000.0_r8, 'cm^-1')
132 2304 : idx_lw_cloudsim = get_band_index_by_value('lw', 10.5_r8, 'micron')
133 :
134 2304 : wavenumber_boundaries_set = .true.
135 :
136 2304 : end subroutine set_wavenumber_bands
137 :
138 : !=========================================================================================
139 :
140 4608 : subroutine get_sw_spectral_boundaries(low_boundaries, high_boundaries, units)
141 :
142 : ! provide spectral boundaries of each shortwave band
143 :
144 : real(r8), intent(out) :: low_boundaries(nswbands), high_boundaries(nswbands)
145 : character(*), intent(in) :: units ! requested units
146 :
147 : character(len=*), parameter :: sub = 'get_sw_spectral_boundaries'
148 : !----------------------------------------------------------------------------
149 :
150 4608 : if (.not. wavenumber_boundaries_set) then
151 0 : call endrun(sub//': ERROR, wavenumber boundaries not set. ')
152 : end if
153 :
154 2304 : select case (units)
155 : case ('inv_cm','cm^-1','cm-1')
156 2304 : low_boundaries = wavenumber_low_shortwave
157 2304 : high_boundaries = wavenumber_high_shortwave
158 : case('m','meter','meters')
159 0 : low_boundaries = 1.e-2_r8/wavenumber_high_shortwave
160 0 : high_boundaries = 1.e-2_r8/wavenumber_low_shortwave
161 : case('nm','nanometer','nanometers')
162 34560 : low_boundaries = 1.e7_r8/wavenumber_high_shortwave
163 34560 : high_boundaries = 1.e7_r8/wavenumber_low_shortwave
164 : case('um','micrometer','micrometers','micron','microns')
165 0 : low_boundaries = 1.e4_r8/wavenumber_high_shortwave
166 0 : high_boundaries = 1.e4_r8/wavenumber_low_shortwave
167 : case('cm','centimeter','centimeters')
168 0 : low_boundaries = 1._r8/wavenumber_high_shortwave
169 0 : high_boundaries = 1._r8/wavenumber_low_shortwave
170 : case default
171 4608 : call endrun(sub//': ERROR, requested spectral units not recognized: '//units)
172 : end select
173 :
174 4608 : end subroutine get_sw_spectral_boundaries
175 :
176 : !=========================================================================================
177 :
178 2304 : subroutine get_lw_spectral_boundaries(low_boundaries, high_boundaries, units)
179 :
180 : ! provide spectral boundaries of each longwave band
181 :
182 : real(r8), intent(out) :: low_boundaries(nlwbands), high_boundaries(nlwbands)
183 : character(*), intent(in) :: units ! requested units
184 :
185 : character(len=*), parameter :: sub = 'get_lw_spectral_boundaries'
186 : !----------------------------------------------------------------------------
187 :
188 2304 : if (.not. wavenumber_boundaries_set) then
189 0 : call endrun(sub//': ERROR, wavenumber boundaries not set. ')
190 : end if
191 :
192 0 : select case (units)
193 : case ('inv_cm','cm^-1','cm-1')
194 0 : low_boundaries = wavenumber_low_longwave
195 0 : high_boundaries = wavenumber_high_longwave
196 : case('m','meter','meters')
197 0 : low_boundaries = 1.e-2_r8/wavenumber_high_longwave
198 0 : high_boundaries = 1.e-2_r8/wavenumber_low_longwave
199 : case('nm','nanometer','nanometers')
200 0 : low_boundaries = 1.e7_r8/wavenumber_high_longwave
201 0 : high_boundaries = 1.e7_r8/wavenumber_low_longwave
202 : case('um','micrometer','micrometers','micron','microns')
203 39168 : low_boundaries = 1.e4_r8/wavenumber_high_longwave
204 39168 : high_boundaries = 1.e4_r8/wavenumber_low_longwave
205 : case('cm','centimeter','centimeters')
206 0 : low_boundaries = 1._r8/wavenumber_high_longwave
207 0 : high_boundaries = 1._r8/wavenumber_low_longwave
208 : case default
209 2304 : call endrun(sub//': ERROR, requested spectral units not recognized: '//units)
210 : end select
211 :
212 2304 : end subroutine get_lw_spectral_boundaries
213 :
214 : !=========================================================================================
215 :
216 620752 : integer function rad_gas_index(gasname)
217 :
218 : ! return the index in the gaslist array of the specified gasname
219 :
220 : character(len=*),intent(in) :: gasname
221 : integer :: igas
222 :
223 620752 : rad_gas_index = -1
224 2793384 : do igas = 1, nradgas
225 2793384 : if (trim(gaslist(igas)).eq.trim(gasname)) then
226 620752 : rad_gas_index = igas
227 620752 : return
228 : endif
229 : enddo
230 0 : call endrun ("rad_gas_index: can not find gas with name "//gasname)
231 0 : end function rad_gas_index
232 :
233 : !=========================================================================================
234 :
235 13824 : function get_band_index_by_value(swlw, targetvalue, units) result(ans)
236 :
237 : ! Find band index for requested wavelength/wavenumber.
238 :
239 : character(len=*), intent(in) :: swlw ! sw or lw bands
240 : real(r8), intent(in) :: targetvalue
241 : character(len=*), intent(in) :: units ! units of targetvalue
242 : integer :: ans
243 :
244 : ! local
245 13824 : real(r8), pointer, dimension(:) :: lowboundaries, highboundaries
246 : real(r8) :: tgt
247 : integer :: nbnds, i
248 :
249 : character(len=128) :: errmsg
250 : character(len=*), parameter :: sub = 'get_band_index_by_value'
251 : !----------------------------------------------------------------------------
252 :
253 9216 : select case (swlw)
254 : case ('sw','SW','shortwave')
255 9216 : nbnds = nswbands
256 9216 : lowboundaries => wavenumber_low_shortwave
257 9216 : highboundaries => wavenumber_high_shortwave
258 : case ('lw', 'LW', 'longwave')
259 4608 : nbnds = nlwbands
260 4608 : lowboundaries => wavenumber_low_longwave
261 4608 : highboundaries => wavenumber_high_longwave
262 : case default
263 0 : call endrun('radconstants.F90: get_band_index_by_value: type of bands not recognized: '//swlw)
264 : end select
265 :
266 : ! band info is in cm^-1 but target value may be other units,
267 : ! so convert targetvalue to cm^-1
268 2304 : select case (units)
269 : case ('inv_cm','cm^-1','cm-1')
270 2304 : tgt = targetvalue
271 : case('m','meter','meters')
272 0 : tgt = 1.0_r8 / (targetvalue * 1.e2_r8)
273 : case('nm','nanometer','nanometers')
274 6912 : tgt = 1.0_r8 / (targetvalue * 1.e-7_r8)
275 : case('um','micrometer','micrometers','micron','microns')
276 4608 : tgt = 1.0_r8 / (targetvalue * 1.e-4_r8)
277 : case('cm','centimeter','centimeters')
278 0 : tgt = 1._r8/targetvalue
279 : case default
280 13824 : call endrun('radconstants.F90: get_band_index_by_value: units not recognized: '//units)
281 : end select
282 :
283 : ! now just loop through the array
284 13824 : ans = 0
285 126720 : do i = 1,nbnds
286 126720 : if ((tgt > lowboundaries(i)) .and. (tgt <= highboundaries(i))) then
287 : ans = i
288 : exit
289 : end if
290 : end do
291 :
292 13824 : if (ans == 0) then
293 0 : write(errmsg,'(f10.3,a,a)') targetvalue, ' ', trim(units)
294 0 : call endrun(sub//': band not found containing wave: '//trim(errmsg))
295 : end if
296 :
297 13824 : end function get_band_index_by_value
298 :
299 : !=========================================================================================
300 :
301 : end module radconstants
|
