Line data Source code
1 : ! This code is part of Radiative Transfer for Energetics (RTE)
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,
8 : ! Trustees of Columbia University in the City of New York
9 : ! All right reserved.
10 : !
11 : ! Use and duplication is permitted under the terms of the
12 : ! BSD 3-clause license, see http://opensource.org/licenses/BSD-3-Clause
13 : ! -------------------------------------------------------------------------------------------------
14 : !
15 : !> Compute longwave radiative fluxes
16 : !>
17 : !> Contains a single routine to compute direct and diffuse fluxes of solar radiation given
18 : !>
19 : !> - atmospheric optical properties, spectrally-resolved via one of the sub-classes of
20 : !> [[mo_optical_props(module):ty_optical_props_arry(type)]] in module [[mo_optical_props]]
21 : ! (ty_optical_props_arry in module mo_optical_props)
22 : !> - information about vertical ordering
23 : !> - internal Planck source functions, defined per g-point on the same spectral grid at the atmosphere,
24 : !> via [[mo_source_functions(module):ty_source_func_lw(type)]] in module [[mo_source_functions]]
25 : ! (ty_source_func_lw in module mo_source_functions)
26 : !> - boundary conditions: surface emissivity defined per band
27 : !> - optionally, a boundary condition for incident diffuse radiation
28 : !> - optionally, an integer number of angles at which to do Gaussian quadrature if scattering is neglected
29 : !>
30 : !> If optical properties are supplied via class ty_optical_props_1scl (absorption optical thickenss only)
31 : !> ([[mo_optical_props(module):ty_optical_props_1scl(type)]] in module [[mo_optical_props]])
32 : !> then an emission/absorption solver is called.
33 : !> If optical properties are supplied via class ty_optical_props_2str
34 : !> ([[mo_optical_props(module):ty_optical_props_2str(type)]] in module [[mo_optical_props]])
35 : !> fluxes are computed via a rescaling by default or, optionally, using two-stream calculations and adding.
36 : !>
37 : !> Users must ensure that emissivity is on the same spectral grid as the optical properties.
38 : !>
39 : !> Final output is via user-extensible ty_fluxes
40 : !> ([[mo_fluxes(module):ty_fluxes(type)]] in module [[mo_fluxes]])
41 : !> which must reduce the detailed spectral fluxes to whatever summary the user needs
42 : !
43 : !> The routine does error checking and choses which lower-level kernel to invoke based on
44 : !> what kinds of optical properties are supplied
45 : !
46 : ! -------------------------------------------------------------------------------------------------
47 : module mo_rte_lw
48 : use mo_rte_kind, only: wp, wl
49 : use mo_rte_config, only: check_extents, check_values
50 : use mo_rte_util_array,only: zero_array
51 : use mo_rte_util_array_validation, &
52 : only: any_vals_less_than, any_vals_outside, extents_are
53 : use mo_optical_props, only: ty_optical_props, &
54 : ty_optical_props_arry, ty_optical_props_1scl, ty_optical_props_2str, ty_optical_props_nstr
55 : use mo_source_functions, &
56 : only: ty_source_func_lw
57 : use mo_fluxes, only: ty_fluxes, ty_fluxes_broadband
58 : use mo_rte_solver_kernels, &
59 : only: lw_solver_noscat, lw_solver_2stream
60 : implicit none
61 : private
62 :
63 : public :: rte_lw
64 : contains
65 : ! --------------------------------------------------
66 : !
67 : ! Interface using only optical properties and source functions as inputs; fluxes as outputs.
68 : !
69 : ! --------------------------------------------------
70 1492272 : function rte_lw(optical_props, top_at_1, &
71 1492272 : sources, sfc_emis, &
72 : fluxes, &
73 1492272 : inc_flux, n_gauss_angles, use_2stream, &
74 2984544 : lw_Ds, flux_up_Jac) result(error_msg)
75 : class(ty_optical_props_arry), intent(in ) :: optical_props
76 : !! Set of optical properties as one or more arrays
77 : logical, intent(in ) :: top_at_1
78 : !! Is the top of the domain at index 1? (if not, ordering is bottom-to-top)
79 : type(ty_source_func_lw), intent(in ) :: sources
80 : !! Derived type with Planck source functions
81 : real(wp), dimension(:,:), intent(in ) :: sfc_emis
82 : !! emissivity at surface [] (nband, ncol)
83 : class(ty_fluxes), intent(inout) :: fluxes
84 : !! Dervied type for computing spectral integrals from g-point fluxes.
85 : !! Default computes broadband fluxes at all levels if output arrays are defined. Can be extended per user desires.
86 : real(wp), dimension(:,:), &
87 : target, optional, intent(in ) :: inc_flux
88 : !! incident flux at domain top [W/m2] (ncol, ngpts)
89 : integer, optional, intent(in ) :: n_gauss_angles
90 : !! Number of angles used in Gaussian quadrature (max 3, no-scattering solution)
91 : logical, optional, intent(in ) :: use_2stream
92 : !! When 2-stream parameters (tau/ssa/g) are provided, use 2-stream methods
93 : !! Default is to use re-scaled longwave transport
94 : real(wp), dimension(:,:), &
95 : optional, intent(in ) :: lw_Ds
96 : !! User-specifed 1/cos of transport angle per col, g-point
97 : real(wp), dimension(:,:), target, &
98 : optional, intent(inout) :: flux_up_Jac
99 : !! surface temperature flux Jacobian [W/m2/K] (ncol, nlay+1)
100 : character(len=128) :: error_msg
101 : !! If empty, calculation was successful
102 : ! --------------------------------
103 : !
104 : ! Local variables
105 : !
106 : integer :: ncol, nlay, ngpt, nband
107 : integer :: icol, ilev, igpt, imu
108 : integer :: n_quad_angs
109 : logical(wl) :: using_2stream, do_Jacobians, do_broadband
110 1492272 : real(wp), dimension(:,:), allocatable :: sfc_emis_gpt
111 1492272 : real(wp), dimension(:,:,:), allocatable :: secants
112 1492272 : real(wp), dimension(:,:), pointer :: jacobian
113 : real(wp), dimension(optical_props%get_ncol(), &
114 : optical_props%get_nlay()+1), target &
115 1492272 : :: decoy2D ! Used for optional outputs - needs to be full size.
116 : ! Memory needs to be allocated for the full g-point fluxes even if they aren't
117 : ! used later because a) the GPU kernels use this memory to work in parallel and
118 : ! b) the fluxes are intent(out) in the solvers
119 : ! Shortwave solver takes a different approach since three fields are needed
120 : real(wp), dimension(optical_props%get_ncol(), &
121 : optical_props%get_nlay()+1, &
122 : optical_props%get_ngpt()) &
123 2984544 : :: gpt_flux_up, gpt_flux_dn
124 1492272 : real(wp), dimension(:,:), pointer :: flux_dn_loc, flux_up_loc
125 1492272 : real(wp), dimension(:,:), pointer :: inc_flux_diffuse
126 : ! --------------------------------------------------
127 : !
128 : ! Weights and angle secants for first order (k=1) Gaussian quadrature.
129 : ! Values from Table 2, Clough et al, 1992, doi:10.1029/92JD01419
130 : ! after Abramowitz & Stegun 1972, page 921
131 : !
132 : integer, parameter :: max_gauss_pts = 4
133 : real(wp), parameter, &
134 : dimension(max_gauss_pts, max_gauss_pts) :: &
135 : gauss_Ds = RESHAPE([1.66_wp, 0._wp, 0._wp, 0._wp, & ! Diffusivity angle, not Gaussian angle
136 : 1.18350343_wp, 2.81649655_wp, 0._wp, 0._wp, &
137 : 1.09719858_wp, 1.69338507_wp, 4.70941630_wp, 0._wp, &
138 : 1.06056257_wp, 1.38282560_wp, 2.40148179_wp, 7.15513024_wp], &
139 : [max_gauss_pts, max_gauss_pts]), &
140 : gauss_wts = RESHAPE([0.5_wp, 0._wp, 0._wp, 0._wp, &
141 : 0.3180413817_wp, 0.1819586183_wp, 0._wp, 0._wp, &
142 : 0.2009319137_wp, 0.2292411064_wp, 0.0698269799_wp, 0._wp, &
143 : 0.1355069134_wp, 0.2034645680_wp, 0.1298475476_wp, 0.0311809710_wp], &
144 : [max_gauss_pts, max_gauss_pts])
145 : ! ------------------------------------------------------------------------------------
146 1492272 : ncol = optical_props%get_ncol()
147 1492272 : nlay = optical_props%get_nlay()
148 1492272 : ngpt = optical_props%get_ngpt()
149 1492272 : nband = optical_props%get_nband()
150 1492272 : do_Jacobians = present(flux_up_Jac)
151 1492272 : error_msg = ""
152 :
153 : ! ------------------------------------------------------------------------------------
154 : !
155 : ! Error checking -- input consistency of sizes and validity of values
156 :
157 1492272 : if(.not. fluxes%are_desired()) &
158 0 : error_msg = "rte_lw: no space allocated for fluxes"
159 :
160 1492272 : if (do_Jacobians .and. check_extents) then
161 0 : if( .not. extents_are(flux_up_Jac, ncol, nlay+1)) &
162 0 : error_msg = "rte_lw: flux Jacobian inconsistently sized"
163 : endif
164 :
165 1492272 : if (check_extents) then
166 : !
167 : ! Source functions
168 : !
169 5969088 : if(any([sources%get_ncol(), sources%get_nlay(), sources%get_ngpt()] /= [ncol, nlay, ngpt])) &
170 0 : error_msg = "rte_lw: sources and optical properties inconsistently sized"
171 : !
172 : ! Surface emissivity
173 : !
174 1492272 : if(.not. extents_are(sfc_emis, nband, ncol)) &
175 0 : error_msg = "rte_lw: sfc_emis inconsistently sized"
176 : !
177 : ! Incident flux, if present
178 : !
179 1492272 : if(present(inc_flux)) then
180 0 : if(.not. extents_are(inc_flux, ncol, ngpt)) &
181 0 : error_msg = "rte_lw: inc_flux inconsistently sized"
182 : end if
183 1492272 : if (present(lw_Ds)) then
184 0 : if(.not. extents_are(lw_Ds, ncol, ngpt)) &
185 0 : error_msg = "rte_lw: lw_Ds inconsistently sized"
186 : end if
187 : end if
188 :
189 1492272 : if(check_values) then
190 1492272 : if(any_vals_outside(sfc_emis, 0._wp, 1._wp)) &
191 0 : error_msg = "rte_lw: sfc_emis has values < 0 or > 1"
192 :
193 1492272 : if(present(inc_flux)) then
194 0 : if(any_vals_less_than(inc_flux, 0._wp)) &
195 0 : error_msg = "rte_lw: inc_flux has values < 0"
196 : end if
197 :
198 1492272 : if (present(lw_Ds)) then
199 0 : if(any_vals_less_than(lw_Ds, 1._wp)) &
200 0 : error_msg = "rte_lw: one or more values of lw_Ds < 1."
201 : end if
202 :
203 1492272 : if(present(n_gauss_angles)) then
204 0 : if(n_gauss_angles > max_gauss_pts) &
205 0 : error_msg = "rte_lw: asking for too many quadrature points for no-scattering calculation"
206 0 : if(n_gauss_angles < 1) &
207 0 : error_msg = "rte_lw: have to ask for at least one quadrature point for no-scattering calculation"
208 : end if
209 : end if
210 1492272 : if(len_trim(error_msg) > 0) return
211 :
212 : !
213 : ! Number of quadrature points for no-scattering calculation
214 : !
215 1492272 : n_quad_angs = 1
216 1492272 : if(present(n_gauss_angles)) n_quad_angs = n_gauss_angles
217 : !
218 : ! Optionally - use 2-stream methods when low-order scattering properties are provided?
219 : !
220 1492272 : using_2stream = .false.
221 1492272 : if(present(use_2stream)) using_2stream = use_2stream
222 :
223 : !
224 : ! Checking that optional arguments are consistent with one another and with optical properties
225 : !
226 : select type (optical_props)
227 : class is (ty_optical_props_1scl)
228 1492272 : if (using_2stream) &
229 0 : error_msg = "rte_lw: can't use two-stream methods with only absorption optical depth"
230 1492272 : if(present(lw_Ds) .and. n_quad_angs /= 1) &
231 0 : error_msg = "rte_lw: providing lw_Ds incompatible with specifying n_gauss_angles"
232 : class is (ty_optical_props_2str)
233 0 : if (present(lw_Ds)) &
234 0 : error_msg = "rte_lw: lw_Ds not valid when providing scattering optical properties"
235 0 : if (using_2stream .and. n_quad_angs /= 1) &
236 0 : error_msg = "rte_lw: using_2stream=true incompatible with specifying n_gauss_angles"
237 0 : if (using_2stream .and. do_Jacobians) &
238 0 : error_msg = "rte_lw: can't provide Jacobian of fluxes w.r.t surface temperature with 2-stream"
239 : class default
240 0 : error_msg = "rte_lw: lw_solver(...ty_optical_props_nstr...) not yet implemented"
241 : end select
242 :
243 1492272 : if(len_trim(error_msg) > 0) then
244 0 : if(len_trim(optical_props%get_name()) > 0) &
245 0 : error_msg = trim(optical_props%get_name()) // ': ' // trim(error_msg)
246 0 : return
247 : end if
248 :
249 : ! ------------------------------------------------------------------------------------
250 : ! Boundary conditions
251 : ! Lower boundary condition -- expand surface emissivity by band to gpoints
252 : !
253 5969088 : allocate(sfc_emis_gpt(ncol, ngpt))
254 :
255 : ! Upper boundary condition - use values in optional arg or be set to 0
256 : !
257 1492272 : if (present(inc_flux)) then
258 0 : inc_flux_diffuse => inc_flux
259 : !$acc enter data copyin( inc_flux_diffuse)
260 : !$omp target enter data map(to: inc_flux_diffuse)
261 : else
262 4476816 : allocate(inc_flux_diffuse(ncol, ngpt))
263 : !$acc enter data create( inc_flux_diffuse)
264 : !$omp target enter data map(alloc:inc_flux_diffuse)
265 1492272 : call zero_array(ncol, ngpt, inc_flux_diffuse)
266 : end if
267 :
268 : ! ------------------------------------------------------------------------------------
269 1492272 : if(do_Jacobians) then
270 0 : jacobian => flux_up_Jac
271 : else
272 1492272 : jacobian => decoy2D
273 : end if
274 :
275 : select type(fluxes)
276 : !
277 : ! Broadband fluxes are treated as a special case within the solvers; memory
278 : ! for both up and down fluxes needs to be available even if the user doesn't
279 : ! want one of them
280 : !
281 : type is (ty_fluxes_broadband)
282 746136 : do_broadband = .true._wl
283 : !
284 : ! Broadband fluxes class has three possible outputs; allocate memory for local use
285 : ! if one or more haven't been requested
286 : !
287 746136 : if(associated(fluxes%flux_up)) then
288 746136 : flux_up_loc => fluxes%flux_up
289 : else
290 0 : allocate(flux_up_loc(ncol, nlay+1))
291 : end if
292 1492272 : if(associated(fluxes%flux_dn)) then
293 746136 : flux_dn_loc => fluxes%flux_dn
294 : else
295 0 : allocate(flux_dn_loc(ncol, nlay+1))
296 : end if
297 : !$acc enter data create( flux_up_loc, flux_dn_loc)
298 : !$omp target enter data map(alloc:flux_up_loc, flux_dn_loc)
299 : class default
300 : !
301 : ! If broadband integrals aren't being computed, allocate working space
302 : ! and decoy addresses for spectrally-integrated fields
303 : !
304 746136 : do_broadband = .false._wl
305 746136 : flux_up_loc => decoy2D
306 746136 : flux_dn_loc => decoy2D
307 : end select
308 :
309 : !
310 : ! Compute the radiative transfer...
311 : !
312 : !$acc data create( sfc_emis_gpt, flux_up_loc, flux_dn_loc, gpt_flux_up, gpt_flux_dn)
313 : !$omp target data map(alloc:sfc_emis_gpt, flux_up_loc, flux_dn_loc, gpt_flux_up, gpt_flux_dn)
314 1492272 : call expand_and_transpose(optical_props, sfc_emis, sfc_emis_gpt)
315 1492272 : if(check_values) error_msg = optical_props%validate()
316 1492272 : if(len_trim(error_msg) == 0) then ! Can't do an early return within OpenACC/MP data regions
317 : select type (optical_props)
318 : class is (ty_optical_props_1scl)
319 : !
320 : ! No scattering two-stream calculation
321 : !
322 : !
323 : ! Secant of radiation angle - either user-supplied, one per g-point, or
324 : ! taken from first-order Gaussian quadrate and applied to all columns a g-points
325 : !
326 7461360 : allocate(secants(ncol, ngpt, n_quad_angs))
327 : !$acc data create( secants)
328 : !$omp target data map(alloc:secants)
329 1492272 : if (present(lw_Ds)) then
330 : !$acc parallel loop collapse(2) copyin(lw_Ds)
331 : !$omp target teams distribute parallel do simd collapse(2)
332 : ! nmu is 1
333 0 : do igpt = 1, ngpt
334 0 : do icol = 1, ncol
335 0 : secants(icol,igpt,1) = lw_Ds(icol,igpt)
336 : end do
337 : end do
338 : else
339 : !
340 : ! Is there an alternative to making ncol x ngpt copies of each value?
341 : !
342 : !$acc parallel loop collapse(3)
343 : !$omp target teams distribute parallel do simd collapse(3)
344 2984544 : do imu = 1, n_quad_angs
345 193995360 : do igpt = 1, ngpt
346 3190928688 : do icol = 1, ncol
347 3189436416 : secants(icol,igpt,imu) = gauss_Ds(imu,n_quad_angs)
348 : end do
349 : end do
350 : end do
351 : end if
352 : call lw_solver_noscat(ncol, nlay, ngpt, &
353 : logical(top_at_1, wl), n_quad_angs, &
354 0 : secants, gauss_wts(1:n_quad_angs,n_quad_angs), &
355 : optical_props%tau, &
356 : sources%lay_source, sources%lev_source_inc, &
357 : sources%lev_source_dec, &
358 : sfc_emis_gpt, sources%sfc_source, &
359 : inc_flux_diffuse, &
360 : gpt_flux_up, gpt_flux_dn, &
361 : do_broadband, flux_up_loc, flux_dn_loc, &
362 : logical(do_Jacobians, wl), sources%sfc_source_Jac, jacobian, &
363 1492272 : logical(.false., wl), optical_props%tau, optical_props%tau)
364 : ! The last two arguments won't be used since the
365 : ! third-to-last is .false. but need valid addresses
366 : !$acc end data
367 : !$omp end target data
368 : class is (ty_optical_props_2str)
369 0 : if (using_2stream) then
370 : !
371 : ! two-stream calculation with scattering
372 : !
373 : call lw_solver_2stream(ncol, nlay, ngpt, logical(top_at_1, wl), &
374 : optical_props%tau, optical_props%ssa, optical_props%g, &
375 : sources%lay_source, sources%lev_source_inc, sources%lev_source_dec, &
376 : sfc_emis_gpt, sources%sfc_source, &
377 : inc_flux_diffuse, &
378 0 : gpt_flux_up, gpt_flux_dn)
379 : else
380 0 : allocate(secants(ncol, ngpt, n_quad_angs))
381 : !$acc data create( secants)
382 : !$omp target data map(alloc:secants)
383 : !$acc parallel loop collapse(3)
384 : !$omp target teams distribute parallel do simd collapse(3)
385 0 : do imu = 1, n_quad_angs
386 0 : do igpt = 1, ngpt
387 0 : do icol = 1, ncol
388 0 : secants(icol,igpt,imu) = gauss_Ds(imu,n_quad_angs)
389 : end do
390 : end do
391 : end do
392 : !
393 : ! Re-scaled solution to account for scattering
394 : !
395 : call lw_solver_noscat(ncol, nlay, ngpt, &
396 : logical(top_at_1, wl), n_quad_angs, &
397 0 : secants, gauss_wts(1:n_quad_angs,n_quad_angs), &
398 : optical_props%tau, &
399 : sources%lay_source, sources%lev_source_inc, &
400 : sources%lev_source_dec, &
401 : sfc_emis_gpt, sources%sfc_source, &
402 : inc_flux_diffuse, &
403 : gpt_flux_up, gpt_flux_dn, &
404 : do_broadband, flux_up_loc, flux_dn_loc, &
405 : logical(do_Jacobians, wl), sources%sfc_source_Jac, jacobian, &
406 0 : logical(.true., wl), optical_props%ssa, optical_props%g)
407 : !$acc end data
408 : !$omp end target data
409 : endif
410 : class is (ty_optical_props_nstr)
411 : !
412 : ! n-stream calculation
413 : !
414 0 : error_msg = 'lw_solver(...ty_optical_props_nstr...) not yet implemented'
415 : end select
416 :
417 : select type(fluxes)
418 : !
419 : ! Tidy up memory for broadband fluxes on GPUs
420 : !
421 : type is (ty_fluxes_broadband)
422 746136 : if(associated(fluxes%flux_net)) then
423 : !
424 : ! FIXME: Do we need the create/copyout here?
425 : !
426 : !$acc parallel loop collapse(2) copyin(fluxes) copyout( fluxes%flux_net)
427 : !$omp target teams distribute parallel do simd collapse(2) map(from:fluxes%flux_net)
428 71629056 : do ilev = 1, nlay+1
429 1184326056 : do icol = 1, ncol
430 1183579920 : fluxes%flux_net(icol,ilev) = flux_dn_loc(icol,ilev) - flux_up_loc(icol,ilev)
431 : end do
432 : end do
433 : end if
434 : class default
435 : !
436 : ! ...or reduce spectral fluxes to desired output quantities
437 : !
438 746136 : error_msg = fluxes%reduce(gpt_flux_up, gpt_flux_dn, optical_props, top_at_1)
439 : end select
440 : end if ! no error message from validation
441 : !$acc end data
442 : !$omp end target data
443 :
444 1492272 : if(.not. present(inc_flux)) then
445 : !$acc exit data delete( inc_flux_diffuse)
446 : !$omp target exit data map(release:inc_flux_diffuse)
447 1492272 : deallocate(inc_flux_diffuse)
448 : end if
449 : select type(fluxes)
450 : type is (ty_fluxes_broadband)
451 : !$acc exit data copyout( flux_up_loc, flux_dn_loc)
452 : !$omp target exit data map(from:flux_up_loc, flux_dn_loc)
453 0 : if(.not. associated(flux_up_loc, fluxes%flux_up)) deallocate(flux_up_loc)
454 746136 : if(.not. associated(flux_dn_loc, fluxes%flux_dn)) deallocate(flux_dn_loc)
455 : end select
456 1492272 : end function rte_lw
457 : !--------------------------------------------------------------------------------------------------------------------
458 : !
459 : ! Expand from band to g-point dimension, transpose dimensions (nband, ncol) -> (ncol,ngpt)
460 : !
461 1492272 : subroutine expand_and_transpose(ops,arr_in,arr_out)
462 : class(ty_optical_props), intent(in ) :: ops
463 : real(wp), dimension(:,:), intent(in ) :: arr_in ! (nband, ncol)
464 : real(wp), dimension(:,:), intent(out) :: arr_out ! (ncol, igpt)
465 : ! -------------
466 : integer :: ncol, nband, ngpt
467 : integer :: icol, iband, igpt
468 1492272 : integer, dimension(2,ops%get_nband()) :: limits
469 :
470 1492272 : ncol = size(arr_in, 2)
471 1492272 : nband = ops%get_nband()
472 1492272 : ngpt = ops%get_ngpt()
473 73121328 : limits = ops%get_band_lims_gpoint()
474 : !$acc parallel loop collapse(2) copyin(arr_in, limits)
475 : !$omp target teams distribute parallel do simd collapse(2) map(to:arr_in, limits)
476 25368624 : do iband = 1, nband
477 400171824 : do icol = 1, ncol
478 3397105152 : do igpt = limits(1, iband), limits(2, iband)
479 3373228800 : arr_out(icol, igpt) = arr_in(iband,icol)
480 : end do
481 : end do
482 : end do
483 :
484 1492272 : end subroutine expand_and_transpose
485 : !--------------------------------------------------------------------------------------------------------------------
486 7461360 : end module mo_rte_lw
|
