Line data Source code
1 : ! This code is part of RRTM for GCM Applications - Parallel (RRTMGP)
2 : !
3 : ! Contacts: Robert Pincus and Eli Mlawer
4 : ! email: rrtmgp@aer.com
5 : !
6 : ! Copyright 2015-, Atmospheric and Environmental Research,
7 : ! Regents of the University of Colorado, Trustees of Columbia University. All right reserved.
8 : !
9 : ! Use and duplication is permitted under the terms of the
10 : ! BSD 3-clause license, see http://opensource.org/licenses/BSD-3-Clause
11 : ! -------------------------------------------------------------------------------------------------
12 : !
13 : !> ## Class implementing the RRTMGP correlated-_k_ distribution
14 : !>
15 : !> Implements a class for computing spectrally-resolved gas optical properties and source functions
16 : !> given atmopsheric physical properties (profiles of temperature, pressure, and gas concentrations)
17 : !> The class must be initialized with data (provided as a netCDF file) before being used.
18 : !>
19 : !> Two variants apply to internal Planck sources (longwave radiation in the Earth's atmosphere) and to
20 : !> external stellar radiation (shortwave radiation in the Earth's atmosphere).
21 : !> The variant is chosen based on what information is supplied during initialization.
22 : ! (It might make more sense to define two sub-classes)
23 : !
24 : ! -------------------------------------------------------------------------------------------------
25 : module mo_gas_optics_rrtmgp
26 : use mo_rte_kind, only: wp, wl
27 : use mo_rte_config, only: check_extents, check_values
28 : use mo_rte_util_array, only: zero_array
29 : use mo_rte_util_array_validation, &
30 : only: any_vals_less_than, any_vals_outside, extents_are
31 : use mo_optical_props, only: ty_optical_props
32 : use mo_source_functions, only: ty_source_func_lw
33 : use mo_gas_optics_rrtmgp_kernels, &
34 : only: interpolation, compute_tau_absorption, compute_tau_rayleigh, compute_Planck_source
35 : use mo_gas_optics_constants, only: avogad, m_dry, m_h2o, grav
36 : use mo_gas_optics_util_string, only: lower_case, string_in_array, string_loc_in_array
37 : use mo_gas_concentrations, only: ty_gas_concs
38 : use mo_optical_props, only: ty_optical_props_arry, ty_optical_props_1scl, ty_optical_props_2str, ty_optical_props_nstr
39 : use mo_gas_optics, only: ty_gas_optics
40 : implicit none
41 : private
42 : real(wp), parameter :: pi = acos(-1._wp)
43 :
44 : ! -------------------------------------------------------------------------------------------------
45 : type, extends(ty_gas_optics), public :: ty_gas_optics_rrtmgp
46 : private
47 : !
48 : ! RRTMGP computes absorption in each band arising from
49 : ! two major species in each band, which are combined to make
50 : ! a relative mixing ratio eta and a total column amount (col_mix)
51 : ! contributions from zero or more minor species whose concentrations
52 : ! may be scaled by other components of the atmosphere
53 : !
54 : ! Absorption coefficients are interpolated from tables on a pressure/temperature/(eta) grid
55 : !
56 : ! ------------------------------------
57 : ! Interpolation variables: Temperature and pressure grids
58 : !
59 : real(wp), dimension(:), allocatable :: press_ref, press_ref_log, temp_ref
60 : !
61 : ! Derived and stored for convenience:
62 : ! Min and max for temperature and pressure intepolation grids
63 : ! difference in ln pressure between consecutive reference levels
64 : ! log of reference pressure separating the lower and upper atmosphere
65 : !
66 : real(wp) :: press_ref_min, press_ref_max, &
67 : temp_ref_min, temp_ref_max
68 : real(wp) :: press_ref_log_delta, temp_ref_delta, press_ref_trop_log
69 : ! ------------------------------------
70 : ! Major absorbers ("key species")
71 : ! Each unique set of major species is called a flavor.
72 : !
73 : ! Names and reference volume mixing ratios of major gases
74 : !
75 : character(32), dimension(:), allocatable :: gas_names ! gas names
76 : real(wp), dimension(:,:,:), allocatable :: vmr_ref ! vmr_ref(lower or upper atmosphere, gas, temp)
77 : !
78 : ! Which two gases are in each flavor? By index
79 : !
80 : integer, dimension(:,:), allocatable :: flavor ! major species pair; (2,nflav)
81 : !
82 : ! Which flavor for each g-point? One each for lower, upper atmosphere
83 : !
84 : integer, dimension(:,:), allocatable :: gpoint_flavor ! flavor = gpoint_flavor(2, g-point)
85 : !
86 : ! Major gas absorption coefficients
87 : !
88 : real(wp), dimension(:,:,:,:), allocatable :: kmajor ! kmajor(g-point,eta,pressure,temperature)
89 : !
90 : ! ------------------------------------
91 : ! Minor species, independently for upper and lower atmospheres
92 : ! Array extents in the n_minor dimension will differ between upper and lower atmospheres
93 : ! Each contribution has starting and ending g-points
94 : !
95 : integer, dimension(:,:), allocatable :: minor_limits_gpt_lower, &
96 : minor_limits_gpt_upper
97 : !
98 : ! Minor gas contributions might be scaled by other gas amounts; if so we need to know
99 : ! the total density and whether the contribution is scaled by the partner gas
100 : ! or its complement (i.e. all other gases)
101 : ! Water vapor self- and foreign continua work like this, as do
102 : ! all collision-induced abosption pairs
103 : !
104 : logical(wl), dimension(:), allocatable :: minor_scales_with_density_lower, &
105 : minor_scales_with_density_upper
106 : logical(wl), dimension(:), allocatable :: scale_by_complement_lower, scale_by_complement_upper
107 : integer, dimension(:), allocatable :: idx_minor_lower, idx_minor_upper
108 : integer, dimension(:), allocatable :: idx_minor_scaling_lower, idx_minor_scaling_upper
109 : !
110 : ! Index into table of absorption coefficients
111 : !
112 : integer, dimension(:), allocatable :: kminor_start_lower, kminor_start_upper
113 : !
114 : ! The absorption coefficients themselves
115 : !
116 : real(wp), dimension(:,:,:), allocatable :: kminor_lower, kminor_upper ! kminor_lower(n_minor,eta,temperature)
117 : !
118 : ! -----------------------------------------------------------------------------------
119 : !
120 : ! Rayleigh scattering coefficients
121 : !
122 : real(wp), dimension(:,:,:,:), allocatable :: krayl ! krayl(g-point,eta,temperature,upper/lower atmosphere)
123 : !
124 : ! -----------------------------------------------------------------------------------
125 : ! Planck function spectral mapping
126 : ! Allocated only when gas optics object is internal-source
127 : !
128 : real(wp), dimension(:,:,:,:), allocatable :: planck_frac ! stored fraction of Planck irradiance in band for given g-point
129 : ! planck_frac(g-point, eta, pressure, temperature)
130 : real(wp), dimension(:,:), allocatable :: totplnk ! integrated Planck irradiance by band; (Planck temperatures,band)
131 : real(wp) :: totplnk_delta ! temperature steps in totplnk
132 : real(wp), dimension(:,:), allocatable :: optimal_angle_fit ! coefficients of linear function
133 : ! of vertical path clear-sky transmittance that is used to
134 : ! determine the secant of single angle used for the
135 : ! no-scattering calculation,
136 : ! optimal_angle_fit(coefficient, band)
137 : ! -----------------------------------------------------------------------------------
138 : ! Solar source function spectral mapping with solar variability capability
139 : ! Allocated when gas optics object is external-source
140 : ! n-solar-terms: quiet sun, facular brightening and sunspot dimming components
141 : ! following the NRLSSI2 model of Coddington et al. 2016, doi:10.1175/BAMS-D-14-00265.1.
142 : !
143 : real(wp), dimension(:), allocatable :: solar_source ! incoming solar irradiance, computed from other three terms (g-point)
144 : real(wp), dimension(:), allocatable :: solar_source_quiet ! incoming solar irradiance, quiet sun term (g-point)
145 : real(wp), dimension(:), allocatable :: solar_source_facular ! incoming solar irradiance, facular term (g-point)
146 : real(wp), dimension(:), allocatable :: solar_source_sunspot ! incoming solar irradiance, sunspot term (g-point)
147 :
148 : !
149 : ! -----------------------------------------------------------------------------------
150 : ! Ancillary
151 : ! -----------------------------------------------------------------------------------
152 : ! Index into %gas_names -- is this a key species in any band?
153 : logical, dimension(:), allocatable :: is_key
154 : ! -----------------------------------------------------------------------------------
155 :
156 : contains
157 : ! Type-bound procedures
158 : ! Public procedures
159 : ! public interface
160 : generic, public :: load => load_int, load_ext
161 : procedure, public :: source_is_internal
162 : procedure, public :: source_is_external
163 : procedure, public :: is_loaded
164 : procedure, public :: finalize
165 : procedure, public :: get_ngas
166 : procedure, public :: get_gases
167 : procedure, public :: get_press_min
168 : procedure, public :: get_press_max
169 : procedure, public :: get_temp_min
170 : procedure, public :: get_temp_max
171 : procedure, public :: compute_optimal_angles
172 : procedure, public :: set_solar_variability
173 : procedure, public :: set_tsi
174 : ! Internal procedures
175 : procedure, private :: load_int
176 : procedure, private :: load_ext
177 : procedure, public :: gas_optics_int
178 : procedure, public :: gas_optics_ext
179 : procedure, private :: check_key_species_present
180 : ! Interpolation table dimensions
181 : procedure, private :: get_nflav
182 : procedure, private :: get_neta
183 : procedure, private :: get_npres
184 : procedure, private :: get_ntemp
185 : procedure, private :: get_nPlanckTemp
186 : end type ty_gas_optics_rrtmgp
187 : ! -------------------------------------------------------------------------------------------------
188 : !
189 : !> col_dry is the number of molecules per cm-2 of dry air
190 : !
191 : public :: get_col_dry ! Utility function, not type-bound
192 :
193 : contains
194 : ! --------------------------------------------------------------------------------------
195 : !
196 : ! Public procedures
197 : !
198 : ! --------------------------------------------------------------------------------------
199 : !
200 : !> Two functions to define array sizes needed by gas_optics()
201 : !
202 3798532 : pure function get_ngas(this)
203 : ! return the number of gases registered in the spectral configuration
204 : class(ty_gas_optics_rrtmgp), intent(in) :: this
205 : integer :: get_ngas
206 :
207 3798532 : get_ngas = size(this%gas_names)
208 3798532 : end function get_ngas
209 : !--------------------------------------------------------------------------------------------------------------------
210 : !
211 : !> return the number of distinct major gas pairs in the spectral bands (referred to as
212 : !> "flavors" - all bands have a flavor even if there is one or no major gas)
213 : !
214 6423131 : pure function get_nflav(this)
215 : class(ty_gas_optics_rrtmgp), intent(in) :: this
216 : integer :: get_nflav
217 :
218 6423131 : get_nflav = size(this%flavor,dim=2)
219 6423131 : end function get_nflav
220 : !--------------------------------------------------------------------------------------------------------------------
221 : !
222 : !> Compute gas optical depth and Planck source functions,
223 : !> given temperature, pressure, and composition
224 : !
225 746136 : function gas_optics_int(this, &
226 1492272 : play, plev, tlay, tsfc, gas_desc, &
227 : optical_props, sources, &
228 1492272 : col_dry, tlev) result(error_msg)
229 : ! inputs
230 : class(ty_gas_optics_rrtmgp), intent(in) :: this
231 : real(wp), dimension(:,:), intent(in ) :: play, & !! layer pressures [Pa, mb]; (ncol,nlay)
232 : plev, & !! level pressures [Pa, mb]; (ncol,nlay+1)
233 : tlay !! layer temperatures [K]; (ncol,nlay)
234 : real(wp), dimension(:), intent(in ) :: tsfc !! surface skin temperatures [K]; (ncol)
235 : type(ty_gas_concs), intent(in ) :: gas_desc !! Gas volume mixing ratios
236 : ! output
237 : class(ty_optical_props_arry), &
238 : intent(inout) :: optical_props !! Optical properties
239 : class(ty_source_func_lw ), &
240 : intent(inout) :: sources !! Planck sources
241 : character(len=128) :: error_msg !! Empty if succssful
242 : ! Optional inputs
243 : real(wp), dimension(:,:), intent(in ), &
244 : optional, target :: col_dry, & !! Column dry amount; dim(ncol,nlay)
245 : tlev !! level temperatures [K]; (ncol,nlay+1)
246 : ! ----------------------------------------------------------
247 : ! Local variables
248 : ! Interpolation coefficients for use in source function
249 1492272 : integer, dimension(size(play,dim=1), size(play,dim=2)) :: jtemp, jpress
250 1492272 : logical(wl), dimension(size(play,dim=1), size(play,dim=2)) :: tropo
251 1492272 : real(wp), dimension(2,2,2,size(play,dim=1),size(play,dim=2), get_nflav(this)) :: fmajor
252 746136 : integer, dimension(2, size(play,dim=1),size(play,dim=2), get_nflav(this)) :: jeta
253 :
254 : integer :: ncol, nlay, ngpt, nband
255 : ! ----------------------------------------------------------
256 746136 : ncol = size(play,dim=1)
257 746136 : nlay = size(play,dim=2)
258 746136 : ngpt = this%get_ngpt()
259 746136 : nband = this%get_nband()
260 : !
261 : ! Gas optics
262 : !
263 : !$acc enter data create(jtemp, jpress, tropo, fmajor, jeta)
264 : !$omp target enter data map(alloc:jtemp, jpress, tropo, fmajor, jeta)
265 : error_msg = compute_gas_taus(this, &
266 : ncol, nlay, ngpt, nband, &
267 : play, plev, tlay, gas_desc, &
268 : optical_props, &
269 : jtemp, jpress, jeta, tropo, fmajor, &
270 1492272 : col_dry)
271 746136 : if(error_msg /= '') return
272 : ! ----------------------------------------------------------
273 : !
274 : ! External source -- check arrays sizes and values
275 : ! input data sizes and values
276 : !
277 : !$acc enter data copyin(tsfc, tlev) ! Should be fine even if tlev is not supplied
278 : !$omp target enter data map(to:tsfc, tlev)
279 :
280 746136 : if(check_extents) then
281 746136 : if(.not. extents_are(tsfc, ncol)) &
282 0 : error_msg = "gas_optics(): array tsfc has wrong size"
283 746136 : if(present(tlev)) then
284 0 : if(.not. extents_are(tlev, ncol, nlay+1)) &
285 0 : error_msg = "gas_optics(): array tlev has wrong size"
286 : end if
287 : !
288 : ! output extents
289 : !
290 2984544 : if(any([sources%get_ncol(), sources%get_nlay(), sources%get_ngpt()] /= [ncol, nlay, ngpt])) &
291 0 : error_msg = "gas_optics%gas_optics: source function arrays inconsistently sized"
292 : end if
293 746136 : if(error_msg /= '') return
294 :
295 746136 : if(check_values) then
296 746136 : if(any_vals_outside(tsfc, this%temp_ref_min, this%temp_ref_max)) &
297 0 : error_msg = "gas_optics(): array tsfc has values outside range"
298 746136 : if(present(tlev)) then
299 0 : if(any_vals_outside(tlev, this%temp_ref_min, this%temp_ref_max)) &
300 0 : error_msg = "gas_optics(): array tlev has values outside range"
301 : end if
302 : end if
303 746136 : if(error_msg /= '') return
304 :
305 : !
306 : ! Interpolate source function
307 : !
308 746136 : if(present(tlev)) then
309 : !
310 : ! present status of optional argument should be passed to source()
311 : ! but isn't with PGI 19.10
312 : !
313 : error_msg = source(this, &
314 : ncol, nlay, nband, ngpt, &
315 : play, plev, tlay, tsfc, &
316 : jtemp, jpress, jeta, tropo, fmajor, &
317 : sources, &
318 0 : tlev)
319 : !$acc exit data delete(tlev)
320 : !$omp target exit data map(release:tlev)
321 : else
322 : error_msg = source(this, &
323 : ncol, nlay, nband, ngpt, &
324 : play, plev, tlay, tsfc, &
325 : jtemp, jpress, jeta, tropo, fmajor, &
326 746136 : sources)
327 : end if
328 : !$acc exit data delete(tsfc)
329 : !$omp target exit data map(release:tsfc)
330 : !$acc exit data delete(jtemp, jpress, tropo, fmajor, jeta)
331 : !$omp target exit data map(release:jtemp, jpress, tropo, fmajor, jeta)
332 : end function gas_optics_int
333 : !------------------------------------------------------------------------------------------
334 : !
335 : !> Compute gas optical depth given temperature, pressure, and composition
336 : !> Top-of-atmosphere stellar insolation is also reported
337 : !
338 389263 : function gas_optics_ext(this, &
339 389263 : play, plev, tlay, gas_desc, & ! mandatory inputs
340 389263 : optical_props, toa_src, & ! mandatory outputs
341 389263 : col_dry) result(error_msg) ! optional input
342 :
343 : class(ty_gas_optics_rrtmgp), intent(in) :: this
344 : real(wp), dimension(:,:), intent(in ) :: play, & !! layer pressures [Pa, mb]; (ncol,nlay)
345 : plev, & !! level pressures [Pa, mb]; (ncol,nlay+1)
346 : tlay !! layer temperatures [K]; (ncol,nlay)
347 : type(ty_gas_concs), intent(in ) :: gas_desc !! Gas volume mixing ratios
348 : ! output
349 : class(ty_optical_props_arry), &
350 : intent(inout) :: optical_props
351 : real(wp), dimension(:,:), intent( out) :: toa_src !! Incoming solar irradiance(ncol,ngpt)
352 : character(len=128) :: error_msg !! Empty if successful
353 :
354 : ! Optional inputs
355 : real(wp), dimension(:,:), intent(in ), &
356 : optional, target :: col_dry ! Column dry amount; dim(ncol,nlay)
357 : ! ----------------------------------------------------------
358 : ! Local variables
359 : ! Interpolation coefficients for use in source function
360 778526 : integer, dimension(size(play,dim=1), size(play,dim=2)) :: jtemp, jpress
361 778526 : logical(wl), dimension(size(play,dim=1), size(play,dim=2)) :: tropo
362 778526 : real(wp), dimension(2,2,2,size(play,dim=1),size(play,dim=2), get_nflav(this)) :: fmajor
363 389263 : integer, dimension(2, size(play,dim=1),size(play,dim=2), get_nflav(this)) :: jeta
364 :
365 : integer :: ncol, nlay, ngpt, nband, ngas, nflav
366 : integer :: igpt, icol
367 : ! ----------------------------------------------------------
368 389263 : ncol = size(play,dim=1)
369 389263 : nlay = size(play,dim=2)
370 389263 : ngpt = this%get_ngpt()
371 389263 : nband = this%get_nband()
372 389263 : ngas = this%get_ngas()
373 : nflav = get_nflav(this)
374 : !
375 : ! Gas optics
376 : !
377 : !$acc enter data create(jtemp, jpress, tropo, fmajor, jeta)
378 : !$omp target enter data map(alloc:jtemp, jpress, tropo, fmajor, jeta)
379 : error_msg = compute_gas_taus(this, &
380 : ncol, nlay, ngpt, nband, &
381 : play, plev, tlay, gas_desc, &
382 : optical_props, &
383 : jtemp, jpress, jeta, tropo, fmajor, &
384 778526 : col_dry)
385 : !$acc exit data delete(jtemp, jpress, tropo, fmajor, jeta)
386 : !$omp target exit data map(release:jtemp, jpress, tropo, fmajor, jeta)
387 389263 : if(error_msg /= '') return
388 :
389 : ! ----------------------------------------------------------
390 : !
391 : ! External source function is constant
392 : !
393 : !$acc enter data create(toa_src)
394 : !$omp target enter data map(alloc:toa_src)
395 389263 : if(check_extents) then
396 389263 : if(.not. extents_are(toa_src, ncol, ngpt)) &
397 0 : error_msg = "gas_optics(): array toa_src has wrong size"
398 : end if
399 389263 : if(error_msg /= '') return
400 :
401 : !$acc parallel loop collapse(2)
402 : !$omp target teams distribute parallel do simd collapse(2)
403 43986719 : do igpt = 1,ngpt
404 699892319 : do icol = 1,ncol
405 699503056 : toa_src(icol,igpt) = this%solar_source(igpt)
406 : end do
407 : end do
408 : !$acc exit data copyout(toa_src)
409 : !$omp target exit data map(from:toa_src)
410 : end function gas_optics_ext
411 : !------------------------------------------------------------------------------------------
412 : !
413 : ! Returns optical properties and interpolation coefficients
414 : !
415 1135399 : function compute_gas_taus(this, &
416 : ncol, nlay, ngpt, nband, &
417 1135399 : play, plev, tlay, gas_desc, &
418 : optical_props, &
419 1135399 : jtemp, jpress, jeta, tropo, fmajor, &
420 1135399 : col_dry) result(error_msg)
421 :
422 : class(ty_gas_optics_rrtmgp), &
423 : intent(in ) :: this
424 : integer, intent(in ) :: ncol, nlay, ngpt, nband
425 : real(wp), dimension(:,:), intent(in ) :: play, & ! layer pressures [Pa, mb]; (ncol,nlay)
426 : plev, & ! level pressures [Pa, mb]; (ncol,nlay+1)
427 : tlay ! layer temperatures [K]; (ncol,nlay)
428 : type(ty_gas_concs), intent(in ) :: gas_desc ! Gas volume mixing ratios
429 : class(ty_optical_props_arry), intent(inout) :: optical_props !inout because components are allocated
430 : ! Interpolation coefficients for use in internal source function
431 : integer, dimension( ncol, nlay), intent( out) :: jtemp, jpress
432 : integer, dimension(2, ncol, nlay,get_nflav(this)), intent( out) :: jeta
433 : logical(wl), dimension( ncol, nlay), intent( out) :: tropo
434 : real(wp), dimension(2,2,2,ncol, nlay,get_nflav(this)), intent( out) :: fmajor
435 : character(len=128) :: error_msg
436 :
437 : ! Optional inputs
438 : real(wp), dimension(:,:), intent(in ), &
439 : optional, target :: col_dry ! Column dry amount; dim(ncol,nlay)
440 : ! ----------------------------------------------------------
441 : ! Local variables
442 2270798 : real(wp), dimension(ncol,nlay,ngpt) :: tau, tau_rayleigh ! absorption, Rayleigh scattering optical depths
443 : ! Number of molecules per cm^2
444 2270798 : real(wp), dimension(ncol,nlay), target :: col_dry_arr
445 1135399 : real(wp), dimension(:,:), pointer :: col_dry_wk
446 : !
447 : ! Interpolation variables used in major gas but not elsewhere, so don't need exporting
448 : !
449 3406197 : real(wp), dimension(ncol,nlay, this%get_ngas()) :: vmr ! volume mixing ratios
450 3406197 : real(wp), dimension(ncol,nlay,0:this%get_ngas()) :: col_gas ! column amounts for each gas, plus col_dry
451 2270798 : real(wp), dimension(2, ncol,nlay,get_nflav(this)) :: col_mix ! combination of major species's column amounts
452 : ! index(1) : reference temperature level
453 : ! index(2) : flavor
454 : ! index(3) : layer
455 2270798 : real(wp), dimension(2,2, ncol,nlay,get_nflav(this)) :: fminor ! interpolation fractions for minor species
456 : ! index(1) : reference eta level (temperature dependent)
457 : ! index(2) : reference temperature level
458 : ! index(3) : flavor
459 : ! index(4) : layer
460 : integer :: ngas, nflav, neta, npres, ntemp
461 : integer :: icol, ilay, igas
462 : integer :: idx_h2o ! index of water vapor
463 : integer :: nminorlower, nminorklower,nminorupper, nminorkupper
464 : logical :: use_rayl
465 : ! ----------------------------------------------------------
466 : !
467 : ! Error checking
468 : !
469 1135399 : use_rayl = allocated(this%krayl)
470 1135399 : error_msg = ''
471 : ! Check for initialization
472 1135399 : if (.not. this%is_loaded()) then
473 0 : error_msg = 'ERROR: spectral configuration not loaded'
474 0 : return
475 : end if
476 : !
477 : ! Check for presence of key species in ty_gas_concs; return error if any key species are not present
478 : !
479 1135399 : error_msg = this%check_key_species_present(gas_desc)
480 1135399 : if (error_msg /= '') return
481 :
482 : !
483 : ! Check input data sizes and values
484 : !
485 : !$acc data copyin(play,plev,tlay) create( vmr,col_gas)
486 : !$omp target data map(to:play,plev,tlay) map(alloc:vmr,col_gas)
487 1135399 : if(check_extents) then
488 1135399 : if(.not. extents_are(play, ncol, nlay )) &
489 0 : error_msg = "gas_optics(): array play has wrong size"
490 1135399 : if(.not. extents_are(tlay, ncol, nlay )) &
491 0 : error_msg = "gas_optics(): array tlay has wrong size"
492 1135399 : if(.not. extents_are(plev, ncol, nlay+1)) &
493 0 : error_msg = "gas_optics(): array plev has wrong size"
494 : if(optical_props%get_ncol() /= ncol .or. &
495 1135399 : optical_props%get_nlay() /= nlay .or. &
496 : optical_props%get_ngpt() /= ngpt) &
497 0 : error_msg = "gas_optics(): optical properties have the wrong extents"
498 1135399 : if(present(col_dry)) then
499 0 : if(.not. extents_are(col_dry, ncol, nlay)) &
500 0 : error_msg = "gas_optics(): array col_dry has wrong size"
501 : end if
502 : end if
503 :
504 1135399 : if(error_msg == '') then
505 1135399 : if(check_values) then
506 1135399 : if(any_vals_outside(play, this%press_ref_min,this%press_ref_max)) &
507 0 : error_msg = "gas_optics(): array play has values outside range"
508 1135399 : if(any_vals_less_than(plev, 0._wp)) &
509 0 : error_msg = "gas_optics(): array plev has values outside range"
510 1135399 : if(any_vals_outside(tlay, this%temp_ref_min, this%temp_ref_max)) &
511 0 : error_msg = "gas_optics(): array tlay has values outside range"
512 1135399 : if(present(col_dry)) then
513 0 : if(any_vals_less_than(col_dry, 0._wp)) &
514 0 : error_msg = "gas_optics(): array col_dry has values outside range"
515 : end if
516 : end if
517 : end if
518 :
519 : ! ----------------------------------------------------------
520 1135399 : if(error_msg == '') then
521 1135399 : ngas = this%get_ngas()
522 1135399 : nflav = get_nflav(this)
523 1135399 : neta = this%get_neta()
524 1135399 : npres = this%get_npres()
525 1135399 : ntemp = this%get_ntemp()
526 : ! number of minor contributors, total num absorption coeffs
527 1135399 : nminorlower = size(this%minor_scales_with_density_lower)
528 1135399 : nminorklower = size(this%kminor_lower, 3)
529 1135399 : nminorupper = size(this%minor_scales_with_density_upper)
530 1135399 : nminorkupper = size(this%kminor_upper, 3)
531 : !
532 : ! Fill out the array of volume mixing ratios
533 : !
534 10218591 : do igas = 1, ngas
535 : !
536 : ! Get vmr if gas is provided in ty_gas_concs
537 : !
538 60176147 : if (any (lower_case(this%gas_names(igas)) == gas_desc%get_gas_names())) then
539 27249576 : error_msg = gas_desc%get_vmr(this%gas_names(igas), vmr(:,:,igas))
540 : endif
541 : end do
542 : end if
543 :
544 1135399 : if(error_msg == '') then
545 :
546 : !
547 : ! Painful hacks to get code to compile with both the CCE-14 and Nvidia 21.3 compiler
548 : !
549 : #ifdef _CRAYFTN
550 : !$acc enter data copyin(optical_props)
551 : #endif
552 : select type(optical_props)
553 : type is (ty_optical_props_1scl)
554 : #ifndef _CRAYFTN
555 : !$acc enter data copyin(optical_props)
556 : #endif
557 : !$acc enter data create( optical_props%tau)
558 : !$omp target enter data map(alloc:optical_props%tau)
559 : type is (ty_optical_props_2str)
560 : #ifndef _CRAYFTN
561 : !$acc enter data copyin(optical_props)
562 : #endif
563 : !$acc enter data create( optical_props%tau, optical_props%ssa, optical_props%g)
564 : !$omp target enter data map(alloc:optical_props%tau, optical_props%ssa, optical_props%g)
565 : type is (ty_optical_props_nstr)
566 : #ifndef _CRAYFTN
567 : !$acc enter data copyin(optical_props)
568 : #endif
569 : !$acc enter data create( optical_props%tau, optical_props%ssa, optical_props%p)
570 : !$omp target enter data map(alloc:optical_props%tau, optical_props%ssa, optical_props%p)
571 : end select
572 :
573 : !
574 : ! Compute dry air column amounts (number of molecule per cm^2) if user hasn't provided them
575 : !
576 1135399 : idx_h2o = string_loc_in_array('h2o', this%gas_names)
577 1135399 : if (present(col_dry)) then
578 : !$acc enter data copyin(col_dry)
579 : !$omp target enter data map(to:col_dry)
580 0 : col_dry_wk => col_dry
581 : else
582 : !$acc enter data create( col_dry_arr)
583 : !$omp target enter data map(alloc:col_dry_arr)
584 1135399 : col_dry_arr = get_col_dry(vmr(:,:,idx_h2o), plev) ! dry air column amounts computation
585 1135399 : col_dry_wk => col_dry_arr
586 : end if
587 : !
588 : ! compute column gas amounts [molec/cm^2]
589 : !
590 : !$acc parallel loop gang vector collapse(2)
591 : !$omp target teams distribute parallel do simd collapse(2)
592 107862905 : do ilay = 1, nlay
593 1759339505 : do icol = 1, ncol
594 1758204106 : col_gas(icol,ilay,0) = col_dry_wk(icol,ilay)
595 : end do
596 : end do
597 : !$acc parallel loop gang vector collapse(3)
598 : !$omp target teams distribute parallel do simd collapse(3)
599 10218591 : do igas = 1, ngas
600 864038639 : do ilay = 1, nlay
601 14074716040 : do icol = 1, ncol
602 14065632848 : col_gas(icol,ilay,igas) = vmr(icol,ilay,igas) * col_dry_wk(icol,ilay)
603 : end do
604 : end do
605 : end do
606 : !
607 : ! ---- calculate gas optical depths ----
608 : !
609 : !$acc data copyout( jtemp, jpress, jeta, tropo, fmajor) create( col_mix, fminor)
610 : !$omp target data map(from:jtemp, jpress, jeta, tropo, fmajor) map(alloc:col_mix, fminor)
611 : call interpolation( &
612 : ncol,nlay, & ! problem dimensions
613 : ngas, nflav, neta, npres, ntemp, & ! interpolation dimensions
614 : this%flavor, &
615 : this%press_ref_log, &
616 : this%temp_ref, &
617 : this%press_ref_log_delta, &
618 : this%temp_ref_min, &
619 : this%temp_ref_delta, &
620 : this%press_ref_trop_log, &
621 : this%vmr_ref, &
622 : play, &
623 : tlay, &
624 : col_gas, &
625 : jtemp, & ! outputs
626 : fmajor,fminor,&
627 : col_mix, &
628 : tropo, &
629 1135399 : jeta,jpress)
630 1135399 : if (allocated(this%krayl)) then
631 : !$acc data copyin(this%gpoint_flavor, this%krayl) create(tau, tau_rayleigh)
632 : !$omp target data map(to:this%gpoint_flavor, this%krayl) map(alloc:tau, tau_rayleigh)
633 389263 : call zero_array(ncol, nlay, ngpt, tau)
634 : call compute_tau_absorption( &
635 : ncol,nlay,nband,ngpt, & ! dimensions
636 : ngas,nflav,neta,npres,ntemp, &
637 : nminorlower, nminorklower, & ! number of minor contributors, total num absorption coeffs
638 : nminorupper, nminorkupper, &
639 : idx_h2o, &
640 : this%gpoint_flavor, &
641 : this%get_band_lims_gpoint(), &
642 : this%kmajor, &
643 : this%kminor_lower, &
644 : this%kminor_upper, &
645 : this%minor_limits_gpt_lower, &
646 : this%minor_limits_gpt_upper, &
647 : this%minor_scales_with_density_lower, &
648 : this%minor_scales_with_density_upper, &
649 : this%scale_by_complement_lower, &
650 : this%scale_by_complement_upper, &
651 : this%idx_minor_lower, &
652 : this%idx_minor_upper, &
653 : this%idx_minor_scaling_lower, &
654 : this%idx_minor_scaling_upper, &
655 : this%kminor_start_lower, &
656 : this%kminor_start_upper, &
657 : tropo, &
658 : col_mix,fmajor,fminor, &
659 : play,tlay,col_gas, &
660 : jeta,jtemp,jpress, &
661 389263 : tau)
662 : call compute_tau_rayleigh( & !Rayleigh scattering optical depths
663 : ncol,nlay,nband,ngpt, &
664 : ngas,nflav,neta,npres,ntemp, & ! dimensions
665 : this%gpoint_flavor, &
666 : this%get_band_lims_gpoint(), &
667 : this%krayl, & ! inputs from object
668 : idx_h2o, col_dry_wk,col_gas, &
669 : fminor,jeta,tropo,jtemp, & ! local input
670 389263 : tau_rayleigh)
671 389263 : call combine_abs_and_rayleigh(tau, tau_rayleigh, optical_props)
672 : !$acc end data
673 : !$omp end target data
674 : else
675 746136 : call zero_array(ncol, nlay, ngpt, optical_props%tau)
676 : call compute_tau_absorption( &
677 : ncol,nlay,nband,ngpt, & ! dimensions
678 : ngas,nflav,neta,npres,ntemp, &
679 : nminorlower, nminorklower, & ! number of minor contributors, total num absorption coeffs
680 : nminorupper, nminorkupper, &
681 : idx_h2o, &
682 : this%gpoint_flavor, &
683 : this%get_band_lims_gpoint(), &
684 : this%kmajor, &
685 : this%kminor_lower, &
686 : this%kminor_upper, &
687 : this%minor_limits_gpt_lower, &
688 : this%minor_limits_gpt_upper, &
689 : this%minor_scales_with_density_lower, &
690 : this%minor_scales_with_density_upper, &
691 : this%scale_by_complement_lower, &
692 : this%scale_by_complement_upper, &
693 : this%idx_minor_lower, &
694 : this%idx_minor_upper, &
695 : this%idx_minor_scaling_lower, &
696 : this%idx_minor_scaling_upper, &
697 : this%kminor_start_lower, &
698 : this%kminor_start_upper, &
699 : tropo, &
700 : col_mix,fmajor,fminor, &
701 : play,tlay,col_gas, &
702 : jeta,jtemp,jpress, &
703 746136 : optical_props%tau) !
704 : select type(optical_props)
705 : type is (ty_optical_props_2str)
706 0 : call zero_array(ncol, nlay, ngpt, optical_props%ssa)
707 0 : call zero_array(ncol, nlay, ngpt, optical_props%g)
708 : type is (ty_optical_props_nstr)
709 0 : call zero_array(ncol, nlay, ngpt, optical_props%ssa)
710 : call zero_array(optical_props%get_nmom(), &
711 0 : ncol, nlay, ngpt, optical_props%p)
712 : end select
713 : end if
714 : !$acc end data
715 : !$omp end target data
716 : if (present(col_dry)) then
717 : !$acc exit data delete( col_dry)
718 : !$omp target exit data map(release:col_dry)
719 : else
720 : !$acc exit data delete( col_dry_arr)
721 : !$omp target exit data map(release:col_dry_arr)
722 : end if
723 :
724 : select type(optical_props)
725 : type is (ty_optical_props_1scl)
726 : !$acc exit data copyout( optical_props%tau)
727 : !$omp target exit data map(from:optical_props%tau)
728 : type is (ty_optical_props_2str)
729 : !$acc exit data copyout( optical_props%tau, optical_props%ssa, optical_props%g)
730 : !$omp target exit data map(from:optical_props%tau, optical_props%ssa, optical_props%g)
731 : type is (ty_optical_props_nstr)
732 : !$acc exit data copyout( optical_props%tau, optical_props%ssa, optical_props%p)
733 : !$omp target exit data map(from:optical_props%tau, optical_props%ssa, optical_props%p)
734 : end select
735 : !$acc exit data delete(optical_props)
736 :
737 : end if
738 : !$acc end data
739 : !$omp end target data
740 :
741 1135399 : end function compute_gas_taus
742 : !------------------------------------------------------------------------------------------
743 : !
744 : !> Compute the spectral solar source function adjusted to account for solar variability
745 : !> following the NRLSSI2 model of Coddington et al. 2016, doi:10.1175/BAMS-D-14-00265.1.
746 : !> as specified by the facular brightening (mg_index) and sunspot dimming (sb_index)
747 : !> indices provided as input.
748 : !>
749 : !> Users provide the NRLSSI2 facular ("Bremen") index and sunspot ("SPOT67") index.
750 : !> Changing either of these indicies will change the total solar irradiance (TSI)
751 : !> Code in extensions/mo_solar_variability may be used to compute the value of these
752 : !> indices through an average solar cycle
753 : !> Users may also specify the TSI, either alone or in conjunction with the facular and sunspot indices
754 : !
755 : !------------------------------------------------------------------------------------------
756 1536 : function set_solar_variability(this, &
757 : mg_index, sb_index, tsi) &
758 : result(error_msg)
759 : !
760 : !! Updates the spectral distribution and, optionally,
761 : !! the integrated value of the solar source function
762 : !! Modifying either index will change the total solar irradiance
763 : !
764 : class(ty_gas_optics_rrtmgp), intent(inout) :: this
765 : !
766 : real(wp), intent(in) :: mg_index !! facular brightening index (NRLSSI2 facular "Bremen" index)
767 : real(wp), intent(in) :: sb_index !! sunspot dimming index (NRLSSI2 sunspot "SPOT67" index)
768 : real(wp), optional, intent(in) :: tsi !! total solar irradiance
769 : character(len=128) :: error_msg !! Empty if successful
770 : ! ----------------------------------------------------------
771 : integer :: igpt
772 : real(wp), parameter :: a_offset = 0.1495954_wp
773 : real(wp), parameter :: b_offset = 0.00066696_wp
774 : ! ----------------------------------------------------------
775 1536 : error_msg = ""
776 1536 : if(mg_index < 0._wp) error_msg = 'mg_index out of range'
777 1536 : if(sb_index < 0._wp) error_msg = 'sb_index out of range'
778 1536 : if(error_msg /= "") return
779 : !
780 : ! Calculate solar source function for provided facular and sunspot indices
781 : !
782 : !$acc parallel loop
783 : !$omp target teams distribute parallel do simd
784 173568 : do igpt = 1, size(this%solar_source_quiet)
785 0 : this%solar_source(igpt) = this%solar_source_quiet(igpt) + &
786 0 : (mg_index - a_offset) * this%solar_source_facular(igpt) + &
787 173568 : (sb_index - b_offset) * this%solar_source_sunspot(igpt)
788 : end do
789 : !
790 : ! Scale solar source to input TSI value
791 : !
792 1536 : if (present(tsi)) error_msg = this%set_tsi(tsi)
793 :
794 : end function set_solar_variability
795 : !------------------------------------------------------------------------------------------
796 0 : function set_tsi(this, tsi) result(error_msg)
797 : !
798 : !> Scale the solar source function without changing the spectral distribution
799 : !
800 : class(ty_gas_optics_rrtmgp), intent(inout) :: this
801 : real(wp), intent(in ) :: tsi !! user-specified total solar irradiance;
802 : character(len=128) :: error_msg !! Empty if successful
803 :
804 : real(wp) :: norm
805 : ! ----------------------------------------------------------
806 0 : error_msg = ""
807 0 : if(tsi < 0._wp) then
808 0 : error_msg = 'tsi out of range'
809 : else
810 : !
811 : ! Scale the solar source function to the input tsi
812 : !
813 : !$acc kernels
814 : !$omp target
815 0 : norm = 1._wp/sum(this%solar_source(:))
816 0 : this%solar_source(:) = this%solar_source(:) * tsi * norm
817 : !$acc end kernels
818 : !$omp end target
819 : end if
820 :
821 0 : end function set_tsi
822 : !------------------------------------------------------------------------------------------
823 : !
824 : ! Compute Planck source functions at layer centers and levels
825 : !
826 746136 : function source(this, &
827 : ncol, nlay, nbnd, ngpt, &
828 746136 : play, plev, tlay, tsfc, &
829 746136 : jtemp, jpress, jeta, tropo, fmajor, &
830 : sources, & ! Planck sources
831 0 : tlev) & ! optional input
832 : result(error_msg)
833 : ! inputs
834 : class(ty_gas_optics_rrtmgp), intent(in ) :: this
835 : integer, intent(in ) :: ncol, nlay, nbnd, ngpt
836 : real(wp), dimension(ncol,nlay), intent(in ) :: play ! layer pressures [Pa, mb]
837 : real(wp), dimension(ncol,nlay+1), intent(in ) :: plev ! level pressures [Pa, mb]
838 : real(wp), dimension(ncol,nlay), intent(in ) :: tlay ! layer temperatures [K]
839 : real(wp), dimension(ncol), intent(in ) :: tsfc ! surface skin temperatures [K]
840 : ! Interplation coefficients
841 : integer, dimension(ncol,nlay), intent(in ) :: jtemp, jpress
842 : logical(wl), dimension(ncol,nlay), intent(in ) :: tropo
843 : real(wp), dimension(2,2,2,ncol,nlay,get_nflav(this)), &
844 : intent(in ) :: fmajor
845 : integer, dimension(2, ncol,nlay,get_nflav(this)), &
846 : intent(in ) :: jeta
847 : class(ty_source_func_lw ), intent(inout) :: sources
848 : real(wp), dimension(ncol,nlay+1), intent(in ), &
849 : optional, target :: tlev ! level temperatures [K]
850 : character(len=128) :: error_msg
851 : ! ----------------------------------------------------------
852 : logical(wl) :: top_at_1
853 : integer :: icol, ilay
854 : ! Variables for temperature at layer edges [K] (ncol, nlay+1)
855 1492272 : real(wp), dimension( ncol,nlay+1), target :: tlev_arr
856 746136 : real(wp), dimension(:,:), pointer :: tlev_wk
857 : ! ----------------------------------------------------------
858 746136 : error_msg = ""
859 : !
860 : ! Source function needs temperature at interfaces/levels and at layer centers
861 : !
862 746136 : if (present(tlev)) then
863 : ! Users might have provided these
864 0 : tlev_wk => tlev
865 : else
866 746136 : tlev_wk => tlev_arr
867 : !
868 : ! Interpolate temperature to levels if not provided
869 : ! Interpolation and extrapolation at boundaries is weighted by pressure
870 : !
871 12458736 : do icol = 1, ncol
872 23425200 : tlev_arr(icol,1) = tlay(icol,1) &
873 11712600 : + (plev(icol,1)-play(icol,1))*(tlay(icol,2)-tlay(icol,1)) &
874 24171336 : & / (play(icol,2)-play(icol,1))
875 : end do
876 70136784 : do ilay = 2, nlay
877 1159408584 : do icol = 1, ncol
878 3267815400 : tlev_arr(icol,ilay) = (play(icol,ilay-1)*tlay(icol,ilay-1)*(plev(icol,ilay )-play(icol,ilay)) &
879 : + play(icol,ilay )*tlay(icol,ilay )*(play(icol,ilay-1)-plev(icol,ilay))) / &
880 4426477848 : (plev(icol,ilay)*(play(icol,ilay-1) - play(icol,ilay)))
881 : end do
882 : end do
883 12458736 : do icol = 1, ncol
884 23425200 : tlev_arr(icol,nlay+1) = tlay(icol,nlay) &
885 11712600 : + (plev(icol,nlay+1)-play(icol,nlay))*(tlay(icol,nlay)-tlay(icol,nlay-1)) &
886 35883936 : / (play(icol,nlay)-play(icol,nlay-1))
887 : end do
888 : end if
889 :
890 : !-------------------------------------------------------------------
891 : ! Compute internal (Planck) source functions at layers and levels,
892 : ! which depend on mapping from spectral space that creates k-distribution.
893 : !$acc data copyin(sources) copyout( sources%lay_source, sources%lev_source_inc, sources%lev_source_dec) &
894 : !$acc copyout( sources%sfc_source, sources%sfc_source_Jac)
895 : !$omp target data map(from:sources%lay_source, sources%lev_source_inc, sources%lev_source_dec) &
896 : !$omp map(from:sources%sfc_source, sources%sfc_source_Jac)
897 :
898 : !$acc kernels copyout(top_at_1)
899 : !$omp target map(from:top_at_1)
900 746136 : top_at_1 = play(1,1) < play(1, nlay)
901 : !$acc end kernels
902 : !$omp end target
903 :
904 : call compute_Planck_source(ncol, nlay, nbnd, ngpt, &
905 : get_nflav(this), this%get_neta(), this%get_npres(), this%get_ntemp(), this%get_nPlanckTemp(), &
906 : tlay, tlev_wk, tsfc, merge(nlay, 1, top_at_1), &
907 : fmajor, jeta, tropo, jtemp, jpress, &
908 : this%get_gpoint_bands(), this%get_band_lims_gpoint(), this%planck_frac, this%temp_ref_min,&
909 : this%totplnk_delta, this%totplnk, this%gpoint_flavor, &
910 : sources%sfc_source, sources%lay_source, sources%lev_source_inc, sources%lev_source_dec, &
911 1492272 : sources%sfc_source_Jac)
912 : !$acc end data
913 : !$omp end target data
914 746136 : end function source
915 : !--------------------------------------------------------------------------------------------------------------------
916 : !
917 : ! Initialization
918 : !
919 : !--------------------------------------------------------------------------------------------------------------------
920 : ! Initialize object based on data read from netCDF file however the user desires.
921 : ! Rayleigh scattering tables may or may not be present; this is indicated with allocation status
922 : ! This interface is for the internal-sources object -- includes Plank functions and fractions
923 : !
924 1536 : function load_int(this, available_gases, gas_names, key_species, &
925 1536 : band2gpt, band_lims_wavenum, &
926 1536 : press_ref, press_ref_trop, temp_ref, &
927 1536 : temp_ref_p, temp_ref_t, vmr_ref, &
928 1536 : kmajor, kminor_lower, kminor_upper, &
929 1536 : gas_minor,identifier_minor, &
930 1536 : minor_gases_lower, minor_gases_upper, &
931 1536 : minor_limits_gpt_lower, minor_limits_gpt_upper, &
932 1536 : minor_scales_with_density_lower, &
933 1536 : minor_scales_with_density_upper, &
934 1536 : scaling_gas_lower, scaling_gas_upper, &
935 1536 : scale_by_complement_lower, &
936 1536 : scale_by_complement_upper, &
937 1536 : kminor_start_lower, &
938 1536 : kminor_start_upper, &
939 1536 : totplnk, planck_frac, &
940 : rayl_lower, rayl_upper, &
941 1536 : optimal_angle_fit) result(err_message)
942 : class(ty_gas_optics_rrtmgp), intent(inout) :: this
943 : class(ty_gas_concs), intent(in ) :: available_gases ! Which gases does the host model have available?
944 : character(len=*), dimension(:), intent(in ) :: gas_names
945 : integer, dimension(:,:,:), intent(in ) :: key_species
946 : integer, dimension(:,:), intent(in ) :: band2gpt
947 : real(wp), dimension(:,:), intent(in ) :: band_lims_wavenum
948 : real(wp), dimension(:), intent(in ) :: press_ref, temp_ref
949 : real(wp), intent(in ) :: press_ref_trop, temp_ref_p, temp_ref_t
950 : real(wp), dimension(:,:,:), intent(in ) :: vmr_ref
951 : real(wp), dimension(:,:,:,:), intent(in ) :: kmajor
952 : real(wp), dimension(:,:,:), intent(in ) :: kminor_lower, kminor_upper
953 : real(wp), dimension(:,:), intent(in ) :: totplnk
954 : real(wp), dimension(:,:,:,:), intent(in ) :: planck_frac
955 : real(wp), dimension(:,:,:), intent(in ), &
956 : allocatable :: rayl_lower, rayl_upper
957 : real(wp), dimension(:,:), intent(in ) :: optimal_angle_fit
958 : character(len=*), dimension(:), intent(in ) :: gas_minor,identifier_minor
959 : character(len=*), dimension(:), intent(in ) :: minor_gases_lower, &
960 : minor_gases_upper
961 : integer, dimension(:,:), intent(in ) :: minor_limits_gpt_lower, &
962 : minor_limits_gpt_upper
963 : logical(wl), dimension(:), intent(in ) :: minor_scales_with_density_lower, &
964 : minor_scales_with_density_upper
965 : character(len=*), dimension(:), intent(in ) :: scaling_gas_lower, &
966 : scaling_gas_upper
967 : logical(wl), dimension(:), intent(in ) :: scale_by_complement_lower,&
968 : scale_by_complement_upper
969 : integer, dimension(:), intent(in ) :: kminor_start_lower,&
970 : kminor_start_upper
971 : character(len = 128) :: err_message
972 : ! ----
973 : !$acc enter data copyin(this)
974 1536 : call this%finalize()
975 : err_message = init_abs_coeffs(this, &
976 : available_gases, &
977 : gas_names, key_species, &
978 : band2gpt, band_lims_wavenum, &
979 : press_ref, temp_ref, &
980 : press_ref_trop, temp_ref_p, temp_ref_t, &
981 : vmr_ref, &
982 : kmajor, kminor_lower, kminor_upper, &
983 : gas_minor,identifier_minor,&
984 : minor_gases_lower, minor_gases_upper, &
985 : minor_limits_gpt_lower, &
986 : minor_limits_gpt_upper, &
987 : minor_scales_with_density_lower, &
988 : minor_scales_with_density_upper, &
989 : scaling_gas_lower, scaling_gas_upper, &
990 : scale_by_complement_lower, &
991 : scale_by_complement_upper, &
992 : kminor_start_lower, &
993 : kminor_start_upper, &
994 1536 : rayl_lower, rayl_upper)
995 : ! Planck function tables
996 : !
997 0 : allocate(this%totplnk (size(totplnk, 1), size(totplnk, 2)), &
998 0 : this%planck_frac (size(planck_frac,4), size(planck_frac,2),size(planck_frac,3), size(planck_frac,1)), &
999 18432 : this%optimal_angle_fit(size(optimal_angle_fit, 1), size(optimal_angle_fit, 2)))
1000 4844544 : this%totplnk = totplnk
1001 : ! this%planck_frac = planck_frac
1002 : this%planck_frac = RESHAPE(planck_frac,(/size(planck_frac,4), size(planck_frac,2), &
1003 1604527104 : size(planck_frac,3), size(planck_frac,1)/),ORDER =(/4,2,3,1/))
1004 76800 : this%optimal_angle_fit = optimal_angle_fit
1005 : !$acc enter data copyin(this%totplnk, this%planck_frac, this%optimal_angle_fit)
1006 : !$omp target enter data map(to:this%totplnk, this%planck_frac, this%optimal_angle_fit)
1007 :
1008 : ! Temperature steps for Planck function interpolation
1009 : ! Assumes that temperature minimum and max are the same for the absorption coefficient grid and the
1010 : ! Planck grid and the Planck grid is equally spaced
1011 1536 : this%totplnk_delta = (this%temp_ref_max-this%temp_ref_min) / (size(this%totplnk,dim=1)-1)
1012 1536 : end function load_int
1013 :
1014 : !--------------------------------------------------------------------------------------------------------------------
1015 : !
1016 : ! Initialize object based on data read from netCDF file however the user desires.
1017 : ! Rayleigh scattering tables may or may not be present; this is indicated with allocation status
1018 : ! This interface is for the external-sources object -- includes TOA source function table
1019 : !
1020 1536 : function load_ext(this, available_gases, gas_names, key_species, &
1021 1536 : band2gpt, band_lims_wavenum, &
1022 1536 : press_ref, press_ref_trop, temp_ref, &
1023 1536 : temp_ref_p, temp_ref_t, vmr_ref, &
1024 1536 : kmajor, kminor_lower, kminor_upper, &
1025 1536 : gas_minor,identifier_minor, &
1026 1536 : minor_gases_lower, minor_gases_upper, &
1027 1536 : minor_limits_gpt_lower, minor_limits_gpt_upper, &
1028 1536 : minor_scales_with_density_lower, &
1029 1536 : minor_scales_with_density_upper, &
1030 1536 : scaling_gas_lower, scaling_gas_upper, &
1031 1536 : scale_by_complement_lower, &
1032 1536 : scale_by_complement_upper, &
1033 1536 : kminor_start_lower, &
1034 1536 : kminor_start_upper, &
1035 1536 : solar_quiet, solar_facular, solar_sunspot, &
1036 : tsi_default, mg_default, sb_default, &
1037 : rayl_lower, rayl_upper) result(err_message)
1038 : class(ty_gas_optics_rrtmgp), intent(inout) :: this
1039 : class(ty_gas_concs), intent(in ) :: available_gases ! Which gases does the host model have available?
1040 : character(len=*), &
1041 : dimension(:), intent(in) :: gas_names
1042 : integer, dimension(:,:,:), intent(in) :: key_species
1043 : integer, dimension(:,:), intent(in) :: band2gpt
1044 : real(wp), dimension(:,:), intent(in) :: band_lims_wavenum
1045 : real(wp), dimension(:), intent(in) :: press_ref, temp_ref
1046 : real(wp), intent(in) :: press_ref_trop, temp_ref_p, temp_ref_t
1047 : real(wp), dimension(:,:,:), intent(in) :: vmr_ref
1048 : real(wp), dimension(:,:,:,:), intent(in) :: kmajor
1049 : real(wp), dimension(:,:,:), intent(in) :: kminor_lower, kminor_upper
1050 : character(len=*), dimension(:), &
1051 : intent(in) :: gas_minor, &
1052 : identifier_minor
1053 : character(len=*), dimension(:), &
1054 : intent(in) :: minor_gases_lower, &
1055 : minor_gases_upper
1056 : integer, dimension(:,:), intent(in) :: &
1057 : minor_limits_gpt_lower, &
1058 : minor_limits_gpt_upper
1059 : logical(wl), dimension(:), intent(in) :: &
1060 : minor_scales_with_density_lower, &
1061 : minor_scales_with_density_upper
1062 : character(len=*),dimension(:),intent(in) :: &
1063 : scaling_gas_lower, &
1064 : scaling_gas_upper
1065 : logical(wl), dimension(:), intent(in) :: &
1066 : scale_by_complement_lower, &
1067 : scale_by_complement_upper
1068 : integer, dimension(:), intent(in) :: &
1069 : kminor_start_lower, &
1070 : kminor_start_upper
1071 : real(wp), dimension(:), intent(in) :: solar_quiet, &
1072 : solar_facular, &
1073 : solar_sunspot
1074 : real(wp), intent(in) :: tsi_default, &
1075 : mg_default, sb_default
1076 : real(wp), dimension(:,:,:), intent(in), &
1077 : allocatable :: rayl_lower, rayl_upper
1078 : character(len = 128) err_message
1079 :
1080 : integer :: ngpt
1081 : ! ----
1082 : !$acc enter data copyin(this)
1083 1536 : call this%finalize()
1084 : err_message = init_abs_coeffs(this, &
1085 : available_gases, &
1086 : gas_names, key_species, &
1087 : band2gpt, band_lims_wavenum, &
1088 : press_ref, temp_ref, &
1089 : press_ref_trop, temp_ref_p, temp_ref_t, &
1090 : vmr_ref, &
1091 : kmajor, kminor_lower, kminor_upper, &
1092 : gas_minor,identifier_minor, &
1093 : minor_gases_lower, minor_gases_upper, &
1094 : minor_limits_gpt_lower, &
1095 : minor_limits_gpt_upper, &
1096 : minor_scales_with_density_lower, &
1097 : minor_scales_with_density_upper, &
1098 : scaling_gas_lower, scaling_gas_upper, &
1099 : scale_by_complement_lower, &
1100 : scale_by_complement_upper, &
1101 : kminor_start_lower, &
1102 : kminor_start_upper, &
1103 1536 : rayl_lower, rayl_upper)
1104 1536 : if(err_message == "") then
1105 : !
1106 : ! Spectral solar irradiance terms init
1107 : !
1108 1536 : ngpt = size(solar_quiet)
1109 0 : allocate(this%solar_source_quiet(ngpt), this%solar_source_facular(ngpt), &
1110 9216 : this%solar_source_sunspot(ngpt), this%solar_source(ngpt))
1111 : !$acc enter data create( this%solar_source_quiet, this%solar_source_facular, this%solar_source_sunspot, this%solar_source)
1112 : !$omp target enter data map(alloc:this%solar_source_quiet, this%solar_source_facular, this%solar_source_sunspot, this%solar_source)
1113 : !$acc kernels
1114 : !$omp target
1115 175104 : this%solar_source_quiet = solar_quiet
1116 175104 : this%solar_source_facular = solar_facular
1117 175104 : this%solar_source_sunspot = solar_sunspot
1118 : !$acc end kernels
1119 : !$omp end target
1120 1536 : err_message = this%set_solar_variability(mg_default, sb_default)
1121 : endif
1122 1536 : end function load_ext
1123 : !--------------------------------------------------------------------------------------------------------------------
1124 : !
1125 : ! Initialize absorption coefficient arrays,
1126 : ! including Rayleigh scattering tables if provided (allocated)
1127 : !
1128 3072 : function init_abs_coeffs(this, &
1129 : available_gases, &
1130 12288 : gas_names, key_species, &
1131 3072 : band2gpt, band_lims_wavenum, &
1132 3072 : press_ref, temp_ref, &
1133 : press_ref_trop, temp_ref_p, temp_ref_t, &
1134 3072 : vmr_ref, &
1135 3072 : kmajor, kminor_lower, kminor_upper, &
1136 3072 : gas_minor,identifier_minor,&
1137 6144 : minor_gases_lower, minor_gases_upper, &
1138 3072 : minor_limits_gpt_lower, &
1139 3072 : minor_limits_gpt_upper, &
1140 3072 : minor_scales_with_density_lower, &
1141 3072 : minor_scales_with_density_upper, &
1142 6144 : scaling_gas_lower, scaling_gas_upper, &
1143 3072 : scale_by_complement_lower, &
1144 3072 : scale_by_complement_upper, &
1145 3072 : kminor_start_lower, &
1146 3072 : kminor_start_upper, &
1147 : rayl_lower, rayl_upper) result(err_message)
1148 : class(ty_gas_optics_rrtmgp), intent(inout) :: this
1149 : class(ty_gas_concs), intent(in ) :: available_gases
1150 : character(len=*), &
1151 : dimension(:), intent(in) :: gas_names
1152 : integer, dimension(:,:,:), intent(in) :: key_species
1153 : integer, dimension(:,:), intent(in) :: band2gpt
1154 : real(wp), dimension(:,:), intent(in) :: band_lims_wavenum
1155 : real(wp), dimension(:), intent(in) :: press_ref, temp_ref
1156 : real(wp), intent(in) :: press_ref_trop, temp_ref_p, temp_ref_t
1157 : real(wp), dimension(:,:,:), intent(in) :: vmr_ref
1158 : real(wp), dimension(:,:,:,:), intent(in) :: kmajor
1159 : real(wp), dimension(:,:,:), intent(in) :: kminor_lower, kminor_upper
1160 : character(len=*), dimension(:), &
1161 : intent(in) :: gas_minor, &
1162 : identifier_minor
1163 : character(len=*), dimension(:), &
1164 : intent(in) :: minor_gases_lower, &
1165 : minor_gases_upper
1166 : integer, dimension(:,:), intent(in) :: minor_limits_gpt_lower, &
1167 : minor_limits_gpt_upper
1168 : logical(wl), dimension(:), intent(in) :: minor_scales_with_density_lower, &
1169 : minor_scales_with_density_upper
1170 : character(len=*), dimension(:),&
1171 : intent(in) :: scaling_gas_lower, &
1172 : scaling_gas_upper
1173 : logical(wl), dimension(:), intent(in) :: scale_by_complement_lower, &
1174 : scale_by_complement_upper
1175 : integer, dimension(:), intent(in) :: kminor_start_lower, &
1176 : kminor_start_upper
1177 : real(wp), dimension(:,:,:), intent(in), &
1178 : allocatable :: rayl_lower, rayl_upper
1179 : character(len=128) :: err_message
1180 : ! --------------------------------------------------------------------------
1181 3072 : logical, dimension(:), allocatable :: gas_is_present
1182 3072 : logical, dimension(:), allocatable :: key_species_present_init
1183 3072 : integer, dimension(:,:,:), allocatable :: key_species_red
1184 3072 : real(wp), dimension(:,:,:), allocatable :: vmr_ref_red
1185 : character(len=256), &
1186 3072 : dimension(:), allocatable :: minor_gases_lower_red, &
1187 3072 : minor_gases_upper_red
1188 : character(len=256), &
1189 3072 : dimension(:), allocatable :: scaling_gas_lower_red, &
1190 3072 : scaling_gas_upper_red
1191 : integer :: i, j, idx
1192 : integer :: ngas
1193 : ! --------------------------------------
1194 3072 : err_message = this%ty_optical_props%init(band_lims_wavenum, band2gpt)
1195 3072 : if(len_trim(err_message) /= 0) return
1196 : !
1197 : ! Which gases known to the gas optics are present in the host model (available_gases)?
1198 : !
1199 3072 : ngas = size(gas_names)
1200 9216 : allocate(gas_is_present(ngas))
1201 61440 : do i = 1, ngas
1202 : ! Next line causes a compiler bug in gfortran 11.0.1 on Mac ARM
1203 : ! Should replace gas_names with get_gas_names() and make gas_names private in ty_gas_concs
1204 61440 : gas_is_present(i) = string_in_array(gas_names(i), available_gases%gas_names)
1205 : end do
1206 : !
1207 : ! Now the number of gases is the union of those known to the k-distribution and provided
1208 : ! by the host model
1209 : !
1210 61440 : ngas = count(gas_is_present)
1211 : !
1212 : ! Initialize the gas optics object, keeping only those gases known to the
1213 : ! gas optics and also present in the host model
1214 : !
1215 30720 : this%gas_names = pack(gas_names,mask=gas_is_present)
1216 : ! Copy-ins below
1217 :
1218 : allocate(vmr_ref_red(size(vmr_ref,dim=1),0:ngas, &
1219 15360 : size(vmr_ref,dim=3)))
1220 : ! Gas 0 is used in single-key species method, set to 1.0 (col_dry)
1221 132096 : vmr_ref_red(:,0,:) = vmr_ref(:,1,:)
1222 27648 : do i = 1, ngas
1223 24576 : idx = string_loc_in_array(this%gas_names(i), gas_names)
1224 1059840 : vmr_ref_red(:,i,:) = vmr_ref(:,idx+1,:)
1225 : enddo
1226 3072 : call move_alloc(vmr_ref_red, this%vmr_ref)
1227 : !$acc enter data copyin(this%vmr_ref, this%gas_names)
1228 : !$omp target enter data map(to:this%vmr_ref, this%gas_names)
1229 : !
1230 : ! Reduce minor arrays so variables only contain minor gases that are available
1231 : ! Reduce size of minor Arrays
1232 : !
1233 : call reduce_minor_arrays(available_gases, &
1234 : gas_minor,identifier_minor, &
1235 : kminor_lower, &
1236 : minor_gases_lower, &
1237 : minor_limits_gpt_lower, &
1238 : minor_scales_with_density_lower, &
1239 : scaling_gas_lower, &
1240 : scale_by_complement_lower, &
1241 : kminor_start_lower, &
1242 : this%kminor_lower, &
1243 : minor_gases_lower_red, &
1244 : this%minor_limits_gpt_lower, &
1245 : this%minor_scales_with_density_lower, &
1246 : scaling_gas_lower_red, &
1247 : this%scale_by_complement_lower, &
1248 3072 : this%kminor_start_lower)
1249 : call reduce_minor_arrays(available_gases, &
1250 : gas_minor,identifier_minor,&
1251 : kminor_upper, &
1252 : minor_gases_upper, &
1253 : minor_limits_gpt_upper, &
1254 : minor_scales_with_density_upper, &
1255 : scaling_gas_upper, &
1256 : scale_by_complement_upper, &
1257 : kminor_start_upper, &
1258 : this%kminor_upper, &
1259 : minor_gases_upper_red, &
1260 : this%minor_limits_gpt_upper, &
1261 : this%minor_scales_with_density_upper, &
1262 : scaling_gas_upper_red, &
1263 : this%scale_by_complement_upper, &
1264 3072 : this%kminor_start_upper)
1265 : !$acc enter data copyin(this%minor_limits_gpt_lower, this%minor_limits_gpt_upper)
1266 : !$omp target enter data map(to:this%minor_limits_gpt_lower, this%minor_limits_gpt_upper)
1267 : !$acc enter data copyin(this%minor_scales_with_density_lower, this%minor_scales_with_density_upper)
1268 : !$omp target enter data map(to:this%minor_scales_with_density_lower, this%minor_scales_with_density_upper)
1269 : !$acc enter data copyin(this%scale_by_complement_lower, this%scale_by_complement_upper)
1270 : !$omp target enter data map(to:this%scale_by_complement_lower, this%scale_by_complement_upper)
1271 : !$acc enter data copyin(this%kminor_start_lower, this%kminor_start_upper)
1272 : !$omp target enter data map(to:this%kminor_start_lower, this%kminor_start_upper)
1273 : !$acc enter data copyin(this%kminor_lower, this%kminor_upper)
1274 : !$omp target enter data map(to:this%kminor_lower, this%kminor_upper)
1275 :
1276 : ! Arrays not reduced by the presence, or lack thereof, of a gas
1277 0 : allocate(this%press_ref(size(press_ref)), this%temp_ref(size(temp_ref)), &
1278 30720 : this%kmajor(size(kmajor,4),size(kmajor,2),size(kmajor,3),size(kmajor,1)))
1279 184320 : this%press_ref(:) = press_ref(:)
1280 46080 : this%temp_ref(:) = temp_ref(:)
1281 3008489472 : this%kmajor = RESHAPE(kmajor,(/size(kmajor,4),size(kmajor,2),size(kmajor,3),size(kmajor,1)/), ORDER= (/4,2,3,1/))
1282 : !$acc enter data copyin(this%press_ref, this%temp_ref, this%kmajor)
1283 : !$omp target enter data map(to:this%press_ref, this%temp_ref, this%kmajor)
1284 :
1285 :
1286 3072 : if(allocated(rayl_lower) .neqv. allocated(rayl_upper)) then
1287 0 : err_message = "rayl_lower and rayl_upper must have the same allocation status"
1288 0 : return
1289 : end if
1290 3072 : if (allocated(rayl_lower)) then
1291 9216 : allocate(this%krayl(size(rayl_lower,dim=3),size(rayl_lower,dim=2),size(rayl_lower,dim=1),2))
1292 0 : this%krayl(:,:,:,1) = RESHAPE(rayl_lower,(/size(rayl_lower,dim=3),size(rayl_lower,dim=2), &
1293 23402496 : size(rayl_lower,dim=1)/),ORDER =(/3,2,1/))
1294 0 : this%krayl(:,:,:,2) = RESHAPE(rayl_upper,(/size(rayl_lower,dim=3),size(rayl_lower,dim=2), &
1295 23402496 : size(rayl_lower,dim=1)/),ORDER =(/3,2,1/))
1296 : !$acc enter data copyin(this%krayl)
1297 : !$omp target enter data map(to:this%krayl)
1298 : end if
1299 :
1300 : ! ---- post processing ----
1301 : ! creates log reference pressure
1302 9216 : allocate(this%press_ref_log(size(this%press_ref)))
1303 184320 : this%press_ref_log(:) = log(this%press_ref(:))
1304 : !$acc enter data copyin(this%press_ref_log)
1305 : !$omp target enter data map(to:this%press_ref_log)
1306 :
1307 : ! log scale of reference pressure
1308 3072 : this%press_ref_trop_log = log(press_ref_trop)
1309 :
1310 : ! Get index of gas (if present) for determining col_gas
1311 3072 : call create_idx_minor(this%gas_names, gas_minor, identifier_minor, minor_gases_lower_red, this%idx_minor_lower)
1312 3072 : call create_idx_minor(this%gas_names, gas_minor, identifier_minor, minor_gases_upper_red, this%idx_minor_upper)
1313 : ! Get index of gas (if present) that has special treatment in density scaling
1314 3072 : call create_idx_minor_scaling(this%gas_names, scaling_gas_lower_red, this%idx_minor_scaling_lower)
1315 3072 : call create_idx_minor_scaling(this%gas_names, scaling_gas_upper_red, this%idx_minor_scaling_upper)
1316 : !$acc enter data copyin(this%idx_minor_lower, this%idx_minor_upper)
1317 : !$omp target enter data map(to:this%idx_minor_lower, this%idx_minor_upper)
1318 : !$acc enter data copyin(this%idx_minor_scaling_lower, this%idx_minor_scaling_upper)
1319 : !$omp target enter data map(to:this%idx_minor_scaling_lower, this%idx_minor_scaling_upper)
1320 :
1321 : ! create flavor list
1322 : ! Reduce (remap) key_species list; checks that all key gases are present in incoming
1323 : call create_key_species_reduce(gas_names,this%gas_names, &
1324 3072 : key_species,key_species_red,key_species_present_init)
1325 3072 : err_message = check_key_species_present_init(gas_names,key_species_present_init)
1326 3072 : if(len_trim(err_message) /= 0) return
1327 : ! create flavor list
1328 3072 : call create_flavor(key_species_red, this%flavor)
1329 : ! create gpoint_flavor list
1330 3072 : call create_gpoint_flavor(key_species_red, this%get_gpoint_bands(), this%flavor, this%gpoint_flavor)
1331 : !Copy-ins at end of subroutine
1332 :
1333 : ! minimum, maximum reference temperature, pressure -- assumes low-to-high ordering
1334 : ! for T, high-to-low ordering for p
1335 3072 : this%temp_ref_min = this%temp_ref (1)
1336 3072 : this%temp_ref_max = this%temp_ref (size(this%temp_ref))
1337 3072 : this%press_ref_min = this%press_ref(size(this%press_ref))
1338 3072 : this%press_ref_max = this%press_ref(1)
1339 :
1340 : ! creates press_ref_log, temp_ref_delta
1341 3072 : this%press_ref_log_delta = (log(this%press_ref_min)-log(this%press_ref_max))/(size(this%press_ref)-1)
1342 3072 : this%temp_ref_delta = (this%temp_ref_max-this%temp_ref_min)/(size(this%temp_ref)-1)
1343 :
1344 : ! Which species are key in one or more bands?
1345 : ! this%flavor is an index into this%gas_names
1346 : !
1347 3072 : if (allocated(this%is_key)) deallocate(this%is_key) ! Shouldn't ever happen...
1348 9216 : allocate(this%is_key(this%get_ngas()))
1349 27648 : this%is_key(:) = .False.
1350 32256 : do j = 1, size(this%flavor, 2)
1351 90624 : do i = 1, size(this%flavor, 1) ! extents should be 2
1352 87552 : if (this%flavor(i,j) /= 0) this%is_key(this%flavor(i,j)) = .true.
1353 : end do
1354 : end do
1355 : !$acc enter data copyin(this%flavor, this%gpoint_flavor, this%is_key)
1356 : !$omp target enter data map(to:this%flavor, this%gpoint_flavor, this%is_key)
1357 :
1358 12288 : end function init_abs_coeffs
1359 : ! ----------------------------------------------------------------------------------------------------
1360 3072 : function check_key_species_present_init(gas_names, key_species_present_init) result(err_message)
1361 : logical, dimension(:), intent(in) :: key_species_present_init
1362 : character(len=*), dimension(:), intent(in) :: gas_names
1363 : character(len=128) :: err_message
1364 :
1365 : integer :: i
1366 :
1367 3072 : err_message=''
1368 61440 : do i = 1, size(key_species_present_init)
1369 58368 : if(.not. key_species_present_init(i)) &
1370 3072 : err_message = ' ' // trim(gas_names(i)) // trim(err_message)
1371 : end do
1372 3072 : if(len_trim(err_message) > 0) err_message = "gas_optics: required gases" // trim(err_message) // " are not provided"
1373 :
1374 3072 : end function check_key_species_present_init
1375 : !------------------------------------------------------------------------------------------
1376 : !
1377 : ! Ensure that every key gas required by the k-distribution is
1378 : ! present in the gas concentration object
1379 : !
1380 1135399 : function check_key_species_present(this, gas_desc) result(error_msg)
1381 : class(ty_gas_optics_rrtmgp), intent(in) :: this
1382 : class(ty_gas_concs), intent(in) :: gas_desc
1383 : character(len=128) :: error_msg
1384 :
1385 : ! Local variables
1386 10218591 : character(len=32), dimension(count(this%is_key(:) )) :: key_gas_names
1387 : integer :: igas
1388 : ! --------------------------------------
1389 1135399 : error_msg = ""
1390 1135399 : key_gas_names = pack(this%gas_names, mask=this%is_key)
1391 6812394 : do igas = 1, size(key_gas_names)
1392 : ! Next line causes a compiler bug in gfortran 11.0.1 on Mac ARM
1393 : ! Should replace gas_names with get_gas_names() and make gas_names private in ty_gas_concs
1394 5676995 : if(.not. string_in_array(key_gas_names(igas), gas_desc%gas_names)) &
1395 1135399 : error_msg = ' ' // trim(lower_case(key_gas_names(igas))) // trim(error_msg)
1396 : end do
1397 1135399 : if(len_trim(error_msg) > 0) error_msg = "gas_optics: required gases" // trim(error_msg) // " are not provided"
1398 :
1399 1135399 : end function check_key_species_present
1400 : !--------------------------------------------------------------------------------------------------------------------
1401 : !
1402 : ! Inquiry functions
1403 : !
1404 : !--------------------------------------------------------------------------------------------------------------------
1405 : !
1406 : !> return true if initialized for internal sources/longwave, false otherwise
1407 : !
1408 0 : pure function source_is_internal(this)
1409 : class(ty_gas_optics_rrtmgp), intent(in) :: this
1410 : logical :: source_is_internal
1411 0 : source_is_internal = allocated(this%totplnk) .and. allocated(this%planck_frac)
1412 0 : end function source_is_internal
1413 : !--------------------------------------------------------------------------------------------------------------------
1414 : !
1415 : !> return true if initialized for external sources/shortwave, false otherwise
1416 : !
1417 0 : pure function source_is_external(this)
1418 : class(ty_gas_optics_rrtmgp), intent(in) :: this
1419 : logical :: source_is_external
1420 0 : source_is_external = allocated(this%solar_source)
1421 0 : end function source_is_external
1422 :
1423 : !--------------------------------------------------------------------------------------------------------------------
1424 : !
1425 : !> return the names of the gases known to the k-distributions
1426 : !
1427 0 : pure function get_gases(this)
1428 : class(ty_gas_optics_rrtmgp), intent(in) :: this
1429 : character(32), dimension(get_ngas(this)) :: get_gases !! names of the gases known to the k-distributions
1430 :
1431 0 : get_gases = this%gas_names
1432 : end function get_gases
1433 : !--------------------------------------------------------------------------------------------------------------------
1434 : !
1435 : !> return the minimum pressure on the interpolation grids
1436 : !
1437 746136 : pure function get_press_min(this)
1438 : class(ty_gas_optics_rrtmgp), intent(in) :: this
1439 : real(wp) :: get_press_min !! minimum pressure for which the k-dsitribution is valid
1440 :
1441 746136 : get_press_min = this%press_ref_min
1442 746136 : end function get_press_min
1443 :
1444 : !--------------------------------------------------------------------------------------------------------------------
1445 : !
1446 : !> return the maximum pressure on the interpolation grids
1447 : !
1448 0 : pure function get_press_max(this)
1449 : class(ty_gas_optics_rrtmgp), intent(in) :: this
1450 : real(wp) :: get_press_max !! maximum pressure for which the k-dsitribution is valid
1451 :
1452 0 : get_press_max = this%press_ref_max
1453 0 : end function get_press_max
1454 :
1455 : !--------------------------------------------------------------------------------------------------------------------
1456 : !
1457 : !> return the minimum temparature on the interpolation grids
1458 : !
1459 746136 : pure function get_temp_min(this)
1460 : class(ty_gas_optics_rrtmgp), intent(in) :: this
1461 : real(wp) :: get_temp_min !! minimum temperature for which the k-dsitribution is valid
1462 :
1463 746136 : get_temp_min = this%temp_ref_min
1464 746136 : end function get_temp_min
1465 :
1466 : !--------------------------------------------------------------------------------------------------------------------
1467 : !
1468 : !> return the maximum temparature on the interpolation grids
1469 : !
1470 746136 : pure function get_temp_max(this)
1471 : class(ty_gas_optics_rrtmgp), intent(in) :: this
1472 : real(wp) :: get_temp_max !! maximum temperature for which the k-dsitribution is valid
1473 :
1474 746136 : get_temp_max = this%temp_ref_max
1475 746136 : end function get_temp_max
1476 : !--------------------------------------------------------------------------------------------------------------------
1477 : !
1478 : !> Utility function, provided for user convenience
1479 : !> computes column amounts of dry air using hydrostatic equation
1480 : !
1481 1135399 : function get_col_dry(vmr_h2o, plev, latitude) result(col_dry)
1482 : ! input
1483 : real(wp), dimension(:,:), intent(in) :: vmr_h2o ! volume mixing ratio of water vapor to dry air; (ncol,nlay)
1484 : real(wp), dimension(:,:), intent(in) :: plev ! Layer boundary pressures [Pa] (ncol,nlay+1)
1485 : real(wp), dimension(:), optional, &
1486 : intent(in) :: latitude ! Latitude [degrees] (ncol)
1487 : ! output
1488 : real(wp), dimension(size(plev,dim=1),size(plev,dim=2)-1) :: col_dry ! Column dry amount (ncol,nlay)
1489 : ! ------------------------------------------------
1490 : ! first and second term of Helmert formula
1491 : real(wp), parameter :: helmert1 = 9.80665_wp
1492 : real(wp), parameter :: helmert2 = 0.02586_wp
1493 : ! local variables
1494 2270798 : real(wp), dimension(size(plev,dim=1)) :: g0 ! (ncol)
1495 : real(wp):: delta_plev, m_air, fact
1496 : integer :: ncol, nlev
1497 : integer :: icol, ilev ! nlay = nlev-1
1498 : ! ------------------------------------------------
1499 1135399 : ncol = size(plev, dim=1)
1500 1135399 : nlev = size(plev, dim=2)
1501 : !$acc data create(g0)
1502 : !$omp target data map(alloc:g0)
1503 1135399 : if(present(latitude)) then
1504 : ! A purely OpenACC implementation would probably compute g0 within the kernel below
1505 : !$acc parallel loop
1506 : !$omp target teams distribute parallel do simd
1507 0 : do icol = 1, ncol
1508 0 : g0(icol) = helmert1 - helmert2 * cos(2.0_wp * pi * latitude(icol) / 180.0_wp) ! acceleration due to gravity [m/s^2]
1509 : end do
1510 : else
1511 : !$acc parallel loop
1512 : !$omp target teams distribute parallel do simd
1513 18704299 : do icol = 1, ncol
1514 18704299 : g0(icol) = grav
1515 : end do
1516 : end if
1517 :
1518 : !$acc parallel loop gang vector collapse(2) copyin(plev,vmr_h2o) copyout(col_dry)
1519 : !$omp target teams distribute parallel do simd collapse(2) map(to:plev,vmr_h2o) map(from:col_dry)
1520 107862905 : do ilev = 1, nlev-1
1521 1759339505 : do icol = 1, ncol
1522 1651476600 : delta_plev = abs(plev(icol,ilev) - plev(icol,ilev+1))
1523 : ! Get average mass of moist air per mole of moist air
1524 1651476600 : fact = 1._wp / (1.+vmr_h2o(icol,ilev))
1525 1651476600 : m_air = (m_dry + m_h2o * vmr_h2o(icol,ilev)) * fact
1526 1758204106 : col_dry(icol,ilev) = 10._wp * delta_plev * avogad * fact/(1000._wp*m_air*100._wp*g0(icol))
1527 : end do
1528 : end do
1529 : !$acc end data
1530 : !$omp end target data
1531 1135399 : end function get_col_dry
1532 : !--------------------------------------------------------------------------------------------------------------------
1533 : !
1534 : !> Compute a transport angle that minimizes flux errors at surface and TOA based on empirical fits
1535 : !
1536 0 : function compute_optimal_angles(this, optical_props, optimal_angles) result(err_msg)
1537 : ! input
1538 : class(ty_gas_optics_rrtmgp), intent(in ) :: this
1539 : class(ty_optical_props_arry), intent(in ) :: optical_props !! Optical properties
1540 : real(wp), dimension(:,:), intent( out) :: optimal_angles !! Secant of optical transport angle
1541 : character(len=128) :: err_msg !! Empty if successful
1542 : !----------------------------
1543 : integer :: ncol, nlay, ngpt
1544 : integer :: icol, ilay, igpt, bnd
1545 : real(wp) :: t, trans_total
1546 : #if defined _CRAYFTN && _RELEASE_MAJOR == 14 && _RELEASE_MINOR == 0 && _RELEASE_PATCHLEVEL == 3
1547 : # define CRAY_WORKAROUND
1548 : #endif
1549 : #ifdef CRAY_WORKAROUND
1550 : integer, allocatable :: bands(:)
1551 : #else
1552 0 : integer :: bands(optical_props%get_ngpt())
1553 : #endif
1554 : !----------------------------
1555 0 : ncol = optical_props%get_ncol()
1556 0 : nlay = optical_props%get_nlay()
1557 0 : ngpt = optical_props%get_ngpt()
1558 : #ifdef CRAY_WORKAROUND
1559 : allocate( bands(ngpt) ) ! In order to work with CCE 14 (it is also better software)
1560 : #endif
1561 :
1562 0 : err_msg=""
1563 0 : if(.not. this%gpoints_are_equal(optical_props)) &
1564 0 : err_msg = "gas_optics%compute_optimal_angles: optical_props has different spectral discretization than gas_optics"
1565 0 : if(.not. extents_are(optimal_angles, ncol, ngpt)) &
1566 0 : err_msg = "gas_optics%compute_optimal_angles: optimal_angles different dimension (ncol)"
1567 0 : if (err_msg /= "") return
1568 :
1569 0 : do igpt = 1, ngpt
1570 0 : bands(igpt) = optical_props%convert_gpt2band(igpt)
1571 : enddo
1572 : !
1573 : ! column transmissivity
1574 : !
1575 : !$acc parallel loop gang vector collapse(2) copyin(bands, optical_props, optical_props%tau) copyout(optimal_angles)
1576 : !$omp target teams distribute parallel do simd collapse(2) map(to:bands, optical_props%tau) map(from:optimal_angles)
1577 0 : do icol = 1, ncol
1578 0 : do igpt = 1, ngpt
1579 : !
1580 : ! Column transmissivity
1581 : !
1582 : t = 0._wp
1583 0 : trans_total = 0._wp
1584 0 : do ilay = 1, nlay
1585 0 : t = t + optical_props%tau(icol,ilay,igpt)
1586 : end do
1587 0 : trans_total = exp(-t)
1588 : !
1589 : ! Optimal transport angle is a linear fit to column transmissivity
1590 : !
1591 0 : optimal_angles(icol,igpt) = this%optimal_angle_fit(1,bands(igpt))*trans_total + &
1592 0 : this%optimal_angle_fit(2,bands(igpt))
1593 : end do
1594 : end do
1595 : end function compute_optimal_angles
1596 : !--------------------------------------------------------------------------------------------------------------------
1597 : !
1598 : ! Internal procedures
1599 : !
1600 : !--------------------------------------------------------------------------------------------------------------------
1601 829440 : pure function rewrite_key_species_pair(key_species_pair)
1602 : ! (0,0) becomes (2,2) -- because absorption coefficients for these g-points will be 0.
1603 : integer, dimension(2) :: rewrite_key_species_pair
1604 : integer, dimension(2), intent(in) :: key_species_pair
1605 2488320 : rewrite_key_species_pair = key_species_pair
1606 943104 : if (all(key_species_pair(:).eq.(/0,0/))) then
1607 170496 : rewrite_key_species_pair(:) = (/2,2/)
1608 : end if
1609 : end function
1610 :
1611 : ! ---------------------------------------------------------------------------------------
1612 : ! true is key_species_pair exists in key_species_list
1613 184320 : pure function key_species_pair_exists(key_species_list, key_species_pair)
1614 : logical :: key_species_pair_exists
1615 : integer, dimension(:,:), intent(in) :: key_species_list
1616 : integer, dimension(2), intent(in) :: key_species_pair
1617 : integer :: i
1618 1892352 : do i=1,size(key_species_list,dim=2)
1619 2697216 : if (all(key_species_list(:,i).eq.key_species_pair(:))) then
1620 125952 : key_species_pair_exists = .true.
1621 125952 : return
1622 : end if
1623 : end do
1624 58368 : key_species_pair_exists = .false.
1625 58368 : end function key_species_pair_exists
1626 : ! ---------------------------------------------------------------------------------------
1627 : ! create flavor list --
1628 : ! an unordered array of extent (2,:) containing all possible pairs of key species
1629 : ! used in either upper or lower atmos
1630 : !
1631 3072 : subroutine create_flavor(key_species, flavor)
1632 : integer, dimension(:,:,:), intent(in) :: key_species
1633 : integer, dimension(:,:), allocatable, intent(out) :: flavor
1634 6144 : integer, dimension(2,size(key_species,3)*2) :: key_species_list
1635 :
1636 : integer :: ibnd, iatm, i, iflavor
1637 : ! prepare list of key_species
1638 3072 : i = 1
1639 49152 : do ibnd=1,size(key_species,3) ! bands
1640 141312 : do iatm=1,size(key_species,2) ! upper/lower atmosphere
1641 276480 : key_species_list(:,i) = key_species(:,iatm,ibnd)
1642 138240 : i = i + 1
1643 : end do
1644 : end do
1645 : ! rewrite single key_species pairs
1646 95232 : do i=1,size(key_species_list,2)
1647 279552 : key_species_list(:,i) = rewrite_key_species_pair(key_species_list(:,i))
1648 : end do
1649 : ! count unique key species pairs
1650 3072 : iflavor = 0
1651 95232 : do i=1,size(key_species_list,2)
1652 95232 : if (.not.key_species_pair_exists(key_species_list(:,1:i-1),key_species_list(:,i))) then
1653 29184 : iflavor = iflavor + 1
1654 : end if
1655 : end do
1656 : ! fill flavors
1657 9216 : allocate(flavor(2,iflavor))
1658 3072 : iflavor = 0
1659 95232 : do i=1,size(key_species_list,2)
1660 95232 : if (.not.key_species_pair_exists(key_species_list(:,1:i-1),key_species_list(:,i))) then
1661 29184 : iflavor = iflavor + 1
1662 87552 : flavor(:,iflavor) = key_species_list(:,i)
1663 : end if
1664 : end do
1665 3072 : end subroutine create_flavor
1666 : ! ---------------------------------------------------------------------------------------
1667 : !
1668 : ! create index list for extracting col_gas needed for minor gas optical depth calculations
1669 : !
1670 18432 : subroutine create_idx_minor(gas_names, &
1671 6144 : gas_minor, identifier_minor, minor_gases_atm, idx_minor_atm)
1672 : character(len=*), dimension(:), intent(in) :: gas_names
1673 : character(len=*), dimension(:), intent(in) :: &
1674 : gas_minor, &
1675 : identifier_minor
1676 : character(len=*), dimension(:), intent(in) :: minor_gases_atm
1677 : integer, dimension(:), allocatable, &
1678 : intent(out) :: idx_minor_atm
1679 :
1680 : ! local
1681 : integer :: imnr
1682 : integer :: idx_mnr
1683 18432 : allocate(idx_minor_atm(size(minor_gases_atm,dim=1)))
1684 185856 : do imnr = 1, size(minor_gases_atm,dim=1) ! loop over minor absorbers in each band
1685 : ! Find identifying string for minor species in list of possible identifiers (e.g. h2o_slf)
1686 179712 : idx_mnr = string_loc_in_array(minor_gases_atm(imnr), identifier_minor)
1687 : ! Find name of gas associated with minor species identifier (e.g. h2o)
1688 185856 : idx_minor_atm(imnr) = string_loc_in_array(gas_minor(idx_mnr), gas_names)
1689 : enddo
1690 6144 : end subroutine create_idx_minor
1691 :
1692 : ! ---------------------------------------------------------------------------------------
1693 : !
1694 : ! create index for special treatment in density scaling of minor gases
1695 : !
1696 18432 : subroutine create_idx_minor_scaling(gas_names, &
1697 6144 : scaling_gas_atm, idx_minor_scaling_atm)
1698 : character(len=*), dimension(:), intent(in) :: gas_names
1699 : character(len=*), dimension(:), intent(in) :: scaling_gas_atm
1700 : integer, dimension(:), allocatable, &
1701 : intent(out) :: idx_minor_scaling_atm
1702 :
1703 : ! local
1704 : integer :: imnr
1705 18432 : allocate(idx_minor_scaling_atm(size(scaling_gas_atm,dim=1)))
1706 185856 : do imnr = 1, size(scaling_gas_atm,dim=1) ! loop over minor absorbers in each band
1707 : ! This will be -1 if there's no interacting gas
1708 185856 : idx_minor_scaling_atm(imnr) = string_loc_in_array(scaling_gas_atm(imnr), gas_names)
1709 : enddo
1710 6144 : end subroutine create_idx_minor_scaling
1711 : !--------------------------------------------------------------------------------------------------------------------
1712 : ! Is the object ready to use?
1713 : !
1714 1138471 : pure function is_loaded(this)
1715 : class(ty_gas_optics_rrtmgp), intent(in) :: this
1716 : logical(wl) :: is_loaded
1717 :
1718 1138471 : is_loaded = allocated(this%kmajor)
1719 1138471 : end function is_loaded
1720 : !--------------------------------------------------------------------------------------------------------------------
1721 : !
1722 : ! Reset the object to un-initialized state
1723 : !
1724 3072 : subroutine finalize(this)
1725 : class(ty_gas_optics_rrtmgp), intent(inout) :: this
1726 : real(wp), dimension(:), allocatable :: press_ref, press_ref_log, temp_ref
1727 :
1728 3072 : if(this%is_loaded()) then
1729 : !$acc exit data delete(this%gas_names, this%vmr_ref, this%flavor) &
1730 : !$acc delete(this%gpoint_flavor, this%kmajor) &
1731 : !$acc delete(this%minor_limits_gpt_lower) &
1732 : !$acc delete(this%minor_scales_with_density_lower, this%scale_by_complement_lower) &
1733 : !$acc delete(this%idx_minor_lower, this%idx_minor_scaling_lower) &
1734 : !$acc delete(this%kminor_start_lower, this%kminor_lower) &
1735 : !$acc delete(this%minor_limits_gpt_upper) &
1736 : !$acc delete(this%minor_scales_with_density_upper, this%scale_by_complement_upper) &
1737 : !$acc delete(this%idx_minor_upper, this%idx_minor_scaling_upper) &
1738 : !$acc delete(this%kminor_start_upper, this%kminor_upper)
1739 : !$omp target exit data map(release:this%gas_names, this%vmr_ref, this%flavor) &
1740 : !$omp map(release:this%gpoint_flavor, this%kmajor) &
1741 : !$omp map(release:this%minor_limits_gpt_lower) &
1742 : !$omp map(release:this%minor_scales_with_density_lower, this%scale_by_complement_lower) &
1743 : !$omp map(release:this%idx_minor_lower, this%idx_minor_scaling_lower) &
1744 : !$omp map(release:this%kminor_start_lower, this%kminor_lower) &
1745 : !$omp map(release:this%minor_limits_gpt_upper) &
1746 : !$omp map(release:this%minor_scales_with_density_upper, this%scale_by_complement_upper) &
1747 : !$omp map(release:this%idx_minor_upper, this%idx_minor_scaling_upper) &
1748 : !$omp map(release:this%kminor_start_upper, this%kminor_upper)
1749 0 : deallocate(this%gas_names, this%vmr_ref, this%flavor, this%gpoint_flavor, this%kmajor)
1750 0 : deallocate(this%minor_limits_gpt_lower, &
1751 0 : this%minor_scales_with_density_lower, this%scale_by_complement_lower, &
1752 0 : this%idx_minor_lower, this%idx_minor_scaling_lower, this%kminor_start_lower, this%kminor_lower)
1753 0 : deallocate(this%minor_limits_gpt_upper, &
1754 0 : this%minor_scales_with_density_upper, this%scale_by_complement_upper, &
1755 0 : this%idx_minor_upper, this%idx_minor_scaling_upper, this%kminor_start_upper, this%kminor_upper)
1756 :
1757 0 : if(allocated(this%krayl)) then
1758 : !$acc exit data delete(this%krayl)
1759 : !$omp target exit data map(release:this%krayl)
1760 0 : deallocate(this%krayl)
1761 : end if
1762 :
1763 0 : if(allocated(this%planck_frac)) then
1764 : !$acc exit data delete(this%planck_frac, this%totplnk, this%optimal_angle_fit)
1765 : !$omp target exit data map(release:this%planck_frac, this%totplnk, this%optimal_angle_fit)
1766 0 : deallocate(this%planck_frac, this%totplnk, this%optimal_angle_fit)
1767 : end if
1768 :
1769 0 : if(allocated(this%solar_source)) then
1770 : !$acc exit data delete(this%solar_source, this%solar_source_quiet) &
1771 : !$acc delete(this%solar_source_facular,this%solar_source_sunspot)
1772 : !$omp target exit data map(release:this%solar_source, this%solar_source_quiet)
1773 : !$omp map(release:this%solar_source_facular,this%solar_source_sunspot)
1774 : deallocate(this%solar_source, &
1775 0 : this%solar_source_quiet, this%solar_source_facular, this%solar_source_sunspot)
1776 : end if
1777 : !$acc exit data delete(this)
1778 : !$omp target exit data map(release:this)
1779 : end if
1780 :
1781 3072 : end subroutine finalize
1782 : ! ---------------------------------------------------------------------------------------
1783 3072 : subroutine create_key_species_reduce(gas_names,gas_names_red, &
1784 3072 : key_species,key_species_red,key_species_present_init)
1785 : character(len=*), &
1786 : dimension(:), intent(in) :: gas_names
1787 : character(len=*), &
1788 : dimension(:), intent(in) :: gas_names_red
1789 : integer, dimension(:,:,:), intent(in) :: key_species
1790 : integer, dimension(:,:,:), allocatable, intent(out) :: key_species_red
1791 :
1792 : logical, dimension(:), allocatable, intent(out) :: key_species_present_init
1793 : integer :: ip, ia, it, np, na, nt
1794 :
1795 3072 : np = size(key_species,dim=1)
1796 3072 : na = size(key_species,dim=2)
1797 3072 : nt = size(key_species,dim=3)
1798 0 : allocate(key_species_red(size(key_species,dim=1), &
1799 : size(key_species,dim=2), &
1800 15360 : size(key_species,dim=3)))
1801 9216 : allocate(key_species_present_init(size(gas_names)))
1802 61440 : key_species_present_init = .true.
1803 :
1804 9216 : do ip = 1, np
1805 21504 : do ia = 1, na
1806 202752 : do it = 1, nt
1807 196608 : if (key_species(ip,ia,it) .ne. 0) then
1808 122880 : key_species_red(ip,ia,it) = string_loc_in_array(gas_names(key_species(ip,ia,it)),gas_names_red)
1809 122880 : if (key_species_red(ip,ia,it) .eq. -1) key_species_present_init(key_species(ip,ia,it)) = .false.
1810 : else
1811 61440 : key_species_red(ip,ia,it) = key_species(ip,ia,it)
1812 : endif
1813 : enddo
1814 : end do
1815 : enddo
1816 :
1817 3072 : end subroutine create_key_species_reduce
1818 :
1819 : ! ---------------------------------------------------------------------------------------
1820 :
1821 6144 : subroutine reduce_minor_arrays(available_gases, &
1822 6144 : gas_minor,identifier_minor,&
1823 12288 : kminor_atm, &
1824 6144 : minor_gases_atm, &
1825 6144 : minor_limits_gpt_atm, &
1826 12288 : minor_scales_with_density_atm, &
1827 6144 : scaling_gas_atm, &
1828 6144 : scale_by_complement_atm, &
1829 6144 : kminor_start_atm, &
1830 : kminor_atm_red, &
1831 : minor_gases_atm_red, &
1832 : minor_limits_gpt_atm_red, &
1833 : minor_scales_with_density_atm_red, &
1834 : scaling_gas_atm_red, &
1835 : scale_by_complement_atm_red, &
1836 : kminor_start_atm_red)
1837 :
1838 : class(ty_gas_concs), intent(in) :: available_gases
1839 : real(wp), dimension(:,:,:), intent(in) :: kminor_atm
1840 : character(len=*), dimension(:), intent(in) :: gas_minor, &
1841 : identifier_minor
1842 : character(len=*), dimension(:), intent(in) :: minor_gases_atm
1843 : integer, dimension(:,:), intent(in) :: minor_limits_gpt_atm
1844 : logical(wl), dimension(:), intent(in) :: minor_scales_with_density_atm
1845 : character(len=*), dimension(:), intent(in) :: scaling_gas_atm
1846 : logical(wl), dimension(:), intent(in) :: scale_by_complement_atm
1847 : integer, dimension(:), intent(in) :: kminor_start_atm
1848 : real(wp), dimension(:,:,:), allocatable, &
1849 : intent(out) :: kminor_atm_red
1850 : character(len=*), dimension(:), allocatable, &
1851 : intent(out) :: minor_gases_atm_red
1852 : integer, dimension(:,:), allocatable, &
1853 : intent(out) :: minor_limits_gpt_atm_red
1854 : logical(wl), dimension(:), allocatable, &
1855 : intent(out) ::minor_scales_with_density_atm_red
1856 : character(len=*), dimension(:), allocatable, &
1857 : intent(out) ::scaling_gas_atm_red
1858 : logical(wl), dimension(:), allocatable, intent(out) :: &
1859 : scale_by_complement_atm_red
1860 : integer, dimension(:), allocatable, intent(out) :: &
1861 : kminor_start_atm_red
1862 :
1863 : ! Local variables
1864 : integer :: i, j, ks
1865 : integer :: idx_mnr, nm, tot_g, red_nm
1866 : integer :: icnt, n_elim, ng
1867 6144 : logical, dimension(:), allocatable :: gas_is_present
1868 6144 : integer, dimension(:), allocatable :: indexes
1869 6144 : real(wp), dimension(:,:,:), allocatable :: kminor_atm_red_t
1870 :
1871 6144 : nm = size(minor_gases_atm)
1872 6144 : tot_g=0
1873 18432 : allocate(gas_is_present(nm))
1874 239616 : do i = 1, size(minor_gases_atm)
1875 233472 : idx_mnr = string_loc_in_array(minor_gases_atm(i), identifier_minor)
1876 : ! Next line causes a compiler bug in gfortran 11.0.1 on Mac ARM
1877 : ! Should replace gas_names with get_gas_names() and make gas_names private in ty_gas_concs
1878 233472 : gas_is_present(i) = string_in_array(gas_minor(idx_mnr),available_gases%gas_names)
1879 239616 : if(gas_is_present(i)) then
1880 179712 : tot_g = tot_g + (minor_limits_gpt_atm(2,i)-minor_limits_gpt_atm(1,i)+1)
1881 : endif
1882 : enddo
1883 239616 : red_nm = count(gas_is_present)
1884 :
1885 0 : allocate(minor_gases_atm_red (red_nm),&
1886 0 : minor_scales_with_density_atm_red(red_nm), &
1887 0 : scaling_gas_atm_red (red_nm), &
1888 0 : scale_by_complement_atm_red (red_nm), &
1889 55296 : kminor_start_atm_red (red_nm))
1890 18432 : allocate(minor_limits_gpt_atm_red(2, red_nm))
1891 30720 : allocate(kminor_atm_red_t(tot_g, size(kminor_atm,2), size(kminor_atm,3)))
1892 30720 : allocate(kminor_atm_red(size(kminor_atm,3),size(kminor_atm,2),tot_g))
1893 :
1894 6144 : if ((red_nm .eq. nm)) then
1895 : ! Character data not allowed in OpenACC regions?
1896 0 : minor_gases_atm_red = minor_gases_atm
1897 0 : scaling_gas_atm_red = scaling_gas_atm
1898 0 : kminor_atm_red_t = kminor_atm
1899 0 : minor_limits_gpt_atm_red = minor_limits_gpt_atm
1900 0 : minor_scales_with_density_atm_red = minor_scales_with_density_atm
1901 0 : scale_by_complement_atm_red = scale_by_complement_atm
1902 0 : kminor_start_atm_red = kminor_start_atm
1903 : else
1904 12288 : allocate(indexes(red_nm))
1905 : ! Find the integer indexes for the gases that are present
1906 473088 : indexes = pack([(i, i = 1, size(minor_gases_atm))], mask=gas_is_present)
1907 :
1908 192000 : minor_gases_atm_red = minor_gases_atm (indexes)
1909 192000 : scaling_gas_atm_red = scaling_gas_atm (indexes)
1910 : minor_scales_with_density_atm_red = &
1911 192000 : minor_scales_with_density_atm(indexes)
1912 : scale_by_complement_atm_red = &
1913 192000 : scale_by_complement_atm(indexes)
1914 192000 : kminor_start_atm_red = kminor_start_atm (indexes)
1915 :
1916 : icnt = 0
1917 : n_elim = 0
1918 239616 : do i = 1, nm
1919 233472 : ng = minor_limits_gpt_atm(2,i)-minor_limits_gpt_atm(1,i)+1
1920 239616 : if(gas_is_present(i)) then
1921 179712 : icnt = icnt + 1
1922 539136 : minor_limits_gpt_atm_red(1:2,icnt) = minor_limits_gpt_atm(1:2,i)
1923 179712 : kminor_start_atm_red(icnt) = kminor_start_atm(i)-n_elim
1924 179712 : ks = kminor_start_atm_red(icnt)
1925 1715712 : do j = 1, ng
1926 0 : kminor_atm_red_t(kminor_start_atm_red(icnt)+j-1,:,:) = &
1927 216755712 : kminor_atm(kminor_start_atm(i)+j-1,:,:)
1928 : enddo
1929 : else
1930 53760 : n_elim = n_elim + ng
1931 : endif
1932 : enddo
1933 : endif
1934 :
1935 : kminor_atm_red = RESHAPE(kminor_atm_red_t,(/size(kminor_atm_red_t,dim=3), &
1936 24576 : size(kminor_atm_red_t,dim=2),size(kminor_atm_red_t,dim=1)/), ORDER=(/3,2,1/))
1937 6144 : deallocate(kminor_atm_red_t)
1938 12288 : end subroutine reduce_minor_arrays
1939 :
1940 : ! ---------------------------------------------------------------------------------------
1941 : ! returns flavor index; -1 if not found
1942 737280 : pure function key_species_pair2flavor(flavor, key_species_pair)
1943 : integer :: key_species_pair2flavor
1944 : integer, dimension(:,:), intent(in) :: flavor
1945 : integer, dimension(2), intent(in) :: key_species_pair
1946 : integer :: iflav
1947 2632704 : do iflav=1,size(flavor,2)
1948 4690944 : if (all(key_species_pair(:).eq.flavor(:,iflav))) then
1949 737280 : key_species_pair2flavor = iflav
1950 737280 : return
1951 : end if
1952 : end do
1953 0 : key_species_pair2flavor = -1
1954 0 : end function key_species_pair2flavor
1955 :
1956 : ! ---------------------------------------------------------------------------------------
1957 : !
1958 : ! create gpoint_flavor list
1959 : ! a map pointing from each g-point to the corresponding entry in the "flavor list"
1960 : !
1961 3072 : subroutine create_gpoint_flavor(key_species, gpt2band, flavor, gpoint_flavor)
1962 : integer, dimension(:,:,:), intent(in) :: key_species
1963 : integer, dimension(:), intent(in) :: gpt2band
1964 : integer, dimension(:,:), intent(in) :: flavor
1965 : integer, dimension(:,:), intent(out), allocatable :: gpoint_flavor
1966 : integer :: ngpt, igpt, iatm
1967 3072 : ngpt = size(gpt2band)
1968 9216 : allocate(gpoint_flavor(2,ngpt))
1969 371712 : do igpt=1,ngpt
1970 1108992 : do iatm=1,2
1971 0 : gpoint_flavor(iatm,igpt) = key_species_pair2flavor( &
1972 : flavor, &
1973 737280 : rewrite_key_species_pair(key_species(:,iatm,gpt2band(igpt))) &
1974 1843200 : )
1975 : end do
1976 : end do
1977 3072 : end subroutine create_gpoint_flavor
1978 :
1979 : !--------------------------------------------------------------------------------------------------------------------
1980 : !
1981 : ! Utility function to combine optical depths from gas absorption and Rayleigh scattering
1982 : ! It may be more efficient to combine scattering and absorption optical depths in place
1983 : ! rather than storing and processing two large arrays
1984 : !
1985 389263 : subroutine combine_abs_and_rayleigh(tau, tau_rayleigh, optical_props)
1986 : real(wp), dimension(:,:,:), intent(in ) :: tau
1987 : real(wp), dimension(:,:,:), intent(in ) :: tau_rayleigh
1988 : class(ty_optical_props_arry), intent(inout) :: optical_props
1989 :
1990 : integer :: icol, ilay, igpt, ncol, nlay, ngpt, nmom
1991 : real(wp) :: t
1992 :
1993 389263 : ncol = size(tau, 1)
1994 389263 : nlay = size(tau, 2)
1995 389263 : ngpt = size(tau, 3)
1996 : select type(optical_props)
1997 : type is (ty_optical_props_1scl)
1998 : !
1999 : ! Extinction optical depth
2000 : !
2001 : !$acc parallel loop gang vector collapse(3) default(present)
2002 : !$omp target teams distribute parallel do simd collapse(3)
2003 0 : do igpt = 1, ngpt
2004 0 : do ilay = 1, nlay
2005 0 : do icol = 1, ncol
2006 0 : optical_props%tau(icol,ilay,igpt) = tau(icol,ilay,igpt) + &
2007 0 : tau_rayleigh(icol,ilay,igpt)
2008 : end do
2009 : end do
2010 : end do
2011 : !
2012 : ! asymmetry factor or phase function moments
2013 : !
2014 : type is (ty_optical_props_2str)
2015 : !
2016 : ! Extinction optical depth and single scattering albedo
2017 : !
2018 : !$acc parallel loop gang vector collapse(3) default(present)
2019 : !$omp target teams distribute parallel do simd collapse(3)
2020 43986719 : do igpt = 1, ngpt
2021 4142147583 : do ilay = 1, nlay
2022 65796884720 : do icol = 1, ncol
2023 61655126400 : t = tau(icol,ilay,igpt) + tau_rayleigh(icol,ilay,igpt)
2024 61655126400 : if(t > 2._wp * tiny(t)) then
2025 61655126400 : optical_props%ssa(icol,ilay,igpt) = tau_rayleigh(icol,ilay,igpt) / t
2026 : else
2027 0 : optical_props%ssa(icol,ilay,igpt) = 0._wp
2028 : end if
2029 65753287264 : optical_props%tau(icol,ilay,igpt) = t
2030 : end do
2031 : end do
2032 : end do
2033 778526 : call zero_array(ncol, nlay, ngpt, optical_props%g)
2034 : type is (ty_optical_props_nstr)
2035 : !
2036 : ! Extinction optical depth and single scattering albedo
2037 : !
2038 : !$acc parallel loop gang vector collapse(3) default(present)
2039 : !$omp target teams distribute parallel do simd collapse(3)
2040 0 : do igpt = 1, ngpt
2041 0 : do ilay = 1, nlay
2042 0 : do icol = 1, ncol
2043 0 : t = tau(icol,ilay,igpt) + tau_rayleigh(icol,ilay,igpt)
2044 0 : if(t > 2._wp * tiny(t)) then
2045 0 : optical_props%ssa(icol,ilay,igpt) = tau_rayleigh(icol,ilay,igpt) / t
2046 : else
2047 0 : optical_props%ssa(icol,ilay,igpt) = 0._wp
2048 : end if
2049 0 : optical_props%tau(icol,ilay,igpt) = t
2050 : end do
2051 : end do
2052 : end do
2053 0 : nmom = size(optical_props%p, 1)
2054 0 : call zero_array(nmom, ncol, nlay, ngpt, optical_props%p)
2055 0 : if(nmom >= 2) then
2056 : !$acc parallel loop gang vector collapse(3) default(present)
2057 : !$omp target teams distribute parallel do simd collapse(3)
2058 0 : do igpt = 1, ngpt
2059 0 : do ilay = 1, nlay
2060 0 : do icol = 1, ncol
2061 0 : optical_props%p(2,icol,ilay,igpt) = 0.1_wp
2062 : end do
2063 : end do
2064 : end do
2065 : end if
2066 : end select
2067 389263 : end subroutine combine_abs_and_rayleigh
2068 : !--------------------------------------------------------------------------------------------------------------------
2069 : ! Sizes of tables: pressure, temperate, eta (mixing fraction)
2070 : ! Equivalent routines for the number of gases and flavors (get_ngas(), get_nflav()) are defined above because they're
2071 : ! used in function defintions
2072 : ! Table kmajor has dimensions (ngpt, neta, npres, ntemp)
2073 : !--------------------------------------------------------------------------------------------------------------------
2074 : !
2075 : ! return extent of eta dimension
2076 : !
2077 1881535 : pure function get_neta(this)
2078 : class(ty_gas_optics_rrtmgp), intent(in) :: this
2079 : integer :: get_neta
2080 :
2081 1881535 : get_neta = size(this%kmajor,dim=2)
2082 1881535 : end function
2083 : ! --------------------------------------------------------------------------------------
2084 : !
2085 : ! return the number of pressures in reference profile
2086 : ! absorption coefficient table is one bigger since a pressure is repeated in upper/lower atmos
2087 : !
2088 1881535 : pure function get_npres(this)
2089 : class(ty_gas_optics_rrtmgp), intent(in) :: this
2090 : integer :: get_npres
2091 :
2092 1881535 : get_npres = size(this%kmajor,dim=3)-1
2093 1881535 : end function get_npres
2094 : ! --------------------------------------------------------------------------------------
2095 : !
2096 : ! return the number of temperatures
2097 : !
2098 1881535 : pure function get_ntemp(this)
2099 : class(ty_gas_optics_rrtmgp), intent(in) :: this
2100 : integer :: get_ntemp
2101 :
2102 1881535 : get_ntemp = size(this%kmajor,dim=1)
2103 1881535 : end function get_ntemp
2104 : ! --------------------------------------------------------------------------------------
2105 : !
2106 : ! return the number of temperatures for Planck function
2107 : !
2108 746136 : pure function get_nPlanckTemp(this)
2109 : class(ty_gas_optics_rrtmgp), intent(in) :: this
2110 : integer :: get_nPlanckTemp
2111 :
2112 746136 : get_nPlanckTemp = size(this%totplnk,dim=1) ! dimensions are Planck-temperature, band
2113 746136 : end function get_nPlanckTemp
2114 2666205 : end module mo_gas_optics_rrtmgp
|
