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, 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 : !> Compute shortwave radiative fluxes
14 :
15 : !> Contains a single routine to compute direct and diffuse fluxes of solar radiation given
16 : !>
17 : !> - atmospheric optical properties on a spectral grid
18 : !> - information about vertical ordering
19 : !> - boundary conditions
20 : !> - solar zenith angle, spectrally-resolved incident colimated flux, surface albedos for direct and diffuse radiation
21 : !> - optionally, a boundary condition for incident diffuse radiation
22 : !>
23 : !> It is the user's responsibility to ensure that boundary conditions (incident fluxes, surface albedos) are on the same
24 : !> spectral grid as the optical properties.
25 : !>
26 : !> Final output is via user-extensible ty_fluxes
27 : !> ([[mo_fluxes(module):ty_fluxes(type)]] in module [[mo_fluxes]])
28 : !> which must reduce the detailed spectral fluxes to whatever summary the user needs
29 : !>
30 : !>
31 : !> The routine does error checking and choses which lower-level kernel to invoke based on
32 : !> what kinds of optical properties are supplied
33 : !
34 : ! -------------------------------------------------------------------------------------------------
35 : module mo_rte_sw
36 : use mo_rte_kind, only: wp, wl
37 : use mo_rte_config, only: check_extents, check_values
38 : use mo_rte_util_array,only: zero_array
39 : use mo_rte_util_array_validation, &
40 : only: any_vals_less_than, any_vals_outside, extents_are
41 : use mo_optical_props, only: ty_optical_props, &
42 : ty_optical_props_arry, ty_optical_props_1scl, ty_optical_props_2str, ty_optical_props_nstr
43 : use mo_fluxes, only: ty_fluxes, ty_fluxes_broadband
44 : use mo_rte_solver_kernels, &
45 : only: sw_solver_noscat, sw_solver_2stream
46 : implicit none
47 : private
48 :
49 : interface rte_sw
50 : module procedure rte_sw_mu0_bycol, rte_sw_mu0_full
51 : end interface rte_sw
52 : public :: rte_sw
53 :
54 : contains
55 : ! -------------------------------------------------------------------------------------------------
56 778526 : function rte_sw_mu0_bycol(atmos, top_at_1, &
57 778526 : mu0, inc_flux, &
58 778526 : sfc_alb_dir, sfc_alb_dif, &
59 778526 : fluxes, inc_flux_dif) result(error_msg)
60 : class(ty_optical_props_arry), intent(in ) :: atmos
61 : !! Optical properties provided as arrays
62 : logical, intent(in ) :: top_at_1
63 : !! Is the top of the domain at index 1? (if not, ordering is bottom-to-top)
64 : real(wp), dimension(:), intent(in ) :: mu0
65 : !! cosine of solar zenith angle (ncol) - will be assumed constant with height
66 : real(wp), dimension(:,:), intent(in ) :: inc_flux
67 : !! incident flux at top of domain [W/m2] (ncol, ngpt)
68 : real(wp), dimension(:,:), intent(in ) :: sfc_alb_dir
69 : !! surface albedo for direct and
70 : real(wp), dimension(:,:), intent(in ) :: sfc_alb_dif
71 : !! diffuse radiation (nband, ncol)
72 : class(ty_fluxes), intent(inout) :: fluxes
73 : !! Class describing output calculations
74 : real(wp), dimension(:,:), optional, target, &
75 : intent(in ) :: inc_flux_dif
76 : !! incident diffuse flux at top of domain [W/m2] (ncol, ngpt)
77 : character(len=128) :: error_msg
78 : !! If empty, calculation was successful
79 : ! --------------------------------
80 778526 : real(wp), dimension(size(mu0), atmos%get_nlay()) :: mu0_bylay
81 : integer :: i, j, ncol, nlay
82 :
83 778526 : ncol = size(mu0)
84 778526 : nlay = atmos%get_nlay()
85 : ! Solar zenith angle cosine is constant with height
86 : !$acc data copyin(mu0) create(mu0_bylay)
87 : !$omp target data map(to:mu0) map(alloc:mu0_bylay)
88 :
89 : !$acc parallel loop collapse(2)
90 : !$omp target teams distribute parallel do simd collapse(2)
91 73959970 : do j = 1, nlay
92 1174944370 : do i = 1, ncol
93 1174165844 : mu0_bylay(i,j) = mu0(i)
94 : end do
95 : end do
96 :
97 : error_msg = rte_sw_mu0_full(atmos, top_at_1, &
98 : mu0_bylay, inc_flux, &
99 : sfc_alb_dir, sfc_alb_dif, &
100 1557052 : fluxes, inc_flux_dif)
101 : !$acc end data
102 : !$omp end target data
103 778526 : end function rte_sw_mu0_bycol
104 : ! -------------------------------------------------------------------------------------------------
105 778526 : function rte_sw_mu0_full(atmos, top_at_1, &
106 778526 : mu0, inc_flux, &
107 778526 : sfc_alb_dir, sfc_alb_dif, &
108 778526 : fluxes, inc_flux_dif) result(error_msg)
109 : class(ty_optical_props_arry), intent(in ) :: atmos
110 : !! Optical properties provided as arrays
111 : logical, intent(in ) :: top_at_1
112 : !! Is the top of the domain at index 1? (if not, ordering is bottom-to-top)
113 : real(wp), dimension(:,:), intent(in ) :: mu0
114 : !! cosine of solar zenith angle (ncol, nlay)
115 : real(wp), dimension(:,:), intent(in ) :: inc_flux
116 : !! incident flux at top of domain [W/m2] (ncol, ngpt)
117 : real(wp), dimension(:,:), intent(in ) :: sfc_alb_dir
118 : !! surface albedo for direct and
119 : real(wp), dimension(:,:), intent(in ) :: sfc_alb_dif
120 : !! diffuse radiation (nband, ncol)
121 : class(ty_fluxes), intent(inout) :: fluxes
122 : !! Class describing output calculations
123 : real(wp), dimension(:,:), optional, target, &
124 : intent(in ) :: inc_flux_dif
125 : !! incident diffuse flux at top of domain [W/m2] (ncol, ngpt)
126 : character(len=128) :: error_msg
127 : !! If empty, calculation was successful
128 : ! --------------------------------
129 : !
130 : ! Local variables
131 : !
132 : integer :: ncol, nlay, ngpt, nband
133 : integer :: icol, ilev
134 : logical(wl) :: has_dif_bc, do_broadband
135 :
136 778526 : real(wp), dimension(:,:,:), pointer :: gpt_flux_up, gpt_flux_dn, gpt_flux_dir
137 778526 : real(wp), dimension(:,:), allocatable :: sfc_alb_dir_gpt, sfc_alb_dif_gpt
138 778526 : real(wp), dimension(:,:), pointer :: flux_dn_loc, flux_up_loc, flux_dir_loc
139 778526 : real(wp), dimension(:,:), pointer :: inc_flux_diffuse
140 778526 : real(wp), dimension(:,:,:), allocatable, target :: decoy3D
141 778526 : real(wp), dimension(:,:), allocatable, target :: decoy2D
142 :
143 : ! ------------------------------------------------------------------------------------
144 778526 : ncol = atmos%get_ncol()
145 778526 : nlay = atmos%get_nlay()
146 778526 : ngpt = atmos%get_ngpt()
147 778526 : nband = atmos%get_nband()
148 778526 : error_msg = ""
149 : ! ------------------------------------------------------------------------------------
150 : !
151 : ! Error checking -- consistency of sizes and validity of values
152 : !
153 : ! --------------------------------
154 778526 : if(.not. fluxes%are_desired()) &
155 0 : error_msg = "rte_sw: no space allocated for fluxes"
156 :
157 778526 : has_dif_bc = logical(present(inc_flux_dif), wl)
158 :
159 : !
160 : ! Sizes of input arrays
161 : !
162 : ! Copy variables whose sizes and values are checked to the GPU so the checks can happen there.
163 : ! No harm done if checks are not performed (?)
164 : !$acc data copyin(mu0, inc_flux, sfc_alb_dir, sfc_alb_dif)
165 : !$omp target data map(to:mu0, inc_flux, sfc_alb_dir, sfc_alb_dif)
166 : !$acc data copyin(inc_flux_dif) if (has_dif_bc)
167 : !$omp target data map(to:inc_flux_dif) if (has_dif_bc)
168 778526 : if(check_extents) then
169 778526 : if(.not. extents_are(mu0, ncol, nlay)) &
170 0 : error_msg = "rte_sw: mu0 inconsistently sized"
171 778526 : if(.not. extents_are(inc_flux, ncol, ngpt)) &
172 0 : error_msg = "rte_sw: inc_flux inconsistently sized"
173 778526 : if(.not. extents_are(sfc_alb_dir, nband, ncol)) &
174 0 : error_msg = "rte_sw: sfc_alb_dir inconsistently sized"
175 778526 : if(.not. extents_are(sfc_alb_dif, nband, ncol)) &
176 0 : error_msg = "rte_sw: sfc_alb_dif inconsistently sized"
177 778526 : if(has_dif_bc) then
178 0 : if(.not. extents_are(inc_flux_dif, ncol, ngpt)) &
179 0 : error_msg = "rte_sw: inc_flux_dif inconsistently sized"
180 : end if
181 : end if
182 : !
183 : ! Values of input arrays
184 : !
185 778526 : if(check_values) then
186 778526 : if(any_vals_outside(mu0, -1._wp, 1._wp)) &
187 0 : error_msg = "rte_sw: one or more mu0 < -1 or > 1"
188 778526 : if(any_vals_less_than(inc_flux, 0._wp)) &
189 0 : error_msg = "rte_sw: one or more inc_flux < 0"
190 778526 : if(any_vals_outside(sfc_alb_dir, 0._wp, 1._wp)) &
191 0 : error_msg = "rte_sw: sfc_alb_dir out of bounds [0,1]"
192 778526 : if(any_vals_outside(sfc_alb_dif, 0._wp, 1._wp)) &
193 0 : error_msg = "rte_sw: sfc_alb_dif out of bounds [0,1]"
194 778526 : if(has_dif_bc) then
195 0 : if(any_vals_less_than(inc_flux_dif, 0._wp)) &
196 0 : error_msg = "rte_sw: one or more inc_flux_dif < 0"
197 : end if
198 : end if
199 :
200 : ! ------------------------------------------------------------------------------------
201 : select type(fluxes)
202 : type is (ty_fluxes_broadband)
203 389263 : do_broadband = .true._wl
204 : !
205 : ! Solvers will integrate in place (one g-point at a time on CPUs)
206 : ! so won't need big working arrays
207 : !
208 1946315 : allocate(decoy3D(ncol, nlay+1, ngpt))
209 389263 : gpt_flux_up => decoy3D
210 389263 : gpt_flux_dn => decoy3D
211 389263 : gpt_flux_dir => decoy3D
212 : !
213 : ! Broadband fluxes class has three possible outputs; allocate memory for local use
214 : ! if one or more haven't been requested
215 : !
216 389263 : if(associated(fluxes%flux_up)) then
217 389263 : flux_up_loc => fluxes%flux_up
218 : else
219 0 : allocate(flux_up_loc(ncol, nlay+1))
220 : end if
221 389263 : if(associated(fluxes%flux_dn)) then
222 389263 : flux_dn_loc => fluxes%flux_dn
223 : else
224 0 : allocate(flux_dn_loc(ncol, nlay+1))
225 : end if
226 778526 : if(associated(fluxes%flux_dn_dir)) then
227 389263 : flux_dir_loc => fluxes%flux_dn_dir
228 : else
229 0 : allocate(flux_dir_loc(ncol, nlay+1))
230 : end if
231 : !$acc enter data create( flux_up_loc, flux_dn_loc, flux_dir_loc)
232 : !$omp target enter data map(alloc:flux_up_loc, flux_dn_loc, flux_dir_loc)
233 : class default
234 : !
235 : ! If broadband integrals aren't being computed, allocate working space
236 : ! and decoy addresses for spectrally-integrated fields
237 : !
238 389263 : do_broadband = .false._wl
239 1557052 : allocate(decoy2D(ncol, nlay+1))
240 389263 : flux_up_loc => decoy2D
241 389263 : flux_dn_loc => decoy2D
242 389263 : flux_dir_loc => decoy2D
243 : allocate(gpt_flux_up (ncol,nlay+1,ngpt), &
244 : gpt_flux_dn (ncol,nlay+1,ngpt), &
245 4281893 : gpt_flux_dir(ncol,nlay+1,ngpt))
246 : end select
247 :
248 4671156 : allocate(sfc_alb_dir_gpt(ncol, ngpt), sfc_alb_dif_gpt(ncol, ngpt))
249 778526 : if(len_trim(error_msg) > 0) then
250 0 : if(len_trim(atmos%get_name()) > 0) &
251 0 : error_msg = trim(atmos%get_name()) // ': ' // trim(error_msg)
252 : end if
253 :
254 : ! Fluxes need to be copied out only if do_broadband is .true.
255 : !$acc data copyin( flux_up_loc,flux_dn_loc,flux_dir_loc) if ( do_broadband)
256 : !$omp target data map(to: flux_up_loc,flux_dn_loc,flux_dir_loc) if ( do_broadband)
257 : !$acc data create( flux_up_loc,flux_dn_loc,flux_dir_loc) if (.not. do_broadband)
258 : !$omp target data map(alloc:flux_up_loc,flux_dn_loc,flux_dir_loc) if (.not. do_broadband)
259 :
260 : !$acc data create( gpt_flux_up,gpt_flux_dn,gpt_flux_dir) &
261 : !$acc create( sfc_alb_dir_gpt, sfc_alb_dif_gpt)
262 : !$omp target data map(alloc:gpt_flux_up,gpt_flux_dn,gpt_flux_dir) &
263 : !$omp map(alloc:sfc_alb_dir_gpt, sfc_alb_dif_gpt)
264 :
265 :
266 : ! ------------------------------------------------------------------------------------
267 : ! Boundary conditions
268 : ! Lower boundary condition -- expand surface albedos by band to gpoints
269 : ! and switch dimension ordering
270 778526 : call expand_and_transpose(atmos, sfc_alb_dir, sfc_alb_dir_gpt)
271 778526 : call expand_and_transpose(atmos, sfc_alb_dif, sfc_alb_dif_gpt)
272 : !
273 : ! Diffuse flux boundary condition - will use values in optional arg or be set to 0
274 : !
275 778526 : if (has_dif_bc) then
276 0 : inc_flux_diffuse => inc_flux_dif
277 : !$acc enter data copyin( inc_flux_diffuse)
278 : !$omp target enter data map(to: inc_flux_diffuse)
279 : else
280 3114104 : allocate(inc_flux_diffuse(ncol, ngpt))
281 : !$acc enter data create( inc_flux_diffuse)
282 : !$omp target enter data map(alloc:inc_flux_diffuse)
283 778526 : call zero_array(ncol, ngpt, inc_flux_diffuse)
284 : end if
285 : ! ------------------------------------------------------------------------------------
286 778526 : if(check_values) error_msg = atmos%validate()
287 : !
288 : ! Compute the radiative transfer...
289 : !
290 778526 : if(len_trim(error_msg) == 0) then
291 : select type (atmos)
292 : class is (ty_optical_props_1scl)
293 : !
294 : ! Direct beam only - for completeness, unlikely to be used in practice
295 : !
296 : call sw_solver_noscat(ncol, nlay, ngpt, logical(top_at_1, wl), &
297 : atmos%tau, mu0, inc_flux, &
298 0 : gpt_flux_dir)
299 0 : call zero_array(ncol, nlay+1, ngpt, gpt_flux_up)
300 : !
301 : !$acc kernels
302 : !$omp target
303 0 : gpt_flux_dn(:,:,:) = gpt_flux_dir(:,:,:)
304 : !$acc end kernels
305 : !$omp end target
306 :
307 : class is (ty_optical_props_2str)
308 : !
309 : ! two-stream calculation with scattering
310 : !
311 : call sw_solver_2stream(ncol, nlay, ngpt, logical(top_at_1, wl), &
312 : atmos%tau, atmos%ssa, atmos%g, mu0, &
313 : sfc_alb_dir_gpt, sfc_alb_dif_gpt, &
314 : inc_flux, &
315 : gpt_flux_up, gpt_flux_dn, gpt_flux_dir, &
316 : has_dif_bc, inc_flux_diffuse, &
317 778526 : do_broadband, flux_up_loc, flux_dn_loc, flux_dir_loc)
318 : class is (ty_optical_props_nstr)
319 : !
320 : ! n-stream calculation
321 : !
322 : ! not yet implemented so fail
323 : !
324 0 : error_msg = 'sw_solver(...ty_optical_props_nstr...) not yet implemented'
325 : end select
326 778526 : if(len_trim(error_msg) > 0) then
327 0 : if(len_trim(atmos%get_name()) > 0) &
328 0 : error_msg = trim(atmos%get_name()) // ': ' // trim(error_msg)
329 : end if
330 : !
331 : ! Flux reduction (summarizing for output)
332 : !
333 : select type(fluxes)
334 : !
335 : ! Tidy up memory for broadband fluxes
336 : !
337 : type is (ty_fluxes_broadband)
338 389263 : if(associated(fluxes%flux_net)) then
339 : !$acc parallel loop collapse(2) copyout(fluxes%flux_net)
340 : !$omp target teams distribute parallel do simd collapse(2)
341 37369248 : do ilev = 1, nlay+1
342 593717748 : do icol = 1, ncol
343 593328485 : fluxes%flux_net(icol,ilev) = flux_dn_loc(icol,ilev) - flux_up_loc(icol,ilev)
344 : end do
345 : end do
346 : end if
347 : class default
348 : !
349 : ! ...or reduce spectral fluxes to desired output quantities
350 : !
351 389263 : error_msg = fluxes%reduce(gpt_flux_up, gpt_flux_dn, atmos, top_at_1, gpt_flux_dir)
352 : end select
353 : end if ! In case of an error we exit here
354 :
355 : !$acc end data
356 : !$omp end target data
357 : !$acc end data
358 : !$omp end target data
359 : !$acc end data
360 : !$omp end target data
361 : !$acc end data
362 : !$omp end target data
363 : !$acc end data
364 : !$omp end target data
365 :
366 : !
367 : ! Deallocate any memory allocated locally to pointer variables
368 : !
369 : select type(fluxes)
370 : type is (ty_fluxes_broadband)
371 : !$acc exit data copyout( flux_up_loc, flux_dn_loc, flux_dir_loc)
372 : !$omp target exit data map(from:flux_up_loc, flux_dn_loc, flux_dir_loc)
373 0 : if(.not. associated(fluxes%flux_up )) deallocate(flux_up_loc)
374 389263 : if(.not. associated(fluxes%flux_dn )) deallocate(flux_dn_loc)
375 778526 : if(.not. associated(fluxes%flux_dn_dir)) deallocate(flux_dir_loc)
376 : class default
377 389263 : deallocate(gpt_flux_up, gpt_flux_dn, gpt_flux_dir)
378 : end select
379 778526 : if(.not. has_dif_bc) then
380 : !$acc exit data delete( inc_flux_diffuse)
381 : !$omp target exit data map(release:inc_flux_diffuse)
382 778526 : deallocate(inc_flux_diffuse)
383 : end if
384 :
385 778526 : end function rte_sw_mu0_full
386 : !--------------------------------------------------------------------------------------------------------------------
387 : !
388 : ! Expand from band to g-point dimension, transpose dimensions (nband, ncol) -> (ncol,ngpt)
389 : !
390 1557052 : subroutine expand_and_transpose(ops,arr_in,arr_out)
391 : class(ty_optical_props), intent(in ) :: ops
392 : real(wp), dimension(:,:), intent(in ) :: arr_in ! (nband, ncol)
393 : real(wp), dimension(:,:), intent(out) :: arr_out ! (ncol, igpt)
394 : ! -------------
395 : integer :: ncol, nband, ngpt
396 : integer :: icol, iband, igpt
397 1557052 : integer, dimension(2,ops%get_nband()) :: limits
398 :
399 1557052 : ncol = size(arr_in, 2)
400 1557052 : nband = ops%get_nband()
401 1557052 : ngpt = ops%get_ngpt()
402 66953236 : limits = ops%get_band_lims_gpoint()
403 : !$acc parallel loop collapse(2) copyin(arr_in, limits)
404 : !$omp target teams distribute parallel do simd collapse(2) map(to:arr_in, limits)
405 23355780 : do iband = 1, nband
406 351308580 : do icol = 1, ncol
407 2973373928 : do igpt = limits(1, iband), limits(2, iband)
408 2951575200 : arr_out(icol, igpt) = arr_in(iband,icol)
409 : end do
410 : end do
411 : end do
412 :
413 1557052 : end subroutine expand_and_transpose
414 3114104 : end module mo_rte_sw
|
