Line data Source code
1 : ! This code is part of
2 : ! RRTM for GCM Applications - Parallel (RRTMGP)
3 : !
4 : ! Eli Mlawer and Robert Pincus
5 : ! Andre Wehe and Jennifer Delamere
6 : ! email: rrtmgp@aer.com
7 : !
8 : ! Copyright 2015-, Atmospheric and Environmental Research,
9 : ! Regents of the University of Colorado, Trustees of Columbia University. All right reserved.
10 : !
11 : ! Use and duplication is permitted under the terms of the
12 : ! BSD 3-clause license, see http://opensource.org/licenses/BSD-3-Clause
13 : ! -------------------------------------------------------------------------------------------------
14 : !>
15 : !> ## Numeric calculations for gas optics. Absorption and Rayleigh optical depths, Planck source functions.
16 : !>
17 : !> - Interpolation coefficients are computed, then used in subsequent routines.
18 : !> - All applications will call compute_tau_absorption();
19 : !> compute_tau_rayleigh() and/or compute_Planck_source() will be called depending on the
20 : !> configuration of the k-distribution.
21 : !> - The details of the interpolation scheme are not particaulrly important as long as arrays including
22 : !> tables are passed consisently between kernels.
23 : !>
24 : ! -------------------------------------------------------------------------------------------------
25 :
26 : module mo_gas_optics_rrtmgp_kernels
27 : use mo_rte_kind, only : wp, wl
28 : use mo_rte_util_array,only : zero_array
29 : implicit none
30 : private
31 : public :: interpolation, compute_tau_absorption, compute_tau_rayleigh, compute_Planck_source
32 : contains
33 : ! --------------------------------------------------------------------------------------
34 : !> Compute interpolation coefficients
35 : !> for calculations of major optical depths, minor optical depths, Rayleigh,
36 : !> and Planck fractions
37 1135399 : subroutine interpolation( &
38 : ncol,nlay,ngas,nflav,neta, npres, ntemp, &
39 1135399 : flavor, &
40 1135399 : press_ref_log, temp_ref,press_ref_log_delta, &
41 : temp_ref_min,temp_ref_delta,press_ref_trop_log, &
42 1135399 : vmr_ref, &
43 1135399 : play,tlay,col_gas, &
44 1135399 : jtemp,fmajor,fminor,col_mix,tropo,jeta,jpress) bind(C, name="rrtmgp_interpolation")
45 : ! input dimensions
46 : integer, intent(in) :: ncol,nlay
47 : !! physical domain size
48 : integer, intent(in) :: ngas,nflav,neta,npres,ntemp
49 : !! k-distribution table dimensions
50 : integer, dimension(2,nflav), intent(in) :: flavor
51 : !! index into vmr_ref of major gases for each flavor
52 : real(wp), dimension(npres), intent(in) :: press_ref_log
53 : !! log of pressure dimension in RRTMGP tables
54 : real(wp), dimension(ntemp), intent(in) :: temp_ref
55 : !! temperature dimension in RRTMGP tables
56 : real(wp), intent(in) :: press_ref_log_delta, &
57 : temp_ref_min, temp_ref_delta, &
58 : press_ref_trop_log
59 : !! constants related to RRTMGP tables
60 : real(wp), dimension(2,0:ngas,ntemp), intent(in) :: vmr_ref
61 : !! reference volume mixing ratios used in compute "binary species parameter" eta
62 :
63 : ! inputs from profile or parent function
64 : real(wp), dimension(ncol,nlay), intent(in) :: play, tlay
65 : !! input pressure (Pa?) and temperature (K)
66 : real(wp), dimension(ncol,nlay,0:ngas), intent(in) :: col_gas
67 : !! input column gas amount - molecules/cm^2
68 : ! outputs
69 : integer, dimension(ncol,nlay), intent(out) :: jtemp, jpress
70 : !! temperature and pressure interpolation indexes
71 : logical(wl), dimension(ncol,nlay), intent(out) :: tropo
72 : !! use lower (or upper) atmosphere tables
73 : integer, dimension(2, ncol,nlay,nflav), intent(out) :: jeta
74 : !! Index for binary species interpolation
75 : #if !defined(__INTEL_LLVM_COMPILER) && __INTEL_COMPILER >= 1910
76 : ! A performance-hitting workaround for the vectorization problem reported in
77 : ! https://github.com/earth-system-radiation/rte-rrtmgp/issues/159
78 : ! The known affected compilers are Intel Fortran Compiler Classic
79 : ! 2021.4, 2021.5 and 2022.1. We do not limit the workaround to these
80 : ! versions because it is not clear when the compiler bug will be fixed, see
81 : ! https://community.intel.com/t5/Intel-Fortran-Compiler/Compiler-vectorization-bug/m-p/1362591.
82 : ! We, however, limit the workaround to the Classic versions only since the
83 : ! problem is not confirmed for the Intel Fortran Compiler oneAPI (a.k.a
84 : ! 'ifx'), which does not mean there is none though.
85 : real(wp), dimension(:, :, :, :), intent(out) :: col_mix
86 : #else
87 : real(wp), dimension(2, ncol,nlay,nflav), intent(out) :: col_mix
88 : !! combination of major species's column amounts (first index is strat/trop)
89 : #endif
90 : real(wp), dimension(2,2,2,ncol,nlay,nflav), intent(out) :: fmajor
91 : !! Interpolation weights in pressure, eta, strat/trop
92 : real(wp), dimension(2,2, ncol,nlay,nflav), intent(out) :: fminor
93 : !! Interpolation fraction in eta, strat/trop
94 : ! -----------------
95 : ! local
96 2270798 : real(wp), dimension(ncol,nlay) :: ftemp, fpress ! interpolation fraction for temperature, pressure
97 : real(wp) :: locpress ! needed to find location in pressure grid
98 : real(wp) :: ratio_eta_half ! ratio of vmrs of major species that defines eta=0.5
99 : ! for given flavor and reference temperature level
100 : real(wp) :: eta, feta ! binary_species_parameter, interpolation variable for eta
101 : real(wp) :: loceta ! needed to find location in eta grid
102 : real(wp) :: ftemp_term
103 : ! -----------------
104 : ! local indexes
105 : integer :: icol, ilay, iflav, igases(2), itropo, itemp
106 :
107 107862905 : do ilay = 1, nlay
108 1759339505 : do icol = 1, ncol
109 : ! index and factor for temperature interpolation
110 1651476600 : jtemp(icol,ilay) = int((tlay(icol,ilay) - (temp_ref_min - temp_ref_delta)) / temp_ref_delta)
111 1651476600 : jtemp(icol,ilay) = min(ntemp - 1, max(1, jtemp(icol,ilay))) ! limit the index range
112 1651476600 : ftemp(icol,ilay) = (tlay(icol,ilay) - temp_ref(jtemp(icol,ilay))) / temp_ref_delta
113 :
114 : ! index and factor for pressure interpolation
115 1651476600 : locpress = 1._wp + (log(play(icol,ilay)) - press_ref_log(1)) / press_ref_log_delta
116 1651476600 : jpress(icol,ilay) = min(npres-1, max(1, int(locpress)))
117 1651476600 : fpress(icol,ilay) = locpress - float(jpress(icol,ilay))
118 :
119 : ! determine if in lower or upper part of atmosphere
120 1758204106 : tropo(icol,ilay) = log(play(icol,ilay)) > press_ref_trop_log
121 : end do
122 : end do
123 :
124 12100126 : do iflav = 1, nflav
125 32894181 : igases(:) = flavor(:,iflav)
126 1042784464 : do ilay = 1, nlay
127 17005922865 : do icol = 1, ncol
128 : ! itropo = 1 lower atmosphere; itropo = 2 upper atmosphere
129 15964273800 : itropo = merge(1,2,tropo(icol,ilay))
130 : ! loop over implemented combinations of major species
131 48923505738 : do itemp = 1, 2
132 : ! compute interpolation fractions needed for lower, then upper reference temperature level
133 : ! compute binary species parameter (eta) for flavor and temperature and
134 : ! associated interpolation index and factors
135 95785642800 : ratio_eta_half = vmr_ref(itropo,igases(1),(jtemp(icol,ilay)+itemp-1)) / &
136 >12771*10^7 : vmr_ref(itropo,igases(2),(jtemp(icol,ilay)+itemp-1))
137 31928547600 : col_mix(itemp,icol,ilay,iflav) = col_gas(icol,ilay,igases(1)) + ratio_eta_half * col_gas(icol,ilay,igases(2))
138 : ! Keep this commented lines. Fortran does allow for
139 : ! substantial optimizations and in this merge cases may
140 : ! happen that all expressions are evaluated and so create
141 : ! a division by zero. In the if construct this should be
142 : ! save. Merge is the way to do it in general inside of
143 : ! loops, but sometimes it may not work.
144 : !
145 : ! eta = merge(col_gas(icol,ilay,igases(1)) / col_mix(itemp,icol,ilay,iflav), 0.5_wp, &
146 : ! col_mix(itemp,icol,ilay,iflav) > 2._wp * tiny(col_mix))
147 : !
148 : ! In essence: do not turn it back to merge(...)!
149 31928547600 : if (col_mix(itemp,icol,ilay,iflav) > 2._wp * tiny(col_mix)) then
150 31928547600 : eta = col_gas(icol,ilay,igases(1)) / col_mix(itemp,icol,ilay,iflav)
151 : else
152 : eta = 0.5_wp
153 : endif
154 31928547600 : loceta = eta * float(neta-1)
155 31928547600 : jeta(itemp,icol,ilay,iflav) = min(int(loceta)+1, neta-1)
156 31928547600 : feta = mod(loceta, 1.0_wp)
157 : ! compute interpolation fractions needed for minor species
158 : ! ftemp_term = (1._wp-ftemp(icol,ilay)) for itemp = 1, ftemp(icol,ilay) for itemp=2
159 31928547600 : ftemp_term = (real(2-itemp, wp) + real(2*itemp-3, wp) * ftemp(icol,ilay))
160 31928547600 : fminor(1,itemp,icol,ilay,iflav) = (1._wp-feta) * ftemp_term
161 31928547600 : fminor(2,itemp,icol,ilay,iflav) = feta * ftemp_term
162 : ! compute interpolation fractions needed for major species
163 31928547600 : fmajor(1,1,itemp,icol,ilay,iflav) = (1._wp-fpress(icol,ilay)) * fminor(1,itemp,icol,ilay,iflav)
164 31928547600 : fmajor(2,1,itemp,icol,ilay,iflav) = (1._wp-fpress(icol,ilay)) * fminor(2,itemp,icol,ilay,iflav)
165 31928547600 : fmajor(1,2,itemp,icol,ilay,iflav) = fpress(icol,ilay) * fminor(1,itemp,icol,ilay,iflav)
166 47892821400 : fmajor(2,2,itemp,icol,ilay,iflav) = fpress(icol,ilay) * fminor(2,itemp,icol,ilay,iflav)
167 : end do ! reference temperatures
168 : end do ! icol
169 : end do ! ilay
170 : end do ! iflav
171 :
172 1135399 : end subroutine interpolation
173 : ! --------------------------------------------------------------------------------------
174 : !
175 : !> Compute minor and major species optical depth using pre-computed interpolation coefficients
176 : !> (jeta,jtemp,jpress) and weights (fmajor, fminor)
177 : !
178 1135399 : subroutine compute_tau_absorption( &
179 : ncol,nlay,nbnd,ngpt, & ! dimensions
180 : ngas,nflav,neta,npres,ntemp, &
181 : nminorlower, nminorklower, & ! number of minor contributors, total num absorption coeffs
182 : nminorupper, nminorkupper, &
183 : idx_h2o, &
184 1135399 : gpoint_flavor, &
185 1135399 : band_lims_gpt, &
186 1135399 : kmajor, &
187 1135399 : kminor_lower, &
188 1135399 : kminor_upper, &
189 1135399 : minor_limits_gpt_lower, &
190 1135399 : minor_limits_gpt_upper, &
191 1135399 : minor_scales_with_density_lower, &
192 1135399 : minor_scales_with_density_upper, &
193 1135399 : scale_by_complement_lower, &
194 1135399 : scale_by_complement_upper, &
195 1135399 : idx_minor_lower, &
196 1135399 : idx_minor_upper, &
197 1135399 : idx_minor_scaling_lower, &
198 1135399 : idx_minor_scaling_upper, &
199 1135399 : kminor_start_lower, &
200 1135399 : kminor_start_upper, &
201 1135399 : tropo, &
202 1135399 : col_mix,fmajor,fminor, &
203 1135399 : play,tlay,col_gas, &
204 1135399 : jeta,jtemp,jpress, &
205 1135399 : tau) bind(C, name="rrtmgp_compute_tau_absorption")
206 : ! ---------------------
207 : ! input dimensions
208 : integer, intent(in) :: ncol,nlay,nbnd,ngpt !! array sizes
209 : integer, intent(in) :: ngas,nflav,neta,npres,ntemp !! tables sizes
210 : integer, intent(in) :: nminorlower, nminorklower,nminorupper, nminorkupper
211 : !! table sizes
212 : integer, intent(in) :: idx_h2o !! index of water vapor in col_gas
213 : ! ---------------------
214 : ! inputs from object
215 : integer, dimension(2,ngpt), intent(in) :: gpoint_flavor
216 : !! major gas flavor (pair) by upper/lower, g-point
217 : integer, dimension(2,nbnd), intent(in) :: band_lims_gpt
218 : !! beginning and ending g-point for each band
219 : real(wp), dimension(ntemp,neta,npres+1,ngpt), intent(in) :: kmajor
220 : !! absorption coefficient table - major gases
221 : real(wp), dimension(ntemp,neta,nminorklower), intent(in) :: kminor_lower
222 : !! absorption coefficient table - minor gases, lower atmosphere
223 : real(wp), dimension(ntemp,neta,nminorkupper), intent(in) :: kminor_upper
224 : !! absorption coefficient table - minor gases, upper atmosphere
225 : integer, dimension(2,nminorlower), intent(in) :: minor_limits_gpt_lower
226 : !! beginning and ending g-point for each minor gas
227 : integer, dimension(2,nminorupper), intent(in) :: minor_limits_gpt_upper
228 : logical(wl), dimension( nminorlower), intent(in) :: minor_scales_with_density_lower
229 : !! generic treatment of minor gases - scales with density (e.g. continuum, collision-induced absorption)?
230 : logical(wl), dimension( nminorupper), intent(in) :: minor_scales_with_density_upper
231 : logical(wl), dimension( nminorlower), intent(in) :: scale_by_complement_lower
232 : !! generic treatment of minor gases - scale by density (e.g. self-continuum) or complement?
233 : logical(wl), dimension( nminorupper), intent(in) :: scale_by_complement_upper
234 : integer, dimension( nminorlower), intent(in) :: idx_minor_lower
235 : !! index of each minor gas in col_gas
236 : integer, dimension( nminorupper), intent(in) :: idx_minor_upper
237 : integer, dimension( nminorlower), intent(in) :: idx_minor_scaling_lower
238 : !! for this minor gas, index of the "scaling gas" in col_gas
239 : integer, dimension( nminorupper), intent(in) :: idx_minor_scaling_upper
240 : integer, dimension( nminorlower), intent(in) :: kminor_start_lower
241 : !! starting g-point index in minor gas absorption table
242 : integer, dimension( nminorupper), intent(in) :: kminor_start_upper
243 : logical(wl), dimension(ncol,nlay), intent(in) :: tropo
244 : !! use upper- or lower-atmospheric tables?
245 : ! ---------------------
246 : ! inputs from profile or parent function
247 : real(wp), dimension(2, ncol,nlay,nflav ), intent(in) :: col_mix
248 : !! combination of major species's column amounts - computed in interpolation()
249 : real(wp), dimension(2,2,2,ncol,nlay,nflav ), intent(in) :: fmajor
250 : !! interpolation weights for major gases - computed in interpolation()
251 : real(wp), dimension(2,2, ncol,nlay,nflav ), intent(in) :: fminor
252 : !! interpolation weights for minor gases - computed in interpolation()
253 : real(wp), dimension( ncol,nlay ), intent(in) :: play, tlay
254 : !! input temperature and pressure
255 : real(wp), dimension( ncol,nlay,0:ngas), intent(in) :: col_gas
256 : !! input column gas amount (molecules/cm^2)
257 : integer, dimension(2, ncol,nlay,nflav ), intent(in) :: jeta
258 : !! interpolation indexes in eta - computed in interpolation()
259 : integer, dimension( ncol,nlay ), intent(in) :: jtemp
260 : !! interpolation indexes in temperature - computed in interpolation()
261 : integer, dimension( ncol,nlay ), intent(in) :: jpress
262 : !! interpolation indexes in pressure - computed in interpolation()
263 : ! ---------------------
264 : ! output - optical depth
265 : real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau !! aborption optional depth
266 : ! ---------------------
267 : ! Local variables
268 : !
269 : logical :: top_at_1
270 2270798 : integer, dimension(ncol,2) :: itropo_lower, itropo_upper
271 : ! ----------------------------------------------------------------
272 :
273 : ! ---------------------
274 : ! Layer limits of upper, lower atmospheres
275 : ! ---------------------
276 1135399 : top_at_1 = play(1,1) < play(1, nlay)
277 1135399 : if(top_at_1) then
278 1135399 : itropo_lower(:, 1) = minloc(play, dim=2, mask=tropo)
279 18704299 : itropo_lower(:, 2) = nlay
280 18704299 : itropo_upper(:, 1) = 1
281 1759339505 : itropo_upper(:, 2) = maxloc(play, dim=2, mask=(.not. tropo))
282 : else
283 0 : itropo_lower(:, 1) = 1
284 0 : itropo_lower(:, 2) = minloc(play, dim=2, mask= tropo)
285 0 : itropo_upper(:, 1) = maxloc(play, dim=2, mask=(.not. tropo))
286 0 : itropo_upper(:, 2) = nlay
287 : end if
288 : ! ---------------------
289 : ! Major Species
290 : ! ---------------------
291 : call gas_optical_depths_major( &
292 : ncol,nlay,nbnd,ngpt, & ! dimensions
293 : nflav,neta,npres,ntemp, &
294 : gpoint_flavor, &
295 : band_lims_gpt, &
296 : kmajor, &
297 : col_mix,fmajor, &
298 : jeta,tropo,jtemp,jpress, &
299 1135399 : tau)
300 : ! ---------------------
301 : ! Minor Species - lower
302 : ! ---------------------
303 : call gas_optical_depths_minor( &
304 : ncol,nlay,ngpt, & ! dimensions
305 : ngas,nflav,ntemp,neta, &
306 : nminorlower,nminorklower, &
307 : idx_h2o, &
308 : gpoint_flavor(1,:), &
309 : kminor_lower, &
310 : minor_limits_gpt_lower, &
311 : minor_scales_with_density_lower, &
312 : scale_by_complement_lower, &
313 : idx_minor_lower, &
314 : idx_minor_scaling_lower, &
315 : kminor_start_lower, &
316 : play, tlay, &
317 : col_gas,fminor,jeta, &
318 : itropo_lower,jtemp, &
319 140238263 : tau)
320 : ! ---------------------
321 : ! Minor Species - upper
322 : ! ---------------------
323 : call gas_optical_depths_minor( &
324 : ncol,nlay,ngpt, & ! dimensions
325 : ngas,nflav,ntemp,neta, &
326 : nminorupper,nminorkupper, &
327 : idx_h2o, &
328 : gpoint_flavor(2,:), &
329 : kminor_upper, &
330 : minor_limits_gpt_upper, &
331 : minor_scales_with_density_upper, &
332 : scale_by_complement_upper, &
333 : idx_minor_upper, &
334 : idx_minor_scaling_upper, &
335 : kminor_start_upper, &
336 : play, tlay, &
337 : col_gas,fminor,jeta, &
338 : itropo_upper,jtemp, &
339 140238263 : tau)
340 1135399 : end subroutine compute_tau_absorption
341 : ! --------------------------------------------------------------------------------------
342 :
343 : ! --------------------------------------------------------------------------------------
344 : !
345 : ! compute minor species optical depths
346 : !
347 1135399 : subroutine gas_optical_depths_major(ncol,nlay,nbnd,ngpt,&
348 : nflav,neta,npres,ntemp, & ! dimensions
349 1135399 : gpoint_flavor, band_lims_gpt, & ! inputs from object
350 1135399 : kmajor, &
351 1135399 : col_mix,fmajor, &
352 1135399 : jeta,tropo,jtemp,jpress, & ! local input
353 1135399 : tau)
354 : ! input dimensions
355 : integer, intent(in) :: ncol, nlay, nbnd, ngpt, nflav,neta,npres,ntemp ! dimensions
356 :
357 : ! inputs from object
358 : integer, dimension(2,ngpt), intent(in) :: gpoint_flavor
359 : integer, dimension(2,nbnd), intent(in) :: band_lims_gpt ! start and end g-point for each band
360 : real(wp), dimension(ntemp,neta,npres+1,ngpt), intent(in) :: kmajor
361 :
362 : ! inputs from profile or parent function
363 : real(wp), dimension(2, ncol,nlay,nflav), intent(in) :: col_mix
364 : real(wp), dimension(2,2,2,ncol,nlay,nflav), intent(in) :: fmajor
365 : integer, dimension(2, ncol,nlay,nflav), intent(in) :: jeta
366 : logical(wl), dimension(ncol,nlay), intent(in) :: tropo
367 : integer, dimension(ncol,nlay), intent(in) :: jtemp, jpress
368 :
369 : ! outputs
370 : real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau
371 : ! -----------------
372 : ! local variables
373 2270798 : real(wp) :: tau_major(ngpt) ! major species optical depth
374 : ! local index
375 : integer :: icol, ilay, iflav, ibnd, itropo
376 : integer :: gptS, gptE
377 :
378 : ! optical depth calculation for major species
379 18523257 : do ibnd = 1, nbnd
380 17387858 : gptS = band_lims_gpt(1, ibnd)
381 17387858 : gptE = band_lims_gpt(2, ibnd)
382 1652981909 : do ilay = 1, nlay
383 26974487710 : do icol = 1, ncol
384 : ! itropo = 1 lower atmosphere; itropo = 2 upper atmosphere
385 25322641200 : itropo = merge(1,2,tropo(icol,ilay))
386 25322641200 : iflav = gpoint_flavor(itropo, gptS) !eta interpolation depends on band's flavor
387 0 : tau_major(gptS:gptE) = &
388 : ! interpolation in temperature, pressure, and eta
389 : interpolate3D_byflav(col_mix(:,icol,ilay,iflav), &
390 : fmajor(:,:,:,icol,ilay,iflav), kmajor, &
391 : band_lims_gpt(1, ibnd), band_lims_gpt(2, ibnd), &
392 25322641200 : jeta(:,icol,ilay,iflav), jtemp(icol,ilay),jpress(icol,ilay)+itropo)
393 >22953*10^7 : tau(icol,ilay,gptS:gptE) = tau(icol,ilay,gptS:gptE) + tau_major(gptS:gptE)
394 : end do
395 : end do
396 : end do
397 :
398 1135399 : end subroutine gas_optical_depths_major
399 :
400 : ! ----------------------------------------------------------
401 : !
402 : ! compute minor species optical depths
403 : !
404 2270798 : subroutine gas_optical_depths_minor(ncol,nlay,ngpt, &
405 : ngas,nflav,ntemp,neta, &
406 : nminor,nminork, &
407 : idx_h2o, &
408 2270798 : gpt_flv, &
409 2270798 : kminor, &
410 2270798 : minor_limits_gpt, &
411 2270798 : minor_scales_with_density, &
412 2270798 : scale_by_complement, &
413 2270798 : idx_minor, idx_minor_scaling, &
414 2270798 : kminor_start, &
415 2270798 : play, tlay, &
416 2270798 : col_gas,fminor,jeta, &
417 2270798 : layer_limits,jtemp, &
418 2270798 : tau)
419 : integer, intent(in ) :: ncol,nlay,ngpt
420 : integer, intent(in ) :: ngas,nflav
421 : integer, intent(in ) :: ntemp,neta,nminor,nminork
422 : integer, intent(in ) :: idx_h2o
423 : integer, dimension(ngpt), intent(in ) :: gpt_flv
424 : real(wp), dimension(ntemp,neta,nminork), intent(in ) :: kminor
425 : integer, dimension(2,nminor), intent(in ) :: minor_limits_gpt
426 : logical(wl), dimension( nminor), intent(in ) :: minor_scales_with_density
427 : logical(wl), dimension( nminor), intent(in ) :: scale_by_complement
428 : integer, dimension( nminor), intent(in ) :: kminor_start
429 : integer, dimension( nminor), intent(in ) :: idx_minor, idx_minor_scaling
430 : real(wp), dimension(ncol,nlay), intent(in ) :: play, tlay
431 : real(wp), dimension(ncol,nlay,0:ngas), intent(in ) :: col_gas
432 : real(wp), dimension(2,2,ncol,nlay,nflav), intent(in ) :: fminor
433 : integer, dimension(2, ncol,nlay,nflav), intent(in ) :: jeta
434 : integer, dimension(ncol, 2), intent(in ) :: layer_limits
435 : integer, dimension(ncol,nlay), intent(in ) :: jtemp
436 : real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau
437 : ! -----------------
438 : ! local variables
439 : real(wp), parameter :: PaTohPa = 0.01_wp
440 : real(wp) :: vmr_fact, dry_fact ! conversion from column abundance to dry vol. mixing ratio;
441 : real(wp) :: scaling ! optical depth
442 : integer :: icol, ilay, iflav, imnr
443 : integer :: gptS, gptE
444 2270798 : real(wp), dimension(ngpt) :: tau_minor
445 : ! -----------------
446 : !
447 : ! Guard against layer limits being 0 -- that means don't do anything i.e. there are no
448 : ! layers with pressures in the upper or lower atmosphere respectively
449 : ! First check skips the routine entirely if all columns are out of bounds...
450 : !
451 :
452 2270798 : if(any(layer_limits(:,1) > 0)) then
453 70654441 : do imnr = 1, size(scale_by_complement,dim=1) ! loop over minor absorbers in each band
454 1130644741 : do icol = 1, ncol
455 : !
456 : ! This check skips individual columns with no pressures in range
457 : !
458 1128373943 : if(layer_limits(icol,1) > 0) then
459 50868404580 : do ilay = layer_limits(icol,1), layer_limits(icol,2)
460 : !
461 : ! Scaling of minor gas absortion coefficient begins with column amount of minor gas
462 : !
463 49808414280 : scaling = col_gas(icol,ilay,idx_minor(imnr))
464 : !
465 : ! Density scaling (e.g. for h2o continuum, collision-induced absorption)
466 : !
467 49808414280 : if (minor_scales_with_density(imnr)) then
468 : !
469 : ! NOTE: P needed in hPa to properly handle density scaling.
470 : !
471 33845628712 : scaling = scaling * (PaTohPa*play(icol,ilay)/tlay(icol,ilay))
472 33845628712 : if(idx_minor_scaling(imnr) > 0) then ! there is a second gas that affects this gas's absorption
473 31919150504 : vmr_fact = 1._wp / col_gas(icol,ilay,0)
474 31919150504 : dry_fact = 1._wp / (1._wp + col_gas(icol,ilay,idx_h2o) * vmr_fact)
475 : ! scale by density of special gas
476 31919150504 : if (scale_by_complement(imnr)) then ! scale by densities of all gases but the special one
477 15959575252 : scaling = scaling * (1._wp - col_gas(icol,ilay,idx_minor_scaling(imnr)) * vmr_fact * dry_fact)
478 : else
479 15959575252 : scaling = scaling * (col_gas(icol,ilay,idx_minor_scaling(imnr)) * vmr_fact * dry_fact)
480 : endif
481 : endif
482 : endif
483 : !
484 : ! Interpolation of absorption coefficient and calculation of optical depth
485 : !
486 : ! Which gpoint range does this minor gas affect?
487 49808414280 : gptS = minor_limits_gpt(1,imnr)
488 49808414280 : gptE = minor_limits_gpt(2,imnr)
489 49808414280 : iflav = gpt_flv(gptS)
490 0 : tau_minor(gptS:gptE) = scaling * &
491 : interpolate2D_byflav(fminor(:,:,icol,ilay,iflav), &
492 : kminor, &
493 : kminor_start(imnr), kminor_start(imnr)+(gptE-gptS), &
494 >46920*10^7 : jeta(:,icol,ilay,iflav), jtemp(icol,ilay))
495 >47026*10^7 : tau(icol,ilay,gptS:gptE) = tau(icol,ilay,gptS:gptE) + tau_minor(gptS:gptE)
496 : enddo
497 : end if
498 : enddo
499 : enddo
500 : end if
501 :
502 2270798 : end subroutine gas_optical_depths_minor
503 : ! ----------------------------------------------------------
504 : !
505 : ! compute Rayleigh scattering optical depths
506 : !
507 389263 : subroutine compute_tau_rayleigh(ncol,nlay,nbnd,ngpt, &
508 : ngas,nflav,neta,npres,ntemp, &
509 389263 : gpoint_flavor,band_lims_gpt, &
510 389263 : krayl, &
511 389263 : idx_h2o, col_dry,col_gas, &
512 389263 : fminor,jeta,tropo,jtemp, &
513 389263 : tau_rayleigh) bind(C, name="rrtmgp_compute_tau_rayleigh")
514 : integer, intent(in ) :: ncol,nlay,nbnd,ngpt
515 : !! input dimensions
516 : integer, intent(in ) :: ngas,nflav,neta,npres,ntemp
517 : !! table dimensions
518 : integer, dimension(2,ngpt), intent(in ) :: gpoint_flavor
519 : !! major gas flavor (pair) by upper/lower, g-point
520 : integer, dimension(2,nbnd), intent(in ) :: band_lims_gpt
521 : !! start and end g-point for each band
522 : real(wp), dimension(ntemp,neta,ngpt,2), intent(in ) :: krayl
523 : !! Rayleigh scattering coefficients
524 : integer, intent(in ) :: idx_h2o
525 : !! index of water vapor in col_gas
526 : real(wp), dimension(ncol,nlay), intent(in ) :: col_dry
527 : !! column amount of dry air
528 : real(wp), dimension(ncol,nlay,0:ngas), intent(in ) :: col_gas
529 : !! input column gas amount (molecules/cm^2)
530 : real(wp), dimension(2,2,ncol,nlay,nflav), intent(in ) :: fminor
531 : !! interpolation weights for major gases - computed in interpolation()
532 : integer, dimension(2, ncol,nlay,nflav), intent(in ) :: jeta
533 : !! interpolation indexes in eta - computed in interpolation()
534 : logical(wl), dimension(ncol,nlay), intent(in ) :: tropo
535 : !! use upper- or lower-atmospheric tables?
536 : integer, dimension(ncol,nlay), intent(in ) :: jtemp
537 : !! interpolation indexes in temperature - computed in interpolation()
538 : ! outputs
539 : real(wp), dimension(ncol,nlay,ngpt), intent(out) :: tau_rayleigh
540 : !! Rayleigh optical depth
541 : ! -----------------
542 : ! local variables
543 778526 : real(wp) :: k(ngpt) ! rayleigh scattering coefficient
544 : integer :: icol, ilay, iflav, ibnd, gptS, gptE
545 : integer :: itropo
546 : ! -----------------
547 :
548 5838945 : do ibnd = 1, nbnd
549 5449682 : gptS = band_lims_gpt(1, ibnd)
550 5449682 : gptE = band_lims_gpt(2, ibnd)
551 518109053 : do ilay = 1, nlay
552 8224610590 : do icol = 1, ncol
553 7706890800 : itropo = merge(1,2,tropo(icol,ilay)) ! itropo = 1 lower atmosphere;itropo = 2 upper atmosphere
554 7706890800 : iflav = gpoint_flavor(itropo, gptS) !eta interpolation depends on band's flavor
555 0 : k(gptS:gptE) = interpolate2D_byflav(fminor(:,:,icol,ilay,iflav), &
556 : krayl(:,:,:,itropo), &
557 7706890800 : gptS, gptE, jeta(:,icol,ilay,iflav), jtemp(icol,ilay))
558 0 : tau_rayleigh(icol,ilay,gptS:gptE) = k(gptS:gptE) * &
559 69874287308 : (col_gas(icol,ilay,idx_h2o)+col_dry(icol,ilay))
560 : end do
561 : end do
562 : end do
563 :
564 389263 : end subroutine compute_tau_rayleigh
565 :
566 : ! ----------------------------------------------------------
567 746136 : subroutine compute_Planck_source( &
568 : ncol, nlay, nbnd, ngpt, &
569 : nflav, neta, npres, ntemp, nPlanckTemp,&
570 746136 : tlay, tlev, tsfc, sfc_lay, &
571 746136 : fmajor, jeta, tropo, jtemp, jpress, &
572 746136 : gpoint_bands, band_lims_gpt, &
573 746136 : pfracin, temp_ref_min, totplnk_delta, totplnk, gpoint_flavor, &
574 746136 : sfc_src, lay_src, lev_src_inc, lev_src_dec, sfc_source_Jac) bind(C, name="rrtmgp_compute_Planck_source")
575 : integer, intent(in) :: ncol, nlay, nbnd, ngpt
576 : !! input dimensions
577 : integer, intent(in) :: nflav, neta, npres, ntemp, nPlanckTemp
578 : !! table dimensions
579 : real(wp), dimension(ncol,nlay ), intent(in) :: tlay !! temperature at layer centers (K)
580 : real(wp), dimension(ncol,nlay+1), intent(in) :: tlev !! temperature at interfaces (K)
581 : real(wp), dimension(ncol ), intent(in) :: tsfc !! surface temperture
582 : integer, intent(in) :: sfc_lay !! index into surface layer
583 : ! Interpolation variables
584 : real(wp), dimension(2,2,2,ncol,nlay,nflav), intent(in) :: fmajor
585 : !! interpolation weights for major gases - computed in interpolation()
586 : integer, dimension(2, ncol,nlay,nflav), intent(in) :: jeta
587 : !! interpolation indexes in eta - computed in interpolation()
588 : logical(wl), dimension( ncol,nlay), intent(in) :: tropo
589 : !! use upper- or lower-atmospheric tables?
590 : integer, dimension( ncol,nlay), intent(in) :: jtemp, jpress
591 : !! interpolation indexes in temperature and pressure - computed in interpolation()
592 : ! Table-specific
593 : integer, dimension(ngpt), intent(in) :: gpoint_bands !! band to which each g-point belongs
594 : integer, dimension(2, nbnd), intent(in) :: band_lims_gpt !! start and end g-point for each band
595 : real(wp), intent(in) :: temp_ref_min, totplnk_delta !! interpolation constants
596 : real(wp), dimension(ntemp,neta,npres+1,ngpt), intent(in) :: pfracin !! Fraction of the Planck function in each g-point
597 : real(wp), dimension(nPlanckTemp,nbnd), intent(in) :: totplnk !! Total Planck function by band at each temperature
598 : integer, dimension(2,ngpt), intent(in) :: gpoint_flavor !! major gas flavor (pair) by upper/lower, g-point
599 :
600 : real(wp), dimension(ncol, ngpt), intent(out) :: sfc_src !! Planck emssion from the surface
601 : real(wp), dimension(ncol,nlay,ngpt), intent(out) :: lay_src !! Planck emssion from layer centers
602 : real(wp), dimension(ncol,nlay,ngpt), intent(out) :: lev_src_inc, lev_src_dec
603 : !! Planck emission at layer boundaries, using spectral mapping in the direction of propagation
604 : real(wp), dimension(ncol, ngpt), intent(out) :: sfc_source_Jac
605 : !! Jacobian (derivative) of the surface Planck source with respect to surface temperature
606 : ! -----------------
607 : ! local
608 : real(wp), parameter :: delta_Tsurf = 1.0_wp
609 :
610 : integer :: ilay, icol, igpt, ibnd, itropo, iflav
611 : integer :: gptS, gptE
612 : real(wp), dimension(2), parameter :: one = [1._wp, 1._wp]
613 1492272 : real(wp) :: pfrac (ncol,nlay ,ngpt)
614 1492272 : real(wp) :: planck_function(ncol,nlay+1,nbnd)
615 : ! -----------------
616 :
617 : ! Calculation of fraction of band's Planck irradiance associated with each g-point
618 12684312 : do ibnd = 1, nbnd
619 11938176 : gptS = band_lims_gpt(1, ibnd)
620 11938176 : gptE = band_lims_gpt(2, ibnd)
621 1134872856 : do ilay = 1, nlay
622 18749877120 : do icol = 1, ncol
623 : ! itropo = 1 lower atmosphere; itropo = 2 upper atmosphere
624 17615750400 : itropo = merge(1,2,tropo(icol,ilay))
625 17615750400 : iflav = gpoint_flavor(itropo, gptS) !eta interpolation depends on band's flavor
626 0 : pfrac(icol,ilay,gptS:gptE) = &
627 : ! interpolation in temperature, pressure, and eta
628 : interpolate3D_byflav(one, fmajor(:,:,:,icol,ilay,iflav), pfracin, &
629 : band_lims_gpt(1, ibnd), band_lims_gpt(2, ibnd), &
630 >15966*10^7 : jeta(:,icol,ilay,iflav), jtemp(icol,ilay),jpress(icol,ilay)+itropo)
631 : end do ! column
632 : end do ! layer
633 : end do ! band
634 :
635 : !
636 : ! Planck function by band for the surface
637 : ! Compute surface source irradiance for g-point, equals band irradiance x fraction for g-point
638 : !
639 12458736 : do icol = 1, ncol
640 199114200 : planck_function(icol,1,1:nbnd) = interpolate1D(tsfc(icol), temp_ref_min, totplnk_delta, totplnk)
641 199114200 : planck_function(icol,2,1:nbnd) = interpolate1D(tsfc(icol) + delta_Tsurf, temp_ref_min, totplnk_delta, totplnk)
642 : !
643 : ! Map to g-points
644 : !
645 199860336 : do ibnd = 1, nbnd
646 187401600 : gptS = band_lims_gpt(1, ibnd)
647 187401600 : gptE = band_lims_gpt(2, ibnd)
648 1698327000 : do igpt = gptS, gptE
649 1499212800 : sfc_src(icol,igpt) = pfrac(icol,sfc_lay,igpt) * planck_function(icol,1,ibnd)
650 : sfc_source_Jac(icol, igpt) = pfrac(icol,sfc_lay,igpt) * &
651 1686614400 : (planck_function(icol, 2, ibnd) - planck_function(icol,1,ibnd))
652 : end do
653 : end do
654 : end do !icol
655 :
656 70882920 : do ilay = 1, nlay
657 1171867320 : do icol = 1, ncol
658 : ! Compute layer source irradiance for g-point, equals band irradiance x fraction for g-point
659 18786871584 : planck_function(icol,ilay,1:nbnd) = interpolate1D(tlay(icol,ilay), temp_ref_min, totplnk_delta, totplnk)
660 : end do
661 : end do
662 :
663 : !
664 : ! Map to g-points
665 : !
666 12684312 : do ibnd = 1, nbnd
667 11938176 : gptS = band_lims_gpt(1, ibnd)
668 11938176 : gptE = band_lims_gpt(2, ibnd)
669 108189720 : do igpt = gptS, gptE
670 9084951936 : do ilay = 1, nlay
671 >14999*10^7 : do icol = 1, ncol
672 >14990*10^7 : lay_src(icol,ilay,igpt) = pfrac(icol,ilay,igpt) * planck_function(icol,ilay,ibnd)
673 : end do
674 : end do
675 : end do
676 : end do
677 :
678 : ! compute level source irradiances for each g-point, one each for upward and downward paths
679 70882920 : do ilay = 1, nlay
680 1171867320 : do icol = 1, ncol
681 18716734800 : planck_function(icol, 1,1:nbnd) = interpolate1D(tlev(icol, 1),temp_ref_min, totplnk_delta, totplnk)
682 18786871584 : planck_function(icol,ilay+1,1:nbnd) = interpolate1D(tlev(icol,ilay+1),temp_ref_min, totplnk_delta, totplnk)
683 : end do
684 : end do
685 :
686 : !
687 : ! Map to g-points
688 : !
689 12684312 : do ibnd = 1, nbnd
690 11938176 : gptS = band_lims_gpt(1, ibnd)
691 11938176 : gptE = band_lims_gpt(2, ibnd)
692 108189720 : do igpt = gptS, gptE
693 9084951936 : do ilay = 1, nlay
694 >14999*10^7 : do icol = 1, ncol
695 >14092*10^7 : lev_src_inc(icol,ilay,igpt) = pfrac(icol,ilay,igpt) *planck_function(icol,ilay+1,ibnd)
696 >14990*10^7 : lev_src_dec(icol,ilay,igpt) = pfrac(icol,ilay,igpt) *planck_function(icol,ilay ,ibnd)
697 : end do
698 : end do
699 : end do
700 : end do
701 :
702 746136 : end subroutine compute_Planck_source
703 : ! ----------------------------------------------------------
704 : !
705 : ! One dimensional interpolation -- return all values along second table dimension
706 : !
707 3326378400 : pure function interpolate1D(val, offset, delta, table) result(res)
708 : ! input
709 : real(wp), intent(in) :: val, & ! axis value at which to evaluate table
710 : offset, & ! minimum of table axis
711 : delta ! step size of table axis
712 : real(wp), dimension(:,:), &
713 : intent(in) :: table ! dimensions (axis, values)
714 : ! output
715 : real(wp), dimension(size(table,dim=2)) :: res
716 :
717 : ! local
718 : real(wp) :: val0 ! fraction index adjusted by offset and delta
719 : integer :: index ! index term
720 : real(wp) :: frac ! fractional term
721 : ! -------------------------------------
722 3326378400 : val0 = (val - offset) / delta
723 3326378400 : frac = val0 - int(val0) ! get fractional part
724 3326378400 : index = min(size(table,dim=1)-1, max(1, int(val0)+1)) ! limit the index range
725 56548432800 : res(:) = table(index,:) + frac * (table(index+1,:) - table(index,:))
726 3326378400 : end function interpolate1D
727 : ! ----------------------------------------------------------
728 : ! This function returns a range of values from a subset (in gpoint) of the k table
729 : !
730 57515305080 : pure function interpolate2D_byflav(fminor, k, gptS, gptE, jeta, jtemp) result(res)
731 : real(wp), dimension(2,2), intent(in) :: fminor ! interpolation fractions for minor species
732 : ! index(1) : reference eta level (temperature dependent)
733 : ! index(2) : reference temperature level
734 : real(wp), dimension(:,:,:), intent(in) :: k ! (g-point, eta, temp)
735 : integer, intent(in) :: gptS, gptE, jtemp ! interpolation index for temperature
736 : integer, dimension(2), intent(in) :: jeta ! interpolation index for binary species parameter (eta)
737 : real(wp), dimension(gptE-gptS+1) :: res ! the result
738 :
739 : ! Local variable
740 : integer :: igpt
741 : ! each code block is for a different reference temperature
742 :
743 >53856*10^7 : do igpt = 1, gptE-gptS+1
744 >19241*10^8 : res(igpt) = fminor(1,1) * k(jtemp , jeta(1) , gptS+igpt-1) + &
745 >48104*10^7 : fminor(2,1) * k(jtemp , jeta(1)+1, gptS+igpt-1) + &
746 >96209*10^7 : fminor(1,2) * k(jtemp+1, jeta(2) , gptS+igpt-1) + &
747 >39059*10^8 : fminor(2,2) * k(jtemp+1, jeta(2)+1, gptS+igpt-1)
748 : end do
749 :
750 57515305080 : end function interpolate2D_byflav
751 : ! ----------------------------------------------------------
752 42938391600 : pure function interpolate3D_byflav(scaling, fmajor, k, gptS, gptE, jeta, jtemp, jpress) result(res)
753 : real(wp), dimension(2), intent(in) :: scaling
754 : real(wp), dimension(2,2,2), intent(in) :: fmajor ! interpolation fractions for major species
755 : ! index(1) : reference eta level (temperature dependent)
756 : ! index(2) : reference pressure level
757 : ! index(3) : reference temperature level
758 : real(wp), dimension(:,:,:,:),intent(in) :: k ! (temp,eta,press,gpt)
759 : integer, intent(in) :: gptS, gptE
760 : integer, dimension(2), intent(in) :: jeta ! interpolation index for binary species parameter (eta)
761 : integer, intent(in) :: jtemp ! interpolation index for temperature
762 : integer, intent(in) :: jpress ! interpolation index for pressure
763 : real(wp), dimension(gptS:gptE) :: res ! the result
764 :
765 : ! Local variable
766 : integer :: igpt
767 : ! each code block is for a different reference temperature
768 >38644*10^7 : do igpt = gptS, gptE
769 >34350*10^7 : res(igpt) = &
770 : scaling(1) * &
771 >13740*10^8 : ( fmajor(1,1,1) * k(jtemp, jeta(1) , jpress-1, igpt) + &
772 >34350*10^7 : fmajor(2,1,1) * k(jtemp, jeta(1)+1, jpress-1, igpt) + &
773 >34350*10^7 : fmajor(1,2,1) * k(jtemp, jeta(1) , jpress , igpt) + &
774 : fmajor(2,2,1) * k(jtemp, jeta(1)+1, jpress , igpt) ) + &
775 : scaling(2) * &
776 >68701*10^7 : ( fmajor(1,1,2) * k(jtemp+1, jeta(2) , jpress-1, igpt) + &
777 >34350*10^7 : fmajor(2,1,2) * k(jtemp+1, jeta(2)+1, jpress-1, igpt) + &
778 : fmajor(1,2,2) * k(jtemp+1, jeta(2) , jpress , igpt) + &
779 >34780*10^8 : fmajor(2,2,2) * k(jtemp+1, jeta(2)+1, jpress , igpt) )
780 : end do
781 42938391600 : end function interpolate3D_byflav
782 :
783 : end module mo_gas_optics_rrtmgp_kernels
|
