LCOV - code coverage report
Current view: | top level - physics/rrtmgp/ext/rte-frontend - mo_rte_sw.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 92 | 118 | 78.0 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 3 | 3 | 100.0 %

              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

Generated by: LCOV version 1.14