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