diff --git a/README.md b/README.md index e6d624944..15128baa1 100644 --- a/README.md +++ b/README.md @@ -8,6 +8,8 @@ RTE computes fluxes given spectrally-resolved optical descriptions and source fu Example programs and documenation are evolving - please see examples/ in the repo and Wiki on the project's Github page. Suggestions are welcome. Meanwhile for questions please contact Robert Pincus and Eli Mlawer at rrtmgp@aer.com. +In the most recent revision, the default method for solution for longwave problems that include scattering has been changed from 2-stream methods to a re-scaled and refined no-scattering calculation following [Tang et al. 2018](https://doi.org/10.1175/JAS-D-18-0014.1). + ## Building the libraries. 1. `cd build` @@ -17,4 +19,4 @@ Example programs and documenation are evolving - please see examples/ in the rep ## Examples -Two examples are provided, one for clear skies and one including clouds. See the README file and codes in each directory for further information. +Two examples are provided, one for clear skies and one including clouds. See the README file and codes in each directory for further information. diff --git a/azure-pipelines.yml b/azure-pipelines.yml index 62bc80471..c775826e0 100644 --- a/azure-pipelines.yml +++ b/azure-pipelines.yml @@ -1,71 +1,88 @@ jobs: - job: build_daint + pool: CSCS + strategy: + matrix: + pgi_19_9_gpu: + compiler_module: PGI/19.9.0 + accel_module: cudatoolkit + FCFLAGS: "-O3 -ta=tesla:cc60,cuda10.1 -Minfo -Mallocatable=03 -gopt -Minline,reshape,maxsize:40" + RTE_KERNELS: openacc + pgi_default_gpu: + compiler_module: pgi + accel_module: craype-accel-nvidia60 + # Generic accelerator flag + FCFLAGS: "-O3 -acc -Minfo -Mallocatable=03 -gopt" + RTE_KERNELS: openacc + pgi_19_9_cpu: + compiler_module: PGI/19.9.0 + accel_module: + # Error checking flags + FCFLAGS: "-Mallocatable=03 -Mstandard -Mbounds -Mchkptr -Kieee -Mchkstk" + pgi_default_cpu: + compiler_module: pgi + accel_module: + # Error checking flags + FCFLAGS: "-Mallocatable=03 -Mstandard -Mbounds -Mchkptr -Kieee -Mchkstk" + maxParallel: 2 workspace: clean: all steps: - - script: | set -e - # This module will unload some of the build modules, so load the files separately + + echo " + export PATH=$CRAY_BINUTILS_BIN:$PATH module load daint-gpu - module load cray-python/3.6.5.3 netcdf-python + module swap PrgEnv-cray PrgEnv-pgi + module load cray-netcdf cray-hdf5 + module swap pgi $(compiler_module) + module load $(accel_module) + module load cray-python/3.6.5.7 + export LD_LIBRARY_PATH=$CRAY_LD_LIBRARY_PATH:$LD_LIBRARY_PATH + + echo Environment: + module list + " > modules + + displayName: 'Create module environment' + - script: | + set -e + source modules + # This module will unload some of the build modules, so load the files separately + module load netcdf-python cd examples/rfmip-clear-sky python ./stage_files.py displayName: 'Stage files' - - script: | set -e - echo "Initial environment:" - module list - - export PATH=$CRAY_BINUTILS_BIN:$PATH - module load daint-gpu - module swap PrgEnv-cray PrgEnv-pgi - module load cdt/19.08 - module unload pgi - module load PGI/19.9.0 - module load cray-netcdf cray-hdf5 cudatoolkit/9.2.148_3.19-6.0.7.1_2.1__g3d9acc8 - module load cray-python/3.6.5.3 - export LD_LIBRARY_PATH=$CRAY_LD_LIBRARY_PATH:$LD_LIBRARY_PATH - - echo "*********************************************" - echo "Building environment:" - module list - - export RTE_KERNELS="openacc" + source modules export RRTMGP_ROOT=$PWD export RRTMGP_DIR=$PWD/build - export FC=ftn - export FCFLAGS="-O3 -ta=tesla:cc60,cuda9.2 -Minfo -Mallocatable=03 -Mpreprocess -gopt -Minline,reshape,maxsize:40" - make -C build/ -j 8 make -C examples/all-sky -j 8 make -C examples/rfmip-clear-sky -j 8 - + displayName: 'Make' + - script: | + set -e + source modules cd examples/rfmip-clear-sky - srun -C gpu -A c15 -p cscsci python ./run-rfmip-examples.py + srun -C gpu -A c15 -p cscsci python ./run-rfmip-examples.py --block_size 1800 cd ../.. - cd examples/all-sky srun -C gpu -A c15 -p cscsci python ./run-allsky-example.py - displayName: 'Make & run' - + displayName: 'Run' - script: | set -e - - module load daint-gpu + source modules # This module will unload some of the build modules, so do the checks separately - module load cray-python/3.6.5.3 netcdf-python - + module load netcdf-python cd examples/rfmip-clear-sky python ./compare-to-reference.py --fail=1.e-4 cd ../.. - cd examples/all-sky python ./compare-to-reference.py displayName: 'Check results' - - pool: CSCS diff --git a/rte/kernels-openacc/mo_rte_solver_kernels.F90 b/rte/kernels-openacc/mo_rte_solver_kernels.F90 index 30494214b..19384bf6a 100644 --- a/rte/kernels-openacc/mo_rte_solver_kernels.F90 +++ b/rte/kernels-openacc/mo_rte_solver_kernels.F90 @@ -39,6 +39,7 @@ module mo_rte_solver_kernels lw_solver_noscat, lw_solver_noscat_GaussQuad, lw_solver_2stream, & sw_solver_noscat, sw_solver_2stream + public :: lw_solver_1rescl_GaussQuad, lw_solver_1rescl ! These routines don't really need to be visible but making them so is useful for testing. public :: lw_source_noscat, lw_combine_sources, & lw_source_2str, sw_source_2str, & @@ -1197,5 +1198,409 @@ subroutine apply_BC_0(ncol, nlay, ngpt, top_at_1, flux_dn) bind (C, name="apply_ end do end if end subroutine apply_BC_0 +! ------------------------------------------------------------------------------------------------- +! +! Similar to Longwave no-scattering (lw_solver_noscat) +! a) relies on rescaling of the optical parameters based on asymetry factor and single scattering albedo +! scaling can be computed by scaling_1rescl +! b) adds adustment term based on cloud properties (lw_transport_1rescl) +! adustment terms is computed based on solution of the Tang equations +! for "linear-in-tau" internal source (not in the paper) +! +! Attention: +! use must prceompute scaling before colling the function +! +! Implemented based on the paper +! Tang G, et al, 2018: https://doi.org/10.1175/JAS-D-18-0014.1 +! +! ------------------------------------------------------------------------------------------------- + subroutine lw_solver_1rescl(ncol, nlay, ngpt, top_at_1, D, & + tau, scaling, lay_source, lev_source_inc, lev_source_dec, sfc_emis, sfc_src, & + radn_up, radn_dn) bind(C, name="lw_solver_1rescl") + integer, intent(in ) :: ncol, nlay, ngpt ! Number of columns, layers, g-points + logical(wl), intent(in ) :: top_at_1 + real(wp), dimension(ncol, ngpt), intent(in ) :: D ! secant of propagation angle [] + real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: tau ! Absorption optical thickness [] + real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: scaling ! single scattering albedo [] + real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: lay_source ! Planck source at layer average temperature [W/m2] + ! Planck source at layer edge for radiation in increasing/decreasing ilay direction + ! lev_source_dec applies the mapping in layer i to the Planck function at layer i + ! lev_source_inc applies the mapping in layer i to the Planck function at layer i+1 + real(wp), dimension(ncol,nlay, ngpt), target, & + intent(in ) :: lev_source_inc, lev_source_dec + real(wp), dimension(ncol, ngpt), intent(in ) :: sfc_emis ! Surface emissivity [] + real(wp), dimension(ncol, ngpt), intent(in ) :: sfc_src ! Surface source function [W/m2] + real(wp), dimension(ncol,nlay+1,ngpt), intent( out) :: radn_up ! Radiances [W/m2-str] + real(wp), dimension(ncol,nlay+1,ngpt), intent(inout) :: radn_dn ! Top level must contain incident flux boundary condition + + ! Local variables, WITH g-point dependency + real(wp), dimension(ncol,nlay,ngpt) :: tau_loc, & ! path length (tau/mu) + trans ! transmissivity = exp(-tau) + real(wp), dimension(ncol,nlay,ngpt) :: source_dn, source_up + real(wp), dimension(ncol, ngpt) :: source_sfc, sfc_albedo + + real(wp), dimension(:,:,:), pointer :: lev_source_up, lev_source_dn ! Mapping increasing/decreasing indicies to up/down + + real(wp), parameter :: pi = acos(-1._wp) + integer :: ilev, igpt, top_level + ! ------------------------------------ + real(wp) :: fact + real(wp), parameter :: tau_thresh = sqrt(epsilon(tau)) + integer :: icol + real(wp), dimension(ncol ) :: sfcSource + real(wp), dimension(ncol,nlay,ngpt) :: An, Cn + ! ------------------------------------ + + ! Which way is up? + ! Level Planck sources for upward and downward radiation + ! When top_at_1, lev_source_up => lev_source_dec + ! lev_source_dn => lev_source_inc, and vice-versa + if(top_at_1) then + top_level = 1 + lev_source_up => lev_source_dec + lev_source_dn => lev_source_inc + else + top_level = nlay+1 + lev_source_up => lev_source_inc + lev_source_dn => lev_source_dec + end if + + !$acc enter data copyin(d,tau,sfc_src,sfc_emis,lev_source_dec,lev_source_inc,lay_source,radn_dn) + !$acc enter data copyin(scaling) + !$acc enter data create(tau_loc,trans,source_dn,source_up,source_sfc,sfc_albedo,radn_up) + !$acc enter data create(sfcSource, An, Cn) + !$acc enter data attach(lev_source_up,lev_source_dn) + + ! NOTE: This kernel produces small differences between GPU and CPU + ! implementations on Ascent with PGI, we assume due to floating point + ! differences in the exp() function. These differences are small in the + ! RFMIP test case (10^-6). + !$acc parallel loop collapse(3) + do igpt = 1, ngpt + do ilev = 1, nlay + do icol = 1, ncol + ! + ! Optical path and transmission, used in source function and transport calculations + ! + tau_loc(icol,ilev,igpt) = tau(icol,ilev,igpt)*D(icol,igpt) + trans (icol,ilev,igpt) = exp(-tau_loc(icol,ilev,igpt)) + ! here scaling is used to store parameter wb/[(]1-w(1-b)] of Eq.21 of the Tang's paper + ! explanation of factor 0.4 note A of Table + Cn(icol,ilev,igpt) = 0.4_wp*scaling(icol,ilev,igpt) + An(icol,ilev,igpt) = (1._wp-trans(icol,ilev,igpt)*trans(icol,ilev,igpt)) + end do + end do + end do + + !$acc parallel loop collapse(2) + do igpt = 1, ngpt + do icol = 1, ncol + ! + ! Surface albedo, surface source function + ! + sfc_albedo(icol,igpt) = 1._wp - sfc_emis(icol,igpt) + source_sfc(icol,igpt) = sfc_emis(icol,igpt) * sfc_src(icol,igpt) + end do + end do + + + ! + ! Source function for diffuse radiation + ! + call lw_source_noscat(ncol, nlay, ngpt, & + lay_source, lev_source_up, lev_source_dn, & + tau_loc, trans, source_dn, source_up) + + ! + ! Transport + ! + ! compute no-scattering fluxes + call lw_transport_noscat(ncol, nlay, ngpt, top_at_1, & + tau_loc, trans, sfc_albedo, source_dn, source_up, source_sfc, & + radn_up, radn_dn) + ! make adjustment + call lw_transport_1rescl(ncol, nlay, ngpt, top_at_1, & + tau_loc, trans, & + sfc_albedo, source_dn, source_up, source_sfc, & + radn_up, radn_dn, An, Cn) + !$acc exit data copyout(radn_dn,radn_up) + !$acc exit data delete(sfcSource, An, Cn) + !$acc exit data delete(scaling) + !$acc exit data delete(d,tau,sfc_src,sfc_emis,lev_source_dec,lev_source_inc,lay_source,tau_loc,trans,source_dn,source_up,source_sfc,sfc_albedo) + !$acc exit data detach(lev_source_up,lev_source_dn) + + end subroutine lw_solver_1rescl +! ------------------------------------------------------------------------------------------------- +! +! Similar to lw_solver_noscat_GaussQuad. +! It is main solver to use the Tang approximation for fluxes +! In addition to the no scattering input parameters the user must provide +! scattering related properties (ssa and g) that the solver uses to compute scaling +! +! --------------------------------------------------------------- + subroutine lw_solver_1rescl_GaussQuad(ncol, nlay, ngpt, top_at_1, nmus, Ds, weights, & + tau, ssa, g, lay_source, lev_source_inc, lev_source_dec, & + sfc_emis, sfc_src,& + flux_up, flux_dn) & + bind(C, name="lw_solver_1rescl_GaussQuad") + integer, intent(in ) :: ncol, nlay, ngpt ! Number of columns, layers, g-points + logical(wl), intent(in ) :: top_at_1 + integer, intent(in ) :: nmus ! number of quadrature angles + real(wp), dimension(nmus), intent(in ) :: Ds, weights ! quadrature secants, weights + real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: tau ! Optical thickness, + real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: ssa ! single-scattering albedo, + real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: g ! asymmetry parameter [] + real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: lay_source ! Planck source at layer average temperature [W/m2] + real(wp), dimension(ncol,nlay+1,ngpt), intent(in ) :: lev_source_inc + ! Planck source at layer edge for radiation in increasing ilay direction [W/m2] + ! Includes spectral weighting that accounts for state-dependent frequency to g-space mapping + real(wp), dimension(ncol,nlay+1,ngpt), intent(in ) :: lev_source_dec + ! Planck source at layer edge for radiation in decreasing ilay direction [W/m2] + real(wp), dimension(ncol, ngpt), intent(in ) :: sfc_emis ! Surface emissivity [] + real(wp), dimension(ncol, ngpt), intent(in ) :: sfc_src ! Surface source function [W/m2] + real(wp), dimension(ncol,nlay+1,ngpt), intent( out) :: flux_up ! Radiances [W/m2-str] + real(wp), dimension(ncol,nlay+1,ngpt), intent(inout) :: flux_dn ! Top level must contain incident flux boundary condition + ! Local variables + real(wp), dimension(ncol,nlay+1,ngpt) :: radn_dn, radn_up ! Fluxes per quad angle + real(wp), dimension(ncol, ngpt) :: Ds_ncol + + integer :: imu, top_level,icol,ilev,igpt + real :: weight + + real(wp), dimension(ncol,ngpt) :: fluxTOA ! downward flux at TOA + real(wp), dimension(ncol,nlay, ngpt) :: tauLoc ! rescaled Tau + real(wp), dimension(ncol,nlay, ngpt) :: scaling ! scaling + real(wp), parameter :: tresh=1.0_wp - 1e-6_wp + + !$acc enter data copyin(Ds,weights,tau,ssa,g,lay_source,lev_source_inc,lev_source_dec,sfc_emis,sfc_src,flux_dn) + !$acc enter data create(flux_up,radn_dn,radn_up,Ds_ncol, scaling, tauLoc) + + + ! Tang rescaling + if (any(ssa*g >= tresh)) then + call scaling_1rescl_safe(ncol, nlay, ngpt, tauLoc, scaling, tau, ssa, g) + else + call scaling_1rescl(ncol, nlay, ngpt, tauLoc, scaling, tau, ssa, g) + endif + + ! ------------------------------------ + ! + ! For the first angle output arrays store total flux + ! + top_level = MERGE(1, nlay+1, top_at_1) + ! store TOA flux + fluxTOA = flux_dn(1:ncol, top_level, 1:ngpt) + + Ds_ncol(:,:) = Ds(1) + weight = 2._wp*pi*weights(1) + ! Transport is for intensity + ! convert flux at top of domain to intensity assuming azimuthal isotropy + ! + radn_dn(1:ncol, top_level, 1:ngpt) = fluxTOA(1:ncol, 1:ngpt) / weight + + call lw_solver_1rescl(ncol, nlay, ngpt, & + top_at_1, Ds_ncol, tauLoc, scaling, & + lay_source, lev_source_inc, lev_source_dec, sfc_emis, sfc_src, & + flux_up, flux_dn) + !$acc parallel loop collapse(3) + do igpt = 1, ngpt + do ilev = 1, nlay+1 + do icol = 1, ncol + flux_up(icol,ilev,igpt) = weight*flux_up(icol,ilev,igpt) + flux_dn(icol,ilev,igpt) = weight*flux_dn(icol,ilev,igpt) + enddo + enddo + enddo + + do imu = 2, nmus + Ds_ncol(:,:) = Ds(imu) + weight = 2._wp*pi*weights(imu) + radn_dn(1:ncol, top_level, 1:ngpt) = fluxTOA(1:ncol, 1:ngpt) / weight + call lw_solver_1rescl(ncol, nlay, ngpt, & + top_at_1, Ds_ncol, tauLoc, scaling, & + lay_source, lev_source_inc, lev_source_dec, sfc_emis, sfc_src, & + radn_up, radn_dn) + + !$acc parallel loop collapse(3) + do igpt = 1, ngpt + do ilev = 1, nlay+1 + do icol = 1, ncol + flux_up(icol,ilev,igpt) = flux_up(icol,ilev,igpt) + weight*radn_up(icol,ilev,igpt) + flux_dn(icol,ilev,igpt) = flux_dn(icol,ilev,igpt) + weight*radn_dn(icol,ilev,igpt) + enddo + enddo + enddo + + end do + !$acc exit data copyout(flux_up,flux_dn) + !$acc exit data delete(Ds,weights,tau,ssa,g,tauLoc,scaling,lay_source,lev_source_inc,lev_source_dec,sfc_emis,sfc_src,radn_dn,radn_up,Ds_ncol) + end subroutine lw_solver_1rescl_GaussQuad +! ------------------------------------------------------------------------------------------------- +! +! Computes Tang scaling of layer optical thickness and scaling parameter +! unsafe if ssa*g =1. +! +! --------------------------------------------------------------- + pure subroutine scaling_1rescl(ncol, nlay, ngpt, tauLoc, scaling, tau, ssa, g) + integer , intent(in) :: ncol + integer , intent(in) :: nlay + integer , intent(in) :: ngpt + real(wp), dimension(ncol, nlay, ngpt), intent(in) :: tau + real(wp), dimension(ncol, nlay, ngpt), intent(in) :: ssa + real(wp), dimension(ncol, nlay, ngpt), intent(in) :: g + + real(wp), dimension(ncol, nlay, ngpt), intent(inout) :: tauLoc + real(wp), dimension(ncol, nlay, ngpt), intent(inout) :: scaling + + integer :: icol, ilay, igpt + real(wp) :: wb, ssal, scaleTau + !$acc enter data copyin(tau, ssa, g) + !$acc enter data create(tauLoc, scaling) + !$acc parallel loop collapse(3) + do igpt=1,ngpt + do ilay=1,nlay + do icol=1,ncol + ssal = ssa(icol, ilay, igpt) + wb = ssal*(1._wp - g(icol, ilay, igpt)) / 2._wp + scaleTau = (1._wp - ssal + wb ) + + tauLoc(icol, ilay, igpt) = scaleTau * tau(icol, ilay, igpt) ! Eq.15 of the paper + ! + ! here scaling is used to store parameter wb/[1-w(1-b)] of Eq.21 of the Tang's paper + ! actually it is in line of parameter rescaling defined in Eq.7 + ! potentialy if g=ssa=1 then wb/scaleTau = NaN + ! it should not happen + scaling(icol, ilay, igpt) = wb / scaleTau + enddo + enddo + enddo + !$acc exit data copyout(tauLoc, scaling) + !$acc exit data delete(tau, ssa, g) + end subroutine scaling_1rescl +! ------------------------------------------------------------------------------------------------- +! +! Computes Tang scaling of layer optical thickness and scaling parameter +! Safe implementation +! +! --------------------------------------------------------------- + pure subroutine scaling_1rescl_safe(ncol, nlay, ngpt, tauLoc, scaling, tau, ssa, g) + integer , intent(in) :: ncol + integer , intent(in) :: nlay + integer , intent(in) :: ngpt + real(wp), dimension(ncol, nlay, ngpt), intent(in) :: tau + real(wp), dimension(ncol, nlay, ngpt), intent(in) :: ssa + real(wp), dimension(ncol, nlay, ngpt), intent(in) :: g + + real(wp), dimension(ncol, nlay, ngpt), intent(inout) :: tauLoc + real(wp), dimension(ncol, nlay, ngpt), intent(inout) :: scaling + + integer :: icol, ilay, igpt + real(wp) :: wb, ssal, scaleTau + !$acc enter data copyin(tau, ssa, g) + !$acc enter data create(tauLoc, scaling) + !$acc parallel loop collapse(3) + do igpt=1,ngpt + do ilay=1,nlay + do icol=1,ncol + ssal = ssa(icol, ilay, igpt) + wb = ssal*(1._wp - g(icol, ilay, igpt)) / 2._wp + scaleTau = (1._wp - ssal + wb ) + + tauLoc(icol, ilay, igpt) = scaleTau * tau(icol, ilay, igpt) ! Eq.15 of the paper + ! + ! here scaling is used to store parameter wb/[1-w(1-b)] of Eq.21 of the Tang's paper + ! actually it is in line of parameter rescaling defined in Eq.7 + if (scaleTau < 1e-6_wp) then + scaling(icol, ilay, igpt) = 1.0_wp + else + scaling(icol, ilay, igpt) = wb / scaleTau + endif + enddo + enddo + enddo + !$acc exit data copyout(tauLoc, scaling) + !$acc exit data delete(tau, ssa, g) + end subroutine scaling_1rescl_safe +! ------------------------------------------------------------------------------------------------- +! +! Similar to Longwave no-scattering tarnsport (lw_transport_noscat) +! a) adds adjustment factor based on cloud properties +! +! implementation notice: +! the adjustmentFactor computation can be skipped where Cn <= epsilon +! +! ------------------------------------------------------------------------------------------------- + subroutine lw_transport_1rescl(ncol, nlay, ngpt, top_at_1, & + tau, trans, sfc_albedo, source_dn, source_up, source_sfc, & + radn_up, radn_dn, An, Cn) bind(C, name="lw_transport_1rescl") + integer, intent(in ) :: ncol, nlay, ngpt ! Number of columns, layers, g-points + logical(wl), intent(in ) :: top_at_1 ! + real(wp), dimension(ncol,nlay ,ngpt), intent(in ) :: tau, & ! Absorption optical thickness, pre-divided by mu [] + trans ! transmissivity = exp(-tau) + real(wp), dimension(ncol ,ngpt), intent(in ) :: sfc_albedo ! Surface albedo + real(wp), dimension(ncol,nlay ,ngpt), intent(in ) :: source_dn, & + source_up ! Diffuse radiation emitted by the layer + real(wp), dimension(ncol ,ngpt), intent(in ) :: source_sfc ! Surface source function [W/m2] + real(wp), dimension(ncol,nlay+1,ngpt), intent(inout) :: radn_up ! Radiances [W/m2-str] + real(wp), dimension(ncol,nlay+1,ngpt), intent(inout) :: radn_dn !Top level must contain incident flux boundary condition + real(wp), dimension(ncol,nlay ,ngpt), intent(in ) :: An, Cn + ! Local variables + integer :: ilev, icol, igpt + ! --------------------------------------------------- + real(wp) :: adjustmentFactor + if(top_at_1) then + ! + ! Top of domain is index 1 + ! + ! Downward propagation + !$acc parallel loop collapse(2) + do igpt = 1, ngpt + do icol = 1, ncol + ! 1st Upward propagation + do ilev = nlay, 1, -1 + radn_up(icol,ilev,igpt) = trans(icol,ilev,igpt)*radn_up(icol,ilev+1,igpt) + source_up(icol,ilev,igpt) + adjustmentFactor = Cn(icol,ilev,igpt)*& + ( An(icol,ilev,igpt)*radn_dn(icol,ilev,igpt) - & + source_dn(icol,ilev,igpt) *trans(icol,ilev,igpt ) - & + source_up(icol,ilev,igpt)) + radn_up(icol,ilev,igpt) = radn_up(icol,ilev,igpt) + adjustmentFactor + enddo + ! 2nd Downward propagation + do ilev = 1, nlay + radn_dn(icol,ilev+1,igpt) = trans(icol,ilev,igpt)*radn_dn(icol,ilev,igpt) + source_dn(icol,ilev,igpt) + adjustmentFactor = Cn(icol,ilev,igpt)*( & + An(icol,ilev,igpt)*radn_up(icol,ilev,igpt) - & + source_up(icol,ilev,igpt)*trans(icol,ilev,igpt) - & + source_dn(icol,ilev,igpt) ) + radn_dn(icol,ilev+1,igpt) = radn_dn(icol,ilev+1,igpt) + adjustmentFactor + enddo + enddo + enddo + else + !$acc parallel loop collapse(2) + do igpt = 1, ngpt + do icol = 1, ncol + ! Upward propagation + do ilev = 1, nlay + radn_up(icol,ilev+1,igpt) = trans(icol,ilev,igpt) * radn_up(icol,ilev,igpt) + source_up(icol,ilev,igpt) + adjustmentFactor = Cn(icol,ilev,igpt)*& + ( An(icol,ilev,igpt)*radn_dn(icol,ilev+1,igpt) - & + source_dn(icol,ilev,igpt) *trans(icol,ilev ,igpt) - & + source_up(icol,ilev,igpt)) + radn_up(icol,ilev+1,igpt) = radn_up(icol,ilev+1,igpt) + adjustmentFactor + end do + ! 2st Downward propagation + do ilev = nlay, 1, -1 + radn_dn(icol,ilev,igpt) = trans(icol,ilev,igpt)*radn_dn(icol,ilev+1,igpt) + source_dn(icol,ilev,igpt) + adjustmentFactor = Cn(icol,ilev,igpt)*( & + An(icol,ilev,igpt)*radn_up(icol,ilev,igpt) - & + source_up(icol,ilev,igpt)*trans(icol,ilev ,igpt ) - & + source_dn(icol,ilev,igpt) ) + radn_dn(icol,ilev,igpt) = radn_dn(icol,ilev,igpt) + adjustmentFactor + end do + enddo + enddo + end if + end subroutine lw_transport_1rescl end module mo_rte_solver_kernels diff --git a/rte/kernels/mo_rte_solver_kernels.F90 b/rte/kernels/mo_rte_solver_kernels.F90 index a243c261b..156ff027a 100644 --- a/rte/kernels/mo_rte_solver_kernels.F90 +++ b/rte/kernels/mo_rte_solver_kernels.F90 @@ -39,6 +39,8 @@ module mo_rte_solver_kernels lw_solver_noscat, lw_solver_noscat_GaussQuad, lw_solver_2stream, & sw_solver_noscat, sw_solver_2stream + public :: lw_solver_1rescl_GaussQuad, lw_solver_1rescl + ! These routines don't really need to be visible but making them so is useful for testing. public :: lw_source_noscat, lw_combine_sources, & lw_source_2str, sw_source_2str, & @@ -947,5 +949,342 @@ pure subroutine apply_BC_0(ncol, nlay, ngpt, top_at_1, flux_dn) bind (C, name="a flux_dn(1:ncol, nlay+1, 1:ngpt) = 0._wp end if end subroutine apply_BC_0 +! ------------------------------------------------------------------------------------------------- +! +! Similar to Longwave no-scattering (lw_solver_noscat) +! a) relies on rescaling of the optical parameters based on asymetry factor and single scattering albedo +! scaling can be computed by scaling_1rescl +! b) adds adustment term based on cloud properties (lw_transport_1rescl) +! adustment terms is computed based on solution of the Tang equations +! for "linear-in-tau" internal source (not in the paper) +! +! Attention: +! use must prceompute scaling before colling the function +! +! Implemented based on the paper +! Tang G, et al, 2018: https://doi.org/10.1175/JAS-D-18-0014.1 +! +! ------------------------------------------------------------------------------------------------- +subroutine lw_solver_1rescl(ncol, nlay, ngpt, top_at_1, D, & + tau, scaling, lay_source, lev_source_inc, lev_source_dec, sfc_emis, sfc_src, & + radn_up, radn_dn) bind(C, name="lw_solver_1rescl") + integer, intent(in ) :: ncol, nlay, ngpt ! Number of columns, layers, g-points + logical(wl), intent(in ) :: top_at_1 + real(wp), dimension(ncol, ngpt), intent(in ) :: D ! secant of propagation angle [] + real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: tau ! Absorption optical thickness [] + real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: scaling ! single scattering albedo [] + real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: lay_source ! Planck source at layer average temperature [W/m2] + ! Planck source at layer edge for radiation in increasing/decreasing ilay direction + ! lev_source_dec applies the mapping in layer i to the Planck function at layer i + ! lev_source_inc applies the mapping in layer i to the Planck function at layer i+1 + real(wp), dimension(ncol,nlay, ngpt), target, & + intent(in ) :: lev_source_inc, lev_source_dec + real(wp), dimension(ncol, ngpt), intent(in ) :: sfc_emis ! Surface emissivity [] + real(wp), dimension(ncol, ngpt), intent(in ) :: sfc_src ! Surface source function [W/m2] + real(wp), dimension(ncol,nlay+1,ngpt), intent( out) :: radn_up ! Radiances [W/m2-str] + real(wp), dimension(ncol,nlay+1,ngpt), intent(inout) :: radn_dn ! Top level must contain incident flux boundary condition + + ! Local variables, no g-point dependency + real(wp), dimension(ncol,nlay) :: tau_loc, & ! path length (tau/mu) + trans ! transmissivity = exp(-tau) + real(wp), dimension(ncol,nlay) :: source_dn, source_up + real(wp), dimension(ncol ) :: source_sfc, sfc_albedo + real(wp), dimension(ncol,nlay) :: An, Cn + + real(wp), dimension(:,:,:), pointer :: lev_source_up, lev_source_dn ! Mapping increasing/decreasing indicies to up/down + + real(wp), parameter :: pi = acos(-1._wp) + integer :: ilev, igpt, top_level + ! ------------------------------------ + real(wp), parameter :: tau_thresh = sqrt(epsilon(tau)) + ! ------------------------------------ + + ! Which way is up? + ! Level Planck sources for upward and downward radiation + ! When top_at_1, lev_source_up => lev_source_dec + ! lev_source_dn => lev_source_inc, and vice-versa + if(top_at_1) then + top_level = 1 + lev_source_up => lev_source_dec + lev_source_dn => lev_source_inc + else + top_level = nlay+1 + lev_source_up => lev_source_inc + lev_source_dn => lev_source_dec + end if + + do igpt = 1, ngpt + ! + ! Optical path and transmission, used in source function and transport calculations + ! + do ilev = 1, nlay + tau_loc(:,ilev) = tau(:,ilev,igpt)*D(:,igpt) + trans (:,ilev) = exp(-tau_loc(:,ilev)) + ! + ! here scaling is used to store parameter wb/(1-w(1-b)) of Eq.21 of the Tang's paper + ! explanation of factor 0.4 note A of Table + ! + Cn(:,ilev) = 0.4_wp*scaling(:,ilev,igpt) + An(:,ilev) = (1._wp-trans(:,ilev)*trans(:,ilev)) + end do + + ! Source function for diffuse radiation + ! + call lw_source_noscat(ncol, nlay, & + lay_source(:,:,igpt), lev_source_up(:,:,igpt), lev_source_dn(:,:,igpt), & + tau_loc, trans, source_dn, source_up) + + ! + ! Surface albedo, surface source function + ! + sfc_albedo(:) = 1._wp - sfc_emis(:,igpt) + source_sfc(:) = sfc_emis(:,igpt) * sfc_src(:,igpt) + ! + ! Transport + ! + ! compute no-scattering fluxes + call lw_transport_noscat(ncol, nlay, top_at_1, & + tau_loc, trans, sfc_albedo, source_dn, source_up, source_sfc, & + radn_up(:,:,igpt), radn_dn(:,:,igpt)) + + ! make adjustment + call lw_transport_1rescl(ncol, nlay, top_at_1, trans, & + sfc_albedo, source_dn, source_up, source_sfc, & + radn_up(:,:,igpt), radn_dn(:,:,igpt), An, Cn) + + end do ! g point loop +end subroutine lw_solver_1rescl +! ------------------------------------------------------------------------------------------------- +! +! Similar to lw_solver_noscat_GaussQuad. +! It is main solver to use the rescaled-for-scattering approximation for fluxes +! In addition to the no scattering input parameters the user must provide +! scattering related properties (ssa and g) that the solver uses to compute scaling +! +! --------------------------------------------------------------- +subroutine lw_solver_1rescl_GaussQuad(ncol, nlay, ngpt, top_at_1, nmus, Ds, weights, & + tau, ssa, g, lay_source, lev_source_inc, lev_source_dec, & + sfc_emis, sfc_src,& + flux_up, flux_dn) & + bind(C, name="lw_solver_1rescl_GaussQuad") + integer, intent(in ) :: ncol, nlay, ngpt ! Number of columns, layers, g-points + logical(wl), intent(in ) :: top_at_1 + integer, intent(in ) :: nmus ! number of quadrature angles + real(wp), dimension(nmus), intent(in ) :: Ds, weights ! quadrature secants, weights + real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: tau ! Optical thickness, + real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: ssa ! single-scattering albedo, + real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: g ! asymmetry parameter [] + real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: lay_source ! Planck source at layer average temperature [W/m2] + real(wp), dimension(ncol,nlay+1,ngpt), intent(in ) :: lev_source_inc + ! Planck source at layer edge for radiation in increasing ilay direction [W/m2] + ! Includes spectral weighting that accounts for state-dependent frequency to g-space mapping + real(wp), dimension(ncol,nlay+1,ngpt), intent(in ) :: lev_source_dec + ! Planck source at layer edge for radiation in decreasing ilay direction [W/m2] + real(wp), dimension(ncol, ngpt), intent(in ) :: sfc_emis ! Surface emissivity [] + real(wp), dimension(ncol, ngpt), intent(in ) :: sfc_src ! Surface source function [W/m2] + real(wp), dimension(ncol,nlay+1,ngpt), intent( out) :: flux_up ! Radiances [W/m2-str] + real(wp), dimension(ncol,nlay+1,ngpt), intent(inout) :: flux_dn ! Top level must contain incident flux boundary condition + ! Local variables + real(wp), dimension(ncol,nlay+1,ngpt) :: radn_dn, radn_up ! Fluxes per quad angle + real(wp), dimension(ncol, ngpt) :: Ds_ncol + + real(wp), dimension(ncol,nlay, ngpt) :: tauLoc ! rescaled Tau + real(wp), dimension(ncol,nlay, ngpt) :: scaling ! scaling + real(wp), dimension(ncol,ngpt) :: fluxTOA ! downward flux at TOA + + integer :: imu, top_level + real :: weight + real(wp), parameter :: tresh=1.0_wp - 1e-6_wp + + ! Tang rescaling + if (any(ssa*g >= tresh)) then + call scaling_1rescl_safe(ncol, nlay, ngpt, tauLoc, scaling, tau, ssa, g) + else + call scaling_1rescl(ncol, nlay, ngpt, tauLoc, scaling, tau, ssa, g) + endif + ! ------------------------------------ + ! + ! For the first angle output arrays store total flux + ! + top_level = MERGE(1, nlay+1, top_at_1) + fluxTOA = flux_dn(1:ncol, top_level, 1:ngpt) + Ds_ncol(:,:) = Ds(1) + weight = 2._wp*pi*weights(1) + ! Transport is for intensity + ! convert flux at top of domain to intensity assuming azimuthal isotropy + ! + radn_dn(1:ncol, top_level, 1:ngpt) = fluxTOA(1:ncol, 1:ngpt) / weight + call lw_solver_1rescl(ncol, nlay, ngpt, & + top_at_1, Ds_ncol, tauLoc, scaling, & + lay_source, lev_source_inc, lev_source_dec, sfc_emis, sfc_src, & + flux_up, flux_dn) + + flux_up = flux_up * weight + flux_dn = flux_dn * weight + + do imu = 2, nmus + Ds_ncol(:,:) = Ds(imu) + weight = 2._wp*pi*weights(imu) + ! Transport is for intensity + ! convert flux at top of domain to intensity assuming azimuthal isotropy + ! + radn_dn(1:ncol, top_level, 1:ngpt) = fluxTOA(1:ncol, 1:ngpt) / weight + call lw_solver_1rescl(ncol, nlay, ngpt, & + top_at_1, Ds_ncol, tauLoc, scaling, & + lay_source, lev_source_inc, lev_source_dec, sfc_emis, sfc_src, & + radn_up, radn_dn) + + flux_up(:,:,:) = flux_up(:,:,:) + weight*radn_up(:,:,:) + flux_dn(:,:,:) = flux_dn(:,:,:) + weight*radn_dn(:,:,:) + end do + end subroutine lw_solver_1rescl_GaussQuad +! ------------------------------------------------------------------------------------------------- +! +! Computes re-scaled layer optical thickness and scaling parameter +! unsafe if ssa*g =1. +! +! --------------------------------------------------------------- + pure subroutine scaling_1rescl(ncol, nlay, ngpt, tauLoc, scaling, tau, ssa, g) + integer , intent(in) :: ncol + integer , intent(in) :: nlay + integer , intent(in) :: ngpt + real(wp), dimension(ncol, nlay, ngpt), intent(in) :: tau + real(wp), dimension(ncol, nlay, ngpt), intent(in) :: ssa + real(wp), dimension(ncol, nlay, ngpt), intent(in) :: g + + real(wp), dimension(ncol, nlay, ngpt), intent(inout) :: tauLoc + real(wp), dimension(ncol, nlay, ngpt), intent(inout) :: scaling + + integer :: icol, ilay, igpt + real(wp) :: wb, ssal, scaleTau + do igpt=1,ngpt + do ilay=1,nlay + do icol=1,ncol + ssal = ssa(icol, ilay, igpt) + wb = ssal*(1._wp - g(icol, ilay, igpt)) / 2._wp + scaleTau = (1._wp - ssal + wb ) + tauLoc(icol, ilay, igpt) = scaleTau * tau(icol, ilay, igpt) ! Eq.15 of the paper + ! + ! here scaling is used to store parameter wb/(1-w(1-b)) of Eq.21 of the Tang paper + ! actually it is in line of parameter rescaling defined in Eq.7 + ! potentialy if g=ssa=1 then wb/scaleTau = NaN + ! it should not happen + scaling(icol, ilay, igpt) = wb / scaleTau + enddo + enddo + enddo + end subroutine scaling_1rescl +! ------------------------------------------------------------------------------------------------- +! +! Computes re-scaled layer optical thickness and scaling parameter +! safe implementation +! +! --------------------------------------------------------------- + pure subroutine scaling_1rescl_safe(ncol, nlay, ngpt, tauLoc, scaling, tau, ssa, g) + integer , intent(in) :: ncol + integer , intent(in) :: nlay + integer , intent(in) :: ngpt + real(wp), dimension(ncol, nlay, ngpt), intent(in) :: tau + real(wp), dimension(ncol, nlay, ngpt), intent(in) :: ssa + real(wp), dimension(ncol, nlay, ngpt), intent(in) :: g + + real(wp), dimension(ncol, nlay, ngpt), intent(inout) :: tauLoc + real(wp), dimension(ncol, nlay, ngpt), intent(inout) :: scaling + + integer :: icol, ilay, igpt + real(wp) :: wb, ssal, scaleTau + do igpt=1,ngpt + do ilay=1,nlay + do icol=1,ncol + ssal = ssa(icol, ilay, igpt) + wb = ssal*(1._wp - g(icol, ilay, igpt)) / 2._wp + scaleTau = (1._wp - ssal + wb ) + tauLoc(icol, ilay, igpt) = scaleTau * tau(icol, ilay, igpt) ! Eq.15 of the paper + ! + ! here scaling is used to store parameter wb/(1-w(1-b)) of Eq.21 of the Tang paper + ! actually it is in line of parameter rescaling defined in Eq.7 + if (scaleTau < 1e-6_wp) then + scaling(icol, ilay, igpt) = 1.0_wp + else + scaling(icol, ilay, igpt) = wb / scaleTau + endif + enddo + enddo + enddo + end subroutine scaling_1rescl_safe +! ------------------------------------------------------------------------------------------------- +! +! Similar to Longwave no-scattering tarnsport (lw_transport_noscat) +! a) adds adjustment factor based on cloud properties +! +! implementation notice: +! the adjustmentFactor computation can be skipped where Cn <= epsilon +! +! ------------------------------------------------------------------------------------------------- +subroutine lw_transport_1rescl(ncol, nlay, top_at_1, & + trans, sfc_albedo, source_dn, source_up, source_sfc, & + radn_up, radn_dn, An, Cn) bind(C, name="lw_transport_1rescl") + integer, intent(in ) :: ncol, nlay ! Number of columns, layers, g-points + logical(wl), intent(in ) :: top_at_1 ! + real(wp), dimension(ncol,nlay ), intent(in ) :: trans ! transmissivity = exp(-tau) + real(wp), dimension(ncol ), intent(in ) :: sfc_albedo ! Surface albedo + real(wp), dimension(ncol,nlay ), intent(in ) :: source_dn, & + source_up ! Diffuse radiation emitted by the layer + real(wp), dimension(ncol ), intent(in ) :: source_sfc ! Surface source function [W/m2] + real(wp), dimension(ncol,nlay+1), intent(inout) :: radn_up ! Radiances [W/m2-str] + real(wp), dimension(ncol,nlay+1), intent(inout) :: radn_dn !Top level must contain incident flux boundary condition + real(wp), dimension(ncol,nlay), intent(in ) :: An, Cn + ! Local variables + integer :: ilev, icol + ! --------------------------------------------------- + real(wp) :: adjustmentFactor + if(top_at_1) then + ! + ! Top of domain is index 1 + ! + ! 1st Upward propagation + do ilev = nlay, 1, -1 + radn_up(:,ilev) = trans(:,ilev )*radn_up(:,ilev+1) + source_up(:,ilev) + do icol=1,ncol + adjustmentFactor = Cn(icol,ilev)*( An(icol,ilev)*radn_dn(icol,ilev) - & + trans(icol,ilev)*source_dn(icol,ilev) - source_up(icol,ilev) ) + radn_up(icol,ilev) = radn_up(icol,ilev) + adjustmentFactor + enddo + end do + ! 2nd Downward propagation + do ilev = 1, nlay + radn_dn(:,ilev+1) = trans(:,ilev)*radn_dn(:,ilev) + source_dn(:,ilev) + do icol=1,ncol + adjustmentFactor = Cn(icol,ilev)*( An(icol,ilev)*radn_up(icol,ilev) - & + trans(icol,ilev)*source_up(icol,ilev) - source_dn(icol,ilev) ) + + radn_dn(:,ilev+1) = radn_dn(:,ilev+1) + adjustmentFactor + enddo + end do + else + ! + ! Top of domain is index nlay+1 + ! + ! Upward propagation + do ilev = 1, nlay + radn_up(:,ilev+1) = trans(:,ilev) * radn_up(:,ilev) + source_up(:,ilev) + do icol=1,ncol + adjustmentFactor = Cn(icol,ilev)*( An(icol,ilev)*radn_dn(icol,ilev+1) - & + trans(icol,ilev)*source_dn(icol,ilev) - source_up(icol,ilev) ) + radn_up(icol,ilev+1) = radn_up(icol,ilev+1) + adjustmentFactor + enddo + end do + + ! 2st Downward propagation + do ilev = nlay, 1, -1 + radn_dn(:,ilev) = trans(:,ilev )*radn_dn(:,ilev+1) + source_dn(:,ilev) + do icol=1,ncol + adjustmentFactor = Cn(icol,ilev)*( An(icol,ilev)*radn_up(icol,ilev) - & + trans(icol,ilev)*source_up(icol,ilev) - source_dn(icol,ilev) ) + radn_dn(icol,ilev) = radn_dn(icol,ilev) + adjustmentFactor + enddo + end do + end if +end subroutine lw_transport_1rescl end module mo_rte_solver_kernels diff --git a/rte/mo_rte_lw.F90 b/rte/mo_rte_lw.F90 index 468e1ec79..25dde3cad 100644 --- a/rte/mo_rte_lw.F90 +++ b/rte/mo_rte_lw.F90 @@ -42,7 +42,8 @@ module mo_rte_lw only: ty_source_func_lw use mo_fluxes, only: ty_fluxes use mo_rte_solver_kernels, & - only: apply_BC, lw_solver_noscat_GaussQuad, lw_solver_2stream + only: apply_BC, lw_solver_noscat_GaussQuad, lw_solver_2stream,& + lw_solver_1rescl_GaussQuad implicit none private @@ -56,7 +57,7 @@ module mo_rte_lw function rte_lw(optical_props, top_at_1, & sources, sfc_emis, & fluxes, & - inc_flux, n_gauss_angles) result(error_msg) + inc_flux, n_gauss_angles, use_2stream) result(error_msg) class(ty_optical_props_arry), intent(in ) :: optical_props ! Array of ty_optical_props. This type is abstract ! and needs to be made concrete, either as an array ! (class ty_optical_props_arry) or in some user-defined way @@ -67,10 +68,12 @@ function rte_lw(optical_props, top_at_1, & class(ty_fluxes), intent(inout) :: fluxes ! Array of ty_fluxes. Default computes broadband fluxes at all levels ! if output arrays are defined. Can be extended per user desires. real(wp), dimension(:,:), & - target, optional, intent(in ) :: inc_flux ! incident flux at domain top [W/m2] (ncol, ngpts) - integer, optional, intent(in ) :: n_gauss_angles ! Number of angles used in Gaussian quadrature - ! (no-scattering solution) - character(len=128) :: error_msg ! If empty, calculation was successful + target, optional, intent(in ) :: inc_flux ! incident flux at domain top [W/m2] (ncol, ngpts) + integer, optional, intent(in ) :: n_gauss_angles ! Number of angles used in Gaussian quadrature + ! (no-scattering solution) + logical, optional, intent(in ) :: use_2stream ! When 2-stream parameters (tau/ssa/g) are provided, use 2-stream methods + ! Default is to use re-scaled longwave transport + character(len=128) :: error_msg ! If empty, calculation was successful ! -------------------------------- ! ! Local variables @@ -80,6 +83,7 @@ function rte_lw(optical_props, top_at_1, & integer :: icol, iband, igpt real(wp), dimension(:,:,:), allocatable :: gpt_flux_up, gpt_flux_dn real(wp), dimension(:,:), allocatable :: sfc_emis_gpt + logical :: using_2stream ! -------------------------------------------------- ! ! Weights and angle secants for first order (k=1) Gaussian quadrature. @@ -160,6 +164,11 @@ function rte_lw(optical_props, top_at_1, & n_quad_angs = n_gauss_angles end if ! + ! Optionally - use 2-stream methods when low-order scattering properties are provided? + ! + using_2stream = .false. + if(present(use_2stream)) using_2stream = use_2stream + ! ! Ensure values of tau, ssa, and g are reasonable ! error_msg = optical_props%validate() @@ -213,18 +222,34 @@ function rte_lw(optical_props, top_at_1, & gpt_flux_up, gpt_flux_dn) !$acc exit data delete(optical_props%tau) class is (ty_optical_props_2str) - ! - ! two-stream calculation with scattering - ! - !$acc enter data copyin(optical_props%tau, optical_props%ssa, optical_props%g) - error_msg = optical_props%validate() - if(len_trim(error_msg) > 0) return - call lw_solver_2stream(ncol, nlay, ngpt, logical(top_at_1, wl), & - optical_props%tau, optical_props%ssa, optical_props%g, & - sources%lay_source, sources%lev_source_inc, sources%lev_source_dec, & - sfc_emis_gpt, sources%sfc_source, & - gpt_flux_up, gpt_flux_dn) - !$acc exit data delete(optical_props%tau, optical_props%ssa, optical_props%g) + if (using_2stream) then + ! + ! two-stream calculation with scattering + ! + !$acc enter data copyin(optical_props%tau, optical_props%ssa, optical_props%g) + error_msg = optical_props%validate() + if(len_trim(error_msg) > 0) return + call lw_solver_2stream(ncol, nlay, ngpt, logical(top_at_1, wl), & + optical_props%tau, optical_props%ssa, optical_props%g, & + sources%lay_source, sources%lev_source_inc, sources%lev_source_dec, & + sfc_emis_gpt, sources%sfc_source, & + gpt_flux_up, gpt_flux_dn) + !$acc exit data delete(optical_props%tau, optical_props%ssa, optical_props%g) + else + ! + ! Re-scaled solution to account for scattering + ! + !$acc enter data copyin(optical_props%tau, optical_props%ssa, optical_props%g) + call lw_solver_1rescl_GaussQuad(ncol, nlay, ngpt, logical(top_at_1, wl), & + n_quad_angs, gauss_Ds(1:n_quad_angs,n_quad_angs), & + gauss_wts(1:n_quad_angs,n_quad_angs), & + optical_props%tau, optical_props%ssa, optical_props%g, & + sources%lay_source, sources%lev_source_inc, & + sources%lev_source_dec, & + sfc_emis_gpt, sources%sfc_source,& + gpt_flux_up, gpt_flux_dn) + !$acc exit data delete(optical_props%tau, optical_props%ssa, optical_props%g) + endif class is (ty_optical_props_nstr) ! ! n-stream calculation