From 436f7d8108c7c8f5de08914fec3547d2ee8ee782 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Thu, 18 Apr 2024 17:01:27 -0400 Subject: [PATCH 01/11] New longwave integration angles and weights --- rte-frontend/mo_rte_lw.F90 | 25 ++++++++++++--------- rte-kernels/accel/mo_rte_solver_kernels.F90 | 12 +++++----- rte-kernels/mo_rte_solver_kernels.F90 | 12 +++++----- 3 files changed, 26 insertions(+), 23 deletions(-) diff --git a/rte-frontend/mo_rte_lw.F90 b/rte-frontend/mo_rte_lw.F90 index c36c55f8c..47a11ebe6 100644 --- a/rte-frontend/mo_rte_lw.F90 +++ b/rte-frontend/mo_rte_lw.F90 @@ -125,22 +125,25 @@ function rte_lw(optical_props, top_at_1, & real(wp), dimension(:,:), pointer :: inc_flux_diffuse ! -------------------------------------------------- ! - ! Weights and angle secants for first order (k=1) Gaussian quadrature. - ! Values from Table 2, Clough et al, 1992, doi:10.1029/92JD01419 - ! after Abramowitz & Stegun 1972, page 921 + ! Weights and angle secants for "Gauss-Jacobi-5" quadrature. + ! Values from Table 1, R. J. Hogan 2023, doi:10.1002/qj.4598 ! integer, parameter :: max_gauss_pts = 4 real(wp), parameter, & dimension(max_gauss_pts, max_gauss_pts) :: & - gauss_Ds = RESHAPE([1.66_wp, 0._wp, 0._wp, 0._wp, & ! Diffusivity angle, not Gaussian angle - 1.18350343_wp, 2.81649655_wp, 0._wp, 0._wp, & - 1.09719858_wp, 1.69338507_wp, 4.70941630_wp, 0._wp, & - 1.06056257_wp, 1.38282560_wp, 2.40148179_wp, 7.15513024_wp], & + ! + ! Values provided are for mu = cos(theta); we require the inverse + ! + gauss_Ds = 1._wp / & + RESHAPE([0.6096748751_wp, huge(1._wp) , huge(1._wp) , huge(1._wp), & + 0.2509907356_wp, 0.7908473988_wp, huge(1._wp) , huge(1._wp), & + 0.1024922169_wp, 0.4417960320_wp, 0.8633751621_wp, huge(1._wp), & + 0.0454586727_wp, 0.2322334416_wp, 0.5740198775_wp, 0.9030775973_wp], & [max_gauss_pts, max_gauss_pts]), & - gauss_wts = RESHAPE([0.5_wp, 0._wp, 0._wp, 0._wp, & - 0.3180413817_wp, 0.1819586183_wp, 0._wp, 0._wp, & - 0.2009319137_wp, 0.2292411064_wp, 0.0698269799_wp, 0._wp, & - 0.1355069134_wp, 0.2034645680_wp, 0.1298475476_wp, 0.0311809710_wp], & + gauss_wts = RESHAPE([1._wp, 0._wp, 0._wp, 0._wp, & + 0.2300253764_wp, 0.7699746236_wp, 0._wp, 0._wp, & + 0.0437820218_wp, 0.3875796738_wp, 0.5686383044_wp, 0._wp, & + 0.0092068785_wp, 0.1285704278_wp, 0.4323381850_wp, 0.4298845087_wp], & [max_gauss_pts, max_gauss_pts]) ! ------------------------------------------------------------------------------------ ncol = optical_props%get_ncol() diff --git a/rte-kernels/accel/mo_rte_solver_kernels.F90 b/rte-kernels/accel/mo_rte_solver_kernels.F90 index 42d234abd..36d2366ec 100644 --- a/rte-kernels/accel/mo_rte_solver_kernels.F90 +++ b/rte-kernels/accel/mo_rte_solver_kernels.F90 @@ -133,7 +133,7 @@ subroutine lw_solver_noscat_oneangle(ncol, nlay, ngpt, top_at_1, D, weight, & ! Transport is for intensity ! convert flux at top of domain to intensity assuming azimuthal isotropy ! - flux_dn(icol,top_level,igpt) = incident_flux(icol,igpt)/(2._wp * pi * weight) + flux_dn(icol,top_level,igpt) = incident_flux(icol,igpt)/(pi * weight) end do end do @@ -219,16 +219,16 @@ subroutine lw_solver_noscat_oneangle(ncol, nlay, ngpt, top_at_1, D, weight, & ! Broadband reduction including ! conversion from intensity to flux assuming azimuthal isotropy and quadrature weight ! - call sum_broadband_factor(ncol, nlay+1, ngpt, 2._wp * pi * weight, flux_dn, broadband_dn) - call sum_broadband_factor(ncol, nlay+1, ngpt, 2._wp * pi * weight, flux_up, broadband_up) + call sum_broadband_factor(ncol, nlay+1, ngpt, pi * weight, flux_dn, broadband_dn) + call sum_broadband_factor(ncol, nlay+1, ngpt, pi * weight, flux_up, broadband_up) !$acc exit data delete( flux_dn,flux_up) !$omp target exit data map(release:flux_dn,flux_up) else ! ! Convert intensity to flux assuming azimuthal isotropy and quadrature weight ! - call apply_factor_3D(ncol, nlay+1, ngpt, 2._wp*pi*weight, flux_dn) - call apply_factor_3D(ncol, nlay+1, ngpt, 2._wp*pi*weight, flux_up) + call apply_factor_3D(ncol, nlay+1, ngpt, pi * weight, flux_dn) + call apply_factor_3D(ncol, nlay+1, ngpt, pi * weight, flux_up) !$acc exit data copyout( flux_dn,flux_up) !$omp target exit data map(from:flux_dn,flux_up) end if @@ -236,7 +236,7 @@ subroutine lw_solver_noscat_oneangle(ncol, nlay, ngpt, top_at_1, D, weight, & ! Only broadband-integrated Jacobians are provided ! if (do_Jacobians) then - call sum_broadband_factor(ncol, nlay+1, ngpt, 2._wp * pi * weight, gpt_Jac, flux_upJac) + call sum_broadband_factor(ncol, nlay+1, ngpt, pi * weight, gpt_Jac, flux_upJac) end if !$acc end data diff --git a/rte-kernels/mo_rte_solver_kernels.F90 b/rte-kernels/mo_rte_solver_kernels.F90 index e10a0cb36..6a1156224 100644 --- a/rte-kernels/mo_rte_solver_kernels.F90 +++ b/rte-kernels/mo_rte_solver_kernels.F90 @@ -141,7 +141,7 @@ subroutine lw_solver_noscat_oneangle(ncol, nlay, ngpt, top_at_1, D, weight, & ! Transport is for intensity ! convert flux at top of domain to intensity assuming azimuthal isotropy ! - gpt_flux_dn(:,top_level) = incident_flux(:,igpt)/(2._wp * pi * weight) + gpt_flux_dn(:,top_level) = incident_flux(:,igpt)/(pi * weight) ! ! Optical path and transmission, used in source function and transport calculations ! @@ -220,8 +220,8 @@ subroutine lw_solver_noscat_oneangle(ncol, nlay, ngpt, top_at_1, D, weight, & ! ! Convert intensity to flux assuming azimuthal isotropy and quadrature weight ! - gpt_flux_dn(:,:) = 2._wp * pi * weight * gpt_flux_dn(:,:) - gpt_flux_up(:,:) = 2._wp * pi * weight * gpt_flux_up(:,:) + gpt_flux_dn(:,:) = pi * weight * gpt_flux_dn(:,:) + gpt_flux_up(:,:) = pi * weight * gpt_flux_up(:,:) end if ! ! Only broadband-integrated Jacobians are provided @@ -231,11 +231,11 @@ subroutine lw_solver_noscat_oneangle(ncol, nlay, ngpt, top_at_1, D, weight, & end do ! g point loop if(do_broadband) then - broadband_up(:,:) = 2._wp * pi * weight* broadband_up(:,:) - broadband_dn(:,:) = 2._wp * pi * weight* broadband_dn(:,:) + broadband_up(:,:) = pi * weight* broadband_up(:,:) + broadband_dn(:,:) = pi * weight* broadband_dn(:,:) end if if(do_Jacobians) & - flux_upJac(:,:) = 2._wp * pi * weight * flux_upJac(:,:) + flux_upJac(:,:) = pi * weight * flux_upJac(:,:) end subroutine lw_solver_noscat_oneangle ! ------------------------------------------------------------------------------------------------- From 38688a336128fc34336b2c92f509aace96cab6c1 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Thu, 18 Apr 2024 21:19:01 -0400 Subject: [PATCH 02/11] Layer source removed from solver kernels --- rte-frontend/mo_rte_lw.F90 | 8 +-- rte-kernels/accel/mo_rte_solver_kernels.F90 | 68 ++++++++++----------- rte-kernels/api/mo_rte_solver_kernels.F90 | 8 +-- rte-kernels/api/rte_kernels.h | 2 - rte-kernels/mo_rte_solver_kernels.F90 | 43 +++++-------- tests/rte_lw_solver_unit_tests.F90 | 4 +- 6 files changed, 52 insertions(+), 81 deletions(-) diff --git a/rte-frontend/mo_rte_lw.F90 b/rte-frontend/mo_rte_lw.F90 index 47a11ebe6..ca82c2fce 100644 --- a/rte-frontend/mo_rte_lw.F90 +++ b/rte-frontend/mo_rte_lw.F90 @@ -356,7 +356,6 @@ function rte_lw(optical_props, top_at_1, & logical(top_at_1, wl), n_quad_angs, & secants, gauss_wts(1:n_quad_angs,n_quad_angs), & optical_props%tau, & - sources%lay_source, & sources%lev_source, & sfc_emis_gpt, sources%sfc_source, & inc_flux_diffuse, & @@ -375,9 +374,9 @@ function rte_lw(optical_props, top_at_1, & ! 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, & - sfc_emis_gpt, sources%sfc_source, & - inc_flux_diffuse, & + sources%lev_source, & + sfc_emis_gpt, sources%sfc_source, & + inc_flux_diffuse, & gpt_flux_up, gpt_flux_dn) else allocate(secants(ncol, ngpt, n_quad_angs)) @@ -399,7 +398,6 @@ function rte_lw(optical_props, top_at_1, & logical(top_at_1, wl), n_quad_angs, & secants, gauss_wts(1:n_quad_angs,n_quad_angs), & optical_props%tau, & - sources%lay_source, & sources%lev_source, & sfc_emis_gpt, sources%sfc_source, & inc_flux_diffuse, & diff --git a/rte-kernels/accel/mo_rte_solver_kernels.F90 b/rte-kernels/accel/mo_rte_solver_kernels.F90 index 36d2366ec..cbf66616f 100644 --- a/rte-kernels/accel/mo_rte_solver_kernels.F90 +++ b/rte-kernels/accel/mo_rte_solver_kernels.F90 @@ -53,7 +53,7 @@ module mo_rte_solver_kernels ! ! --------------------------------------------------------------- subroutine lw_solver_noscat_oneangle(ncol, nlay, ngpt, top_at_1, D, weight, & - tau, lay_source, lev_source, sfc_emis, sfc_src, & + tau, lev_source, sfc_emis, sfc_src, & incident_flux, & flux_up, flux_dn, & do_broadband, broadband_up, broadband_dn, & @@ -64,7 +64,6 @@ subroutine lw_solver_noscat_oneangle(ncol, nlay, ngpt, top_at_1, D, weight, & real(wp), dimension(ncol, ngpt), intent(in ) :: D ! secant of propagation angle [] real(wp), intent(in ) :: weight ! quadrature weight real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: tau ! Absorption optical thickness [] - 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 @@ -174,10 +173,9 @@ subroutine lw_solver_noscat_oneangle(ncol, nlay, ngpt, top_at_1, D, weight, & trans (icol,ilay,igpt) = exp(-tau_loc(icol,ilay,igpt)) end if call lw_source_noscat(top_at_1, & - lay_source(icol,ilay,igpt), lev_source(icol,ilay,igpt), & - lev_source(icol,ilay+1,igpt), & - tau_loc (icol,ilay,igpt), trans (icol,ilay,igpt), & - source_dn (icol,ilay,igpt), source_up (icol,ilay,igpt)) + lev_source(icol,ilay,igpt), lev_source(icol,ilay+1,igpt),& + tau_loc (icol,ilay,igpt), trans (icol,ilay, igpt), & + source_dn (icol,ilay,igpt), source_up (icol,ilay, igpt)) end do end do end do @@ -253,13 +251,13 @@ end subroutine lw_solver_noscat_oneangle ! Routine sums over single-angle solutions for each sets of angles/weights ! ! --------------------------------------------------------------- - subroutine lw_solver_noscat(ncol, nlay, ngpt, top_at_1, & - nmus, Ds, weights, & - tau, & - lay_source, lev_source, & - sfc_emis, sfc_src, & - inc_flux, & - flux_up, flux_dn, & + subroutine lw_solver_noscat(ncol, nlay, ngpt, top_at_1, & + nmus, Ds, weights, & + tau, & + lev_source, & + sfc_emis, sfc_src, & + inc_flux, & + flux_up, flux_dn, & do_broadband, broadband_up, broadband_dn, & do_Jacobians, sfc_srcJac, flux_upJac, & do_rescaling, ssa, g) bind(C, name="rte_lw_solver_noscat") @@ -270,7 +268,6 @@ subroutine lw_solver_noscat(ncol, nlay, ngpt, top_at_1, & nmus), intent(in ) :: Ds real(wp), dimension(nmus), intent(in ) :: weights ! quadrature secants, weights real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: tau ! Absorption optical thickness [] - 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 ! Planck source at layer edge [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] @@ -301,8 +298,8 @@ subroutine lw_solver_noscat(ncol, nlay, ngpt, top_at_1, & integer :: icol, ilev, igpt, imu ! ------------------------------------ - !$acc data copyin(Ds, tau, lay_source, lev_source, sfc_emis, sfc_src) - !$omp target data map(to:Ds, tau, lay_source, lev_source, sfc_emis, sfc_src) + !$acc data copyin(Ds, tau, lev_source, sfc_emis, sfc_src) + !$omp target data map(to:Ds, tau, lev_source, sfc_emis, sfc_src) !$acc data copyout( flux_up, flux_dn) if (.not. do_broadband) !$omp target data map(from:flux_up, flux_dn) if (.not. do_broadband) !$acc data copyout( broadband_up, broadband_dn) if ( do_broadband) @@ -325,7 +322,7 @@ subroutine lw_solver_noscat(ncol, nlay, ngpt, top_at_1, & !$omp target data map(alloc:this_broadband_up, this_broadband_dn, this_flux_up, this_flux_dn) call lw_solver_noscat_oneangle(ncol, nlay, ngpt, & top_at_1, Ds(:,:,1), weights(1), tau, & - lay_source, lev_source, sfc_emis, sfc_src, & + lev_source, sfc_emis, sfc_src, & inc_flux, & this_flux_up, this_flux_dn, & do_broadband, this_broadband_up, this_broadband_dn, & @@ -359,7 +356,7 @@ subroutine lw_solver_noscat(ncol, nlay, ngpt, top_at_1, & do imu = 2, nmus call lw_solver_noscat_oneangle(ncol, nlay, ngpt, & top_at_1, Ds(:,:,imu), weights(imu), tau, & - lay_source, lev_source, sfc_emis, sfc_src, & + lev_source, sfc_emis, sfc_src, & inc_flux, & this_flux_up, this_flux_dn, & do_broadband, this_broadband_up, this_broadband_dn, & @@ -411,7 +408,7 @@ end subroutine lw_solver_noscat ! ------------------------------------------------------------------------------------------------- subroutine lw_solver_2stream (ncol, nlay, ngpt, top_at_1, & tau, ssa, g, & - lay_source, lev_source, sfc_emis, sfc_src, & + lev_source, sfc_emis, sfc_src, & inc_flux, & flux_up, flux_dn) bind(C, name="rte_lw_solver_2stream") integer, intent(in ) :: ncol, nlay, ngpt ! Number of columns, layers, g-points @@ -419,7 +416,6 @@ subroutine lw_solver_2stream (ncol, nlay, ngpt, top_at_1, & real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: tau, & ! Optical thickness, ssa, & ! single-scattering albedo, 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 ! Planck source at layer edge [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] @@ -433,8 +429,8 @@ subroutine lw_solver_2stream (ncol, nlay, ngpt, top_at_1, & real(wp), dimension(ncol ,ngpt) :: source_sfc ! ------------------------------------ ! ------------------------------------ - !$acc enter data copyin(tau, ssa, g, lay_source, lev_source, sfc_emis, sfc_src, flux_dn) - !$omp target enter data map(to:tau, ssa, g, lay_source, lev_source, sfc_emis, sfc_src, flux_dn) + !$acc enter data copyin(tau, ssa, g, lev_source, sfc_emis, sfc_src, flux_dn) + !$omp target enter data map(to:tau, ssa, g, lev_source, sfc_emis, sfc_src, flux_dn) !$acc enter data create( flux_up, Rdif, Tdif, gamma1, gamma2, sfc_albedo, source_dn, source_up, source_sfc) !$omp target enter data map(alloc:flux_up, Rdif, Tdif, gamma1, gamma2, sfc_albedo, source_dn, source_up, source_sfc) ! @@ -455,7 +451,7 @@ subroutine lw_solver_2stream (ncol, nlay, ngpt, top_at_1, & ! call lw_source_2str(ncol, nlay, ngpt, top_at_1, & sfc_emis, sfc_src, & - lay_source, lev_source, & + lev_source, & gamma1, gamma2, Rdif, Tdif, tau, & source_dn, source_up, source_sfc) @@ -475,8 +471,8 @@ subroutine lw_solver_2stream (ncol, nlay, ngpt, top_at_1, & Rdif, Tdif, & source_dn, source_up, source_sfc, & flux_up, flux_dn) - !$acc exit data delete( tau, ssa, g, lay_source, lev_source, sfc_emis, sfc_src) - !$omp target exit data map(release:tau, ssa, g, lay_source, lev_source, sfc_emis, sfc_src) + !$acc exit data delete( tau, ssa, g, lev_source, sfc_emis, sfc_src) + !$omp target exit data map(release:tau, ssa, g, lev_source, sfc_emis, sfc_src) !$acc exit data delete( Rdif, Tdif, gamma1, gamma2, sfc_albedo, source_dn, source_up, source_sfc) !$omp target exit data map(release:Rdif, Tdif, gamma1, gamma2, sfc_albedo, source_dn, source_up, source_sfc) !$acc exit data copyout(flux_up, flux_dn) @@ -694,14 +690,13 @@ end subroutine sw_solver_2stream ! This routine implements point-wise stencil, and has to be called in a loop ! ! --------------------------------------------------------------- - subroutine lw_source_noscat(top_at_1, lay_source, lev_source, levp1_source, tau, trans, & + subroutine lw_source_noscat(top_at_1, lev_source, levp1_source, tau, trans, & source_dn, source_up) !$acc routine seq !$omp declare target ! logical(wl), intent(in) :: top_at_1 - real(wp), intent(in) :: lay_source, & ! Planck source at layer center - lev_source, & ! Planck source at levels (layer edges), + real(wp), intent(in) :: lev_source, & ! Planck source at levels (layer edges), levp1_source, & ! Planck source at level +1 (layer edges), tau, & ! Optical path (tau/mu) trans ! Transmissivity (exp(-tau)) @@ -725,8 +720,8 @@ subroutine lw_source_noscat(top_at_1, lay_source, lev_source, levp1_source, tau, ! ! Equation below is developed in Clough et al., 1992, doi:10.1029/92JD01419, Eq 13 ! - source_inc = (1._wp - trans) * levp1_source + 2._wp * fact * (lay_source - levp1_source) - source_dec = (1._wp - trans) * lev_source + 2._wp * fact * (lay_source - lev_source) + source_inc = (1._wp - trans) * levp1_source + fact * (lev_source - levp1_source) + source_dec = (1._wp - trans) * lev_source + fact * (levp1_source - lev_source) if (top_at_1) then source_dn = source_inc source_up = source_dec @@ -923,14 +918,13 @@ end subroutine lw_two_stream ! --------------------------------------------------------------- subroutine lw_source_2str(ncol, nlay, ngpt, top_at_1, & sfc_emis, sfc_src, & - lay_source, lev_source, & + lev_source, & gamma1, gamma2, rdif, tdif, tau, source_dn, source_up, source_sfc) & bind (C, name="rte_lw_source_2str") integer, intent(in) :: ncol, nlay, ngpt logical(wl), intent(in) :: top_at_1 real(wp), dimension(ncol , ngpt), intent(in) :: sfc_emis, sfc_src - real(wp), dimension(ncol, nlay, ngpt), intent(in) :: lay_source, & ! Planck source at layer center - tau, & ! Optical depth (tau) + real(wp), dimension(ncol, nlay, ngpt), intent(in) :: tau, & ! Optical depth (tau) gamma1, gamma2,& ! Coupling coefficients rdif, tdif ! Layer reflectance and transmittance real(wp), dimension(ncol, nlay+1, ngpt), target, & @@ -943,8 +937,8 @@ subroutine lw_source_2str(ncol, nlay, ngpt, top_at_1, & real(wp) :: lev_source_bot, lev_source_top ! --------------------------------------------------------------- ! --------------------------------- - !$acc enter data copyin(sfc_emis, sfc_src, lay_source, tau, gamma1, gamma2, rdif, tdif, lev_source) - !$omp target enter data map(to:sfc_emis, sfc_src, lay_source, tau, gamma1, gamma2, rdif, tdif, lev_source) + !$acc enter data copyin(sfc_emis, sfc_src, tau, gamma1, gamma2, rdif, tdif, lev_source) + !$omp target enter data map(to:sfc_emis, sfc_src, tau, gamma1, gamma2, rdif, tdif, lev_source) !$acc enter data create(source_dn, source_up, source_sfc) !$omp target enter data map(alloc:source_dn, source_up, source_sfc) @@ -979,8 +973,8 @@ subroutine lw_source_2str(ncol, nlay, ngpt, top_at_1, & end do end do end do - !$acc exit data delete(sfc_emis, sfc_src, lay_source, tau, gamma1, gamma2, rdif, tdif, lev_source) - !$omp target exit data map(release:sfc_emis, sfc_src, lay_source, tau, gamma1, gamma2, rdif, tdif, lev_source) + !$acc exit data delete(sfc_emis, sfc_src, tau, gamma1, gamma2, rdif, tdif, lev_source) + !$omp target exit data map(release:sfc_emis, sfc_src, tau, gamma1, gamma2, rdif, tdif, lev_source) !$acc exit data copyout(source_dn, source_up, source_sfc) !$omp target exit data map(from:source_dn, source_up, source_sfc) diff --git a/rte-kernels/api/mo_rte_solver_kernels.F90 b/rte-kernels/api/mo_rte_solver_kernels.F90 index da8d0f706..8217b14b4 100644 --- a/rte-kernels/api/mo_rte_solver_kernels.F90 +++ b/rte-kernels/api/mo_rte_solver_kernels.F90 @@ -41,7 +41,7 @@ module mo_rte_solver_kernels subroutine lw_solver_noscat(ncol, nlay, ngpt, top_at_1, & nmus, Ds, weights, & tau, & - lay_source, lev_source, & + lev_source, & sfc_emis, sfc_src, & inc_flux, & flux_up, flux_dn, & @@ -62,8 +62,6 @@ subroutine lw_solver_noscat(ncol, nlay, ngpt, top_at_1, & !! quadrature weights real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: tau !! Absorption optical thickness [] - 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 !! Planck source at layer edge for radiation [W/m2] real(wp), dimension(ncol, ngpt), intent(in ) :: sfc_emis @@ -107,7 +105,7 @@ end subroutine lw_solver_noscat interface subroutine lw_solver_2stream (ncol, nlay, ngpt, top_at_1, & tau, ssa, g, & - lay_source, lev_source, sfc_emis, sfc_src, & + lev_source, sfc_emis, sfc_src, & inc_flux, & flux_up, flux_dn) bind(C, name="rte_lw_solver_2stream") use mo_rte_kind, only: wp, wl @@ -117,8 +115,6 @@ subroutine lw_solver_2stream (ncol, nlay, ngpt, top_at_1, & !! ilay = 1 is the top of the atmosphere? real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: tau, ssa, g !! Optical thickness, single-scattering albedo, 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 !! Planck source at layer edge for radiation [W/m2] real(wp), dimension(ncol, ngpt), intent(in ) :: sfc_emis diff --git a/rte-kernels/api/rte_kernels.h b/rte-kernels/api/rte_kernels.h index 71f86b81e..1f37a2094 100644 --- a/rte-kernels/api/rte_kernels.h +++ b/rte-kernels/api/rte_kernels.h @@ -60,7 +60,6 @@ extern "C" const Float* secants, // (nmus) const Float* weights, // (nmus) const Float* tau, // (ncol,nlay, ngpt) - const Float* lay_source, // (ncol,nlay, ngpt) const Float* lev_source, // (ncol,nlay+1,ngpt) const Float* sfc_emis, // (ncol, ngpt) const Float* sfc_src, // (ncol, ngpt) @@ -90,7 +89,6 @@ extern "C" const Float* tau, // (ncol,nlay, ngpt) const Float* ssa, // (ncol,nlay, ngpt) const Float* g, // (ncol,nlay, ngpt) - const Float* lay_source, // (ncol,nlay, ngpt) const Float* lev_source, // (ncol,nlay+1,ngpt) const Float* sfc_emis, // (ncol, ngpt) const Float* sfc_src, // (ncol, ngpt) diff --git a/rte-kernels/mo_rte_solver_kernels.F90 b/rte-kernels/mo_rte_solver_kernels.F90 index 6a1156224..d26d463a0 100644 --- a/rte-kernels/mo_rte_solver_kernels.F90 +++ b/rte-kernels/mo_rte_solver_kernels.F90 @@ -49,7 +49,7 @@ module mo_rte_solver_kernels ! ! --------------------------------------------------------------- subroutine lw_solver_noscat_oneangle(ncol, nlay, ngpt, top_at_1, D, weight, & - tau, lay_source, lev_source, sfc_emis, sfc_src, & + tau, lev_source, sfc_emis, sfc_src, & incident_flux, & flux_up, flux_dn, & do_broadband, broadband_up, broadband_dn, & @@ -60,7 +60,6 @@ subroutine lw_solver_noscat_oneangle(ncol, nlay, ngpt, top_at_1, D, weight, & real(wp), dimension(ncol, ngpt), intent(in ) :: D ! secant of propagation angle [] real(wp), intent(in ) :: weight ! quadrature weight real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: tau ! Absorption optical thickness [] - 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 ! Planck source at layer edge [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] @@ -186,7 +185,7 @@ subroutine lw_solver_noscat_oneangle(ncol, nlay, ngpt, top_at_1, D, weight, & ! Source function for diffuse radiation ! call lw_source_noscat(ncol, nlay, top_at_1, & - lay_source(:,:,igpt), lev_source(:,:,igpt), & + lev_source(:,:,igpt), & tau_loc, trans, source_dn, source_up) ! ! Transport down @@ -248,7 +247,7 @@ end subroutine lw_solver_noscat_oneangle subroutine lw_solver_noscat(ncol, nlay, ngpt, top_at_1, & nmus, Ds, weights, & tau, & - lay_source, lev_source, & + lev_source, & sfc_emis, sfc_src, & inc_flux, & flux_up, flux_dn, & @@ -268,8 +267,6 @@ subroutine lw_solver_noscat(ncol, nlay, ngpt, top_at_1, & !! quadrature weights real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: tau !! Absorption optical thickness [] - 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 !! Planck source at layer edge for radiation[W/m2] real(wp), dimension(ncol, ngpt), intent(in ) :: sfc_emis @@ -313,7 +310,7 @@ subroutine lw_solver_noscat(ncol, nlay, ngpt, top_at_1, & ! call lw_solver_noscat_oneangle(ncol, nlay, ngpt, & top_at_1, Ds(:,:,1), weights(1), tau, & - lay_source, lev_source, sfc_emis, sfc_src, & + lev_source, sfc_emis, sfc_src, & inc_flux, & flux_up, flux_dn, & do_broadband, broadband_up, broadband_dn, & @@ -343,7 +340,7 @@ subroutine lw_solver_noscat(ncol, nlay, ngpt, top_at_1, & do imu = 2, nmus call lw_solver_noscat_oneangle(ncol, nlay, ngpt, & top_at_1, Ds(:,:,imu), weights(imu), tau, & - lay_source, lev_source, sfc_emis, sfc_src, & + lev_source, sfc_emis, sfc_src, & inc_flux, & this_flux_up, this_flux_dn, & do_broadband, this_broadband_up, this_broadband_dn, & @@ -376,7 +373,7 @@ end subroutine lw_solver_noscat ! ------------------------------------------------------------------------------------------------- subroutine lw_solver_2stream (ncol, nlay, ngpt, top_at_1, & tau, ssa, g, & - lay_source, lev_source, sfc_emis, sfc_src, & + lev_source, sfc_emis, sfc_src, & inc_flux, & flux_up, flux_dn) bind(C, name="rte_lw_solver_2stream") integer, intent(in ) :: ncol, nlay, ngpt @@ -385,8 +382,6 @@ subroutine lw_solver_2stream (ncol, nlay, ngpt, top_at_1, & !! ilay = 1 is the top of the atmosphere? real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: tau, ssa, g !! Optical thickness, single-scattering albedo, 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 !! Planck source at layer edge temperature [W/m2] real(wp), dimension(ncol, ngpt), intent(in ) :: sfc_emis @@ -419,7 +414,7 @@ subroutine lw_solver_2stream (ncol, nlay, ngpt, top_at_1, & ! call lw_source_2str(ncol, nlay, top_at_1, & sfc_emis(:,igpt), sfc_src(:,igpt), & - lay_source(:,:,igpt), lev_source, & + lev_source, & gamma1, gamma2, Rdif, Tdif, tau(:,:,igpt), & source_dn, source_up, source_sfc) ! @@ -617,12 +612,11 @@ end subroutine sw_solver_2stream ! See Clough et al., 1992, doi: 10.1029/92JD01419, Eq 13 ! ! --------------------------------------------------------------- - subroutine lw_source_noscat(ncol, nlay, top_at_1, lay_source, lev_source, tau, trans, & + subroutine lw_source_noscat(ncol, nlay, top_at_1, lev_source, tau, trans, & source_dn, source_up) integer, intent(in) :: ncol, nlay logical(wl), intent(in) :: top_at_1 - real(wp), dimension(ncol, nlay ), intent(in) :: lay_source, & ! Planck source at layer center - tau, & ! Optical path (tau/mu) + real(wp), dimension(ncol, nlay ), intent(in) :: tau, & ! Optical path (tau/mu) trans ! Transmissivity (exp(-tau)) real(wp), dimension(ncol, nlay+1), intent(in) :: lev_source ! Planck source at levels (layer edges) real(wp), dimension(ncol, nlay ), target, & @@ -656,20 +650,12 @@ subroutine lw_source_noscat(ncol, nlay, top_at_1, lay_source, lev_source, tau, t end if ! ! Equation below is developed in Clough et al., 1992, doi:10.1029/92JD01419, Eq 13 + ! Adapted ! source_inc(icol,ilay) = (1._wp - trans(icol,ilay)) * lev_source(icol,ilay+1) + & - 2._wp * fact * (lay_source(icol,ilay) - lev_source(icol,ilay+1)) + fact * (lev_source(icol,ilay ) - lev_source(icol,ilay+1)) source_dec(icol,ilay) = (1._wp - trans(icol,ilay)) * lev_source(icol,ilay ) + & - 2._wp * fact * (lay_source(icol,ilay) - lev_source(icol,ilay )) - ! - ! Even better - omit the layer Planck source (not working so well) - ! - if(.false.) then - source_inc(icol,ilay) = (1._wp - trans(icol,ilay)) * lev_source(icol,ilay+1) + & - fact * (lev_source(icol,ilay ) - lev_source(icol,ilay+1)) - source_dec(icol,ilay) = (1._wp - trans(icol,ilay)) * lev_source(icol,ilay ) + & - fact * (lev_source(icol,ilay+1) - lev_source(icol,ilay )) - end if + fact * (lev_source(icol,ilay+1) - lev_source(icol,ilay )) end do end do end subroutine lw_source_noscat @@ -916,14 +902,13 @@ end subroutine lw_two_stream ! --------------------------------------------------------------- subroutine lw_source_2str(ncol, nlay, top_at_1, & sfc_emis, sfc_src, & - lay_source, lev_source, & + lev_source, & gamma1, gamma2, rdif, tdif, tau, source_dn, source_up, source_sfc) & bind (C, name="rte_lw_source_2str") integer, intent(in) :: ncol, nlay logical(wl), intent(in) :: top_at_1 real(wp), dimension(ncol ), intent(in) :: sfc_emis, sfc_src - real(wp), dimension(ncol, nlay), intent(in) :: lay_source, & ! Planck source at layer center - tau, & ! Optical depth (tau) + real(wp), dimension(ncol, nlay), intent(in) :: tau, & ! Optical depth (tau) gamma1, gamma2,& ! Coupling coefficients rdif, tdif ! Layer reflectance and transmittance real(wp), dimension(ncol, nlay+1), target, & diff --git a/tests/rte_lw_solver_unit_tests.F90 b/tests/rte_lw_solver_unit_tests.F90 index fc18afc44..6f900b5b3 100644 --- a/tests/rte_lw_solver_unit_tests.F90 +++ b/tests/rte_lw_solver_unit_tests.F90 @@ -57,8 +57,8 @@ program rte_lw_solver_unit_tests ! ! Longwave tests - gray radiative equilibrium ! - real(wp), parameter :: sigma = 5.670374419e-8_wp, & ! Stefan-Boltzmann constant - D = 1.66_wp ! Diffusivity angle, from single-angle RRTMGP solver + real(wp), parameter :: sigma = 5.670374419e-8_wp, & ! Stefan-Boltzmann constant + D = 1._wp/0.6096748751_wp ! Hogan 2023 optimal angle real(wp), dimension( ncol), parameter :: sfc_t = [(285._wp, icol = 1,ncol/2), & (310._wp, icol = 1,ncol/2)] real(wp), dimension( ncol), parameter :: total_tau = [0.1_wp, 1._wp, 10._wp, 50._wp, & From 4a12bebce51c0a9997f30ba1a9dbc7f4037552d6 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Thu, 18 Apr 2024 21:27:44 -0400 Subject: [PATCH 03/11] Remove layer source from gas optics kernels --- rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 | 6 ++--- .../accel/mo_gas_optics_rrtmgp_kernels.F90 | 11 +++----- .../api/mo_gas_optics_rrtmgp_kernels.F90 | 3 +-- rrtmgp-kernels/api/rrtmgp_kernels.h | 1 - .../mo_gas_optics_rrtmgp_kernels.F90 | 25 +------------------ 5 files changed, 8 insertions(+), 38 deletions(-) diff --git a/rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 b/rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 index 906c7cccb..3ccfc49ba 100644 --- a/rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 +++ b/rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 @@ -869,10 +869,10 @@ function source(this, & ! Source function needs temperature at interfaces/levels and at layer centers ! Allocate small local array for tlev unconditionally ! - !$acc data copyin(sources) copyout( sources%lay_source, sources%lev_source) & + !$acc data copyin(sources) copyout( sources%lev_source) & !$acc copyout( sources%sfc_source, sources%sfc_source_Jac) & !$acc create(tlev_arr) - !$omp target data map(from:sources%lay_source, sources%lev_source) & + !$omp target data map(from:sources%lev_source) & !$omp map(from:sources%sfc_source, sources%sfc_source_Jac) & !$omp map(alloc:tlev_arr) @@ -922,7 +922,7 @@ function source(this, & fmajor, jeta, tropo, jtemp, jpress, & this%get_gpoint_bands(), this%get_band_lims_gpoint(), this%planck_frac, this%temp_ref_min,& this%totplnk_delta, this%totplnk, this%gpoint_flavor, & - sources%sfc_source, sources%lay_source, sources%lev_source, & + sources%sfc_source, sources%lev_source, & sources%sfc_source_Jac) !$acc end data !$omp end target data diff --git a/rrtmgp-kernels/accel/mo_gas_optics_rrtmgp_kernels.F90 b/rrtmgp-kernels/accel/mo_gas_optics_rrtmgp_kernels.F90 index 103850776..4243aec5a 100644 --- a/rrtmgp-kernels/accel/mo_gas_optics_rrtmgp_kernels.F90 +++ b/rrtmgp-kernels/accel/mo_gas_optics_rrtmgp_kernels.F90 @@ -573,7 +573,7 @@ subroutine compute_Planck_source( & fmajor, jeta, tropo, jtemp, jpress, & gpoint_bands, band_lims_gpt, & pfracin, temp_ref_min, totplnk_delta, totplnk, gpoint_flavor, & - sfc_src, lay_src, lev_src, sfc_source_Jac) bind(C, name="rrtmgp_compute_Planck_source") + sfc_src, lev_src, sfc_source_Jac) bind(C, name="rrtmgp_compute_Planck_source") integer, intent(in) :: ncol, nlay, nbnd, ngpt integer, intent(in) :: nflav, neta, npres, ntemp, nPlanckTemp real(wp), dimension(ncol,nlay ), intent(in) :: tlay @@ -594,7 +594,6 @@ subroutine compute_Planck_source( & integer, dimension(2,ngpt), intent(in) :: gpoint_flavor real(wp), dimension(ncol, ngpt), intent(out) :: sfc_src - real(wp), dimension(ncol,nlay, ngpt), intent(out) :: lay_src real(wp), dimension(ncol,nlay+1,ngpt), intent(out) :: lev_src real(wp), dimension(ncol, ngpt), intent(out) :: sfc_source_Jac ! ----------------- @@ -609,9 +608,9 @@ subroutine compute_Planck_source( & ! ----------------- !$acc data copyin( tlay,tlev,tsfc,fmajor,jeta,tropo,jtemp,jpress,gpoint_bands,pfracin,totplnk,gpoint_flavor) & - !$acc copyout( sfc_src,lay_src,lev_src,sfc_source_Jac) + !$acc copyout( sfc_src,lev_src,sfc_source_Jac) !$omp target data map( to:tlay,tlev,tsfc,fmajor,jeta,tropo,jtemp,jpress,gpoint_bands,pfracin,totplnk,gpoint_flavor) & - !$omp map(from: sfc_src,lay_src,lev_src,sfc_source_Jac) + !$omp map(from: sfc_src,lev_src,sfc_source_Jac) ! Calculation of fraction of band's Planck irradiance associated with each g-point !$acc parallel loop tile(128,2) @@ -630,10 +629,6 @@ subroutine compute_Planck_source( & interpolate3D(one, fmajor(:,:,:,icol,ilay,iflav), pfracin, & igpt, jeta(:,icol,ilay,iflav), jtemp(icol,ilay),jpress(icol,ilay)+itropo) - ! Compute layer source irradiance for g-point, equals band irradiance x fraction for g-point - planck_function_1 = interpolate1D(tlay(icol,ilay), temp_ref_min, totplnk_delta, totplnk(:,ibnd)) - lay_src (icol,ilay,igpt) = pfrac * planck_function_1 - ! Compute level source irradiance for g-point planck_function_1 = interpolate1D(tlev(icol,ilay), temp_ref_min, totplnk_delta, totplnk(:,ibnd)) if (ilay == 1) then diff --git a/rrtmgp-kernels/api/mo_gas_optics_rrtmgp_kernels.F90 b/rrtmgp-kernels/api/mo_gas_optics_rrtmgp_kernels.F90 index 9272bb16e..7301a92c5 100644 --- a/rrtmgp-kernels/api/mo_gas_optics_rrtmgp_kernels.F90 +++ b/rrtmgp-kernels/api/mo_gas_optics_rrtmgp_kernels.F90 @@ -207,7 +207,7 @@ subroutine compute_Planck_source( & fmajor, jeta, tropo, jtemp, jpress, & gpoint_bands, band_lims_gpt, & pfracin, temp_ref_min, totplnk_delta, totplnk, gpoint_flavor, & - sfc_src, lay_src, lev_src, sfc_source_Jac) bind(C, name="rrtmgp_compute_Planck_source") + sfc_src, lev_src, sfc_source_Jac) bind(C, name="rrtmgp_compute_Planck_source") use mo_rte_kind, only : wp, wl integer, intent(in) :: ncol, nlay, nbnd, ngpt !! input dimensions @@ -235,7 +235,6 @@ subroutine compute_Planck_source( & integer, dimension(2,ngpt), intent(in) :: gpoint_flavor !! major gas flavor (pair) by upper/lower, g-point real(wp), dimension(ncol, ngpt), intent(out) :: sfc_src !! Planck emssion from the surface - real(wp), dimension(ncol,nlay, ngpt), intent(out) :: lay_src !! Planck emssion from layer centers real(wp), dimension(ncol,nlay+1,ngpt), intent(out) :: lev_src !! Planck emission at layer boundaries real(wp), dimension(ncol, ngpt), intent(out) :: sfc_source_Jac !! Jacobian (derivative) of the surface Planck source with respect to surface temperature diff --git a/rrtmgp-kernels/api/rrtmgp_kernels.h b/rrtmgp-kernels/api/rrtmgp_kernels.h index 92c36957c..4abc54465 100644 --- a/rrtmgp-kernels/api/rrtmgp_kernels.h +++ b/rrtmgp-kernels/api/rrtmgp_kernels.h @@ -116,7 +116,6 @@ extern "C" const Float* totplnk, // (nPlanckTemp,nbnd) const int* gpoint_flavor, // (2,ngpt) Float* sfc_src, // [out] (ncol, ngpt) - Float* lay_src, // [out] (ncol,nlay, ngpt) Float* lev_src, // [out] (ncol,nlay+1,ngpt) Float* sfc_src_jac // [out] (ncol, ngpt) ); diff --git a/rrtmgp-kernels/mo_gas_optics_rrtmgp_kernels.F90 b/rrtmgp-kernels/mo_gas_optics_rrtmgp_kernels.F90 index 7a050fb84..33fd525e3 100644 --- a/rrtmgp-kernels/mo_gas_optics_rrtmgp_kernels.F90 +++ b/rrtmgp-kernels/mo_gas_optics_rrtmgp_kernels.F90 @@ -571,7 +571,7 @@ subroutine compute_Planck_source( & fmajor, jeta, tropo, jtemp, jpress, & gpoint_bands, band_lims_gpt, & pfracin, temp_ref_min, totplnk_delta, totplnk, gpoint_flavor, & - sfc_src, lay_src, lev_src, sfc_source_Jac) bind(C, name="rrtmgp_compute_Planck_source") + sfc_src, lev_src, sfc_source_Jac) bind(C, name="rrtmgp_compute_Planck_source") integer, intent(in) :: ncol, nlay, nbnd, ngpt !! input dimensions integer, intent(in) :: nflav, neta, npres, ntemp, nPlanckTemp @@ -598,7 +598,6 @@ subroutine compute_Planck_source( & integer, dimension(2,ngpt), intent(in) :: gpoint_flavor !! major gas flavor (pair) by upper/lower, g-point real(wp), dimension(ncol, ngpt), intent(out) :: sfc_src !! Planck emission from the surface - real(wp), dimension(ncol,nlay, ngpt), intent(out) :: lay_src !! Planck emission from layer centers real(wp), dimension(ncol,nlay+1,ngpt), intent(out) :: lev_src !! Planck emission from layer boundaries real(wp), dimension(ncol, ngpt), intent(out) :: sfc_source_Jac !! Jacobian (derivative) of the surface Planck source with respect to surface temperature @@ -652,28 +651,6 @@ subroutine compute_Planck_source( & end do end do !icol - do ilay = 1, nlay - do icol = 1, ncol - ! Compute layer source irradiance for g-point, equals band irradiance x fraction for g-point - planck_function(icol,ilay,1:nbnd) = interpolate1D(tlay(icol,ilay), temp_ref_min, totplnk_delta, totplnk) - end do - end do - - ! - ! Map to g-points - ! - do ibnd = 1, nbnd - gptS = band_lims_gpt(1, ibnd) - gptE = band_lims_gpt(2, ibnd) - do igpt = gptS, gptE - do ilay = 1, nlay - do icol = 1, ncol - lay_src(icol,ilay,igpt) = pfrac(icol,ilay,igpt) * planck_function(icol,ilay,ibnd) - end do - end do - end do - end do - ! compute level source irradiances for each g-point do icol = 1, ncol planck_function (icol, 1,1:nbnd) = interpolate1D(tlev(icol, 1),temp_ref_min, totplnk_delta, totplnk) From 7a3b52ed6796c8cd27bbf537021f28154df4523b Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Thu, 18 Apr 2024 21:35:52 -0400 Subject: [PATCH 04/11] Remove layer source from Fortran front end and testing utilities --- examples/all-sky/rrtmgp_allsky.F90 | 4 ++-- examples/rfmip-clear-sky/rrtmgp_rfmip_lw.F90 | 8 ++++---- rte-frontend/mo_source_functions.F90 | 11 +++-------- tests/mo_testing_utils.F90 | 1 - tests/rte_lw_solver_unit_tests.F90 | 6 ------ 5 files changed, 9 insertions(+), 21 deletions(-) diff --git a/examples/all-sky/rrtmgp_allsky.F90 b/examples/all-sky/rrtmgp_allsky.F90 index e25ec0bf3..d7629155c 100644 --- a/examples/all-sky/rrtmgp_allsky.F90 +++ b/examples/all-sky/rrtmgp_allsky.F90 @@ -310,9 +310,9 @@ program rte_rrtmgp_allsky ! ! Should we allocate these once, rather than once per loop? They're big. ! - !$acc data create( lw_sources, lw_sources%lay_source, lw_sources%lev_source) & + !$acc data create( lw_sources, lw_sources%lev_source) & !$acc create( lw_sources%sfc_source, lw_sources%sfc_source_Jac) - !$omp target data map(alloc: lw_sources%lay_source, lw_sources%lev_source) & + !$omp target data map(alloc: lw_sources%lev_source) & !$omp map(alloc: lw_sources%sfc_source, lw_sources%sfc_source_Jac) call stop_on_err(k_dist%gas_optics(p_lay, p_lev, & t_lay, t_sfc, & diff --git a/examples/rfmip-clear-sky/rrtmgp_rfmip_lw.F90 b/examples/rfmip-clear-sky/rrtmgp_rfmip_lw.F90 index 607e372ef..59c05e8ee 100644 --- a/examples/rfmip-clear-sky/rrtmgp_rfmip_lw.F90 +++ b/examples/rfmip-clear-sky/rrtmgp_rfmip_lw.F90 @@ -219,8 +219,8 @@ program rrtmgp_rfmip_lw !$omp target enter data map(alloc:sfc_emis_spec) !$acc enter data create(optical_props, optical_props%tau) !$omp target enter data map(alloc:optical_props%tau) - !$acc enter data create(source, source%lay_source, source%lev_source, source%sfc_source) - !$omp target enter data map(alloc:source%lay_source, source%lev_source, source%sfc_source) + !$acc enter data create(source, source%lev_source, source%sfc_source) + !$omp target enter data map(alloc: source%lev_source, source%sfc_source) ! -------------------------------------------------- ! ! Loop over blocks @@ -265,8 +265,8 @@ program rrtmgp_rfmip_lw !$omp target exit data map(release:sfc_emis_spec) !$acc exit data delete(optical_props%tau, optical_props) !$omp target exit data map(release:optical_props%tau) - !$acc exit data delete(source%lay_source, source%lev_source, source%sfc_source) - !$omp target exit data map(release:source%lay_source, source%lev_source, source%sfc_source) + !$acc exit data delete(source%lev_source, source%sfc_source) + !$omp target exit data map(release:source%lev_source, source%sfc_source) !$acc exit data delete(source) ! --------------------------------------------------m call unblock_and_write(trim(flxup_file), 'rlu', flux_up) diff --git a/rte-frontend/mo_source_functions.F90 b/rte-frontend/mo_source_functions.F90 index 7df75257b..6a749c98d 100644 --- a/rte-frontend/mo_source_functions.F90 +++ b/rte-frontend/mo_source_functions.F90 @@ -28,8 +28,6 @@ module mo_source_functions !> spectral mapping in each direction separately, and at the surface !> type, extends(ty_optical_props), public :: ty_source_func_lw - real(wp), allocatable, dimension(:,:,:) :: lay_source - !! Planck source at layer average temperature (ncol, nlay, ngpt) real(wp), allocatable, dimension(:,:,:) :: lev_source !! Planck source at layer edge (ncol, nlay+1, ngpt) real(wp), allocatable, dimension(:,: ) :: sfc_source @@ -99,11 +97,10 @@ function alloc_lw(this, ncol, nlay) result(err_message) if(allocated(this%sfc_source)) deallocate(this%sfc_source) if(allocated(this%sfc_source_Jac)) deallocate(this%sfc_source_Jac) - if(allocated(this%lay_source)) deallocate(this%lay_source) if(allocated(this%lev_source)) deallocate(this%lev_source) ngpt = this%get_ngpt() - allocate(this%sfc_source (ncol, ngpt), this%lay_source (ncol,nlay,ngpt), & + allocate(this%sfc_source (ncol, ngpt), & this%lev_source (ncol,nlay+1,ngpt), this%sfc_source_Jac(ncol, ngpt)) end function alloc_lw ! -------------------------------------------------------------- @@ -176,7 +173,6 @@ end function copy_and_alloc_sw subroutine finalize_lw(this) class(ty_source_func_lw), intent(inout) :: this - if(allocated(this%lay_source )) deallocate(this%lay_source) if(allocated(this%lev_source )) deallocate(this%lev_source) if(allocated(this%sfc_source )) deallocate(this%sfc_source) if(allocated(this%sfc_source_Jac)) deallocate(this%sfc_source_Jac) @@ -199,7 +195,7 @@ pure function get_ncol_lw(this) integer :: get_ncol_lw if(this%is_allocated()) then - get_ncol_lw = size(this%lay_source,1) + get_ncol_lw = size(this%lev_source,1) else get_ncol_lw = 0 end if @@ -210,7 +206,7 @@ pure function get_nlay_lw(this) integer :: get_nlay_lw if(this%is_allocated()) then - get_nlay_lw = size(this%lay_source,2) + get_nlay_lw = size(this%lev_source,2)-1 else get_nlay_lw = 0 end if @@ -254,7 +250,6 @@ function get_subset_range_lw(full, start, n, subset) result(err_message) if(err_message /= "") return subset%sfc_source (1:n, :) = full%sfc_source (start:start+n-1, :) subset%sfc_source_Jac(1:n, :) = full%sfc_source_Jac(start:start+n-1, :) - subset%lay_source (1:n,:,:) = full%lay_source (start:start+n-1,:,:) subset%lev_source (1:n,:,:) = full%lev_source (start:start+n-1,:,:) end function get_subset_range_lw ! ------------------------------------------------------------------------------------------ diff --git a/tests/mo_testing_utils.F90 b/tests/mo_testing_utils.F90 index 07adffa39..550ec54d9 100644 --- a/tests/mo_testing_utils.F90 +++ b/tests/mo_testing_utils.F90 @@ -269,7 +269,6 @@ subroutine vr(atmos, sources) if(present(sources)) then sources%lev_source(:,:,:) = sources%lev_source(:,nlay+1:1:-1,:) - sources%lay_source(:,:,:) = sources%lay_source(:,nlay :1:-1,:) end if end subroutine vr ! ---------------------------------------------------------------------------- diff --git a/tests/rte_lw_solver_unit_tests.F90 b/tests/rte_lw_solver_unit_tests.F90 index 6f900b5b3..b04b6a9d4 100644 --- a/tests/rte_lw_solver_unit_tests.F90 +++ b/tests/rte_lw_solver_unit_tests.F90 @@ -276,18 +276,12 @@ subroutine gray_rad_equil(sfc_t, total_tau, nlay, top_at_1, atmos, sources) do ilay = 2, nlay+1 sources%lev_source(:,ilay, 1) = 0.5_wp/pi * olr(:) * & (1._wp + D * sum(atmos%tau(:,:ilay-1,1),dim=2)) - ! - ! The source is linear in optical depth so layer source is average of edges - ! - sources%lay_source(:,ilay-1,1) = 0.5_wp * (sources%lev_source(:,ilay, 1) + & - sources%lev_source(:,ilay-1,1)) end do if (.not. top_at_1) then ! ! Reverse vertical ordering of source functions ! sources%lev_source(:,1:nlay+1,1) = sources%lev_source(:,nlay+1:1:-1,1) - sources%lay_source(:,1:nlay, 1) = sources%lay_source(:,nlay :1:-1,1) end if end subroutine gray_rad_equil ! ------------------------------------------------------------------------------------ From ae6b61da0ecc4fb63ecbce1b3568e5b27e7f28ba Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Fri, 19 Apr 2024 11:42:39 -0400 Subject: [PATCH 05/11] CI uses updated reference data --- .github/workflows/containerized-ci.yml | 3 ++- .github/workflows/continuous-integration.yml | 3 ++- .github/workflows/self-hosted-ci.yml | 1 + 3 files changed, 5 insertions(+), 2 deletions(-) diff --git a/.github/workflows/containerized-ci.yml b/.github/workflows/containerized-ci.yml index 2548e2940..1c1bead9e 100644 --- a/.github/workflows/containerized-ci.yml +++ b/.github/workflows/containerized-ci.yml @@ -57,7 +57,7 @@ jobs: RUN_CMD: # https://github.com/earth-system-radiation/rte-rrtmgp/issues/194 OMP_TARGET_OFFLOAD: DISABLED - FAILURE_THRESHOLD: 5.8e-2 # 7.e-4 + FAILURE_THRESHOLD: 7.e-4 steps: # @@ -72,6 +72,7 @@ jobs: with: repository: earth-system-radiation/rrtmgp-data path: rrtmgp-data + ref: feature-new-lw-quadrature # # Build libraries, examples and tests (expect success) # diff --git a/.github/workflows/continuous-integration.yml b/.github/workflows/continuous-integration.yml index 138994fd4..a38ec53e4 100644 --- a/.github/workflows/continuous-integration.yml +++ b/.github/workflows/continuous-integration.yml @@ -31,7 +31,7 @@ jobs: RRTMGP_ROOT: ${{ github.workspace }} RRTMGP_DATA: ${{ github.workspace }}/rrtmgp-data RUN_CMD: - FAILURE_THRESHOLD: 5.8e-2 # 7.e-4 + FAILURE_THRESHOLD: 7.e-4 steps: # # Relax failure thresholds for single precision @@ -52,6 +52,7 @@ jobs: with: repository: earth-system-radiation/rrtmgp-data path: rrtmgp-data + ref: feature-new-lw-quadrature # # Synchronize the package index # diff --git a/.github/workflows/self-hosted-ci.yml b/.github/workflows/self-hosted-ci.yml index 5ef3cd98b..856812467 100644 --- a/.github/workflows/self-hosted-ci.yml +++ b/.github/workflows/self-hosted-ci.yml @@ -62,6 +62,7 @@ jobs: with: repository: earth-system-radiation/rrtmgp-data path: rrtmgp-data + ref: feature-new-lw-quadrature # # Finalize build environment # From 0cd245a6d00e89950647999abdcbf743db26659c Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Fri, 19 Apr 2024 11:52:34 -0400 Subject: [PATCH 06/11] Gitlab CI uses updated reference data --- .gitlab/levante.yml | 5 +++-- .gitlab/lumi.yml | 5 +++-- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/.gitlab/levante.yml b/.gitlab/levante.yml index 4835a4383..3ca5df136 100644 --- a/.gitlab/levante.yml +++ b/.gitlab/levante.yml @@ -7,6 +7,7 @@ include: file: '.slurm-ci.yml' variables: + RRTMGP_DATA_TAG: feature-update-lw-quadrature SCHEDULER_PARAMETERS: >- --account=mh0287 --time=05:00 @@ -47,7 +48,7 @@ variables: .dp: variables: FPMODEL: DP - FAILURE_THRESHOLD: "5.8e-2" + FAILURE_THRESHOLD: "7.e-4" .sp: variables: @@ -81,7 +82,7 @@ variables: # # Check out data # - - git clone --depth 1 https://github.com/earth-system-radiation/rrtmgp-data.git "${RRTMGP_DATA}" + - git clone --branch ${RRTMGP_DATA_TAG}--depth 1 https://github.com/earth-system-radiation/rrtmgp-data.git "${RRTMGP_DATA}" # # Run examples and tests # diff --git a/.gitlab/lumi.yml b/.gitlab/lumi.yml index 245129d5c..27d74e6a2 100644 --- a/.gitlab/lumi.yml +++ b/.gitlab/lumi.yml @@ -7,6 +7,7 @@ default: - lumi variables: + RRTMGP_DATA_TAG: feature-update-lw-quadrature SCHEDULER_PARAMETERS: >- --account=project_465000454 --nodes=1 @@ -39,7 +40,7 @@ variables: .dp: variables: FPMODEL: DP - FAILURE_THRESHOLD: "5.8e-2" + FAILURE_THRESHOLD: "7.e-4" .sp: variables: @@ -96,7 +97,7 @@ setup-python: # # Check out data # - - git clone --depth 1 https://github.com/earth-system-radiation/rrtmgp-data.git "${RRTMGP_DATA}" + - git clone --branch ${RRTMGP_DATA_TAG}--depth 1 https://github.com/earth-system-radiation/rrtmgp-data.git "${RRTMGP_DATA}" # # Run examples and tests # From 5842bd27f1714b7071cd05e4d369de606a314168 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Fri, 19 Apr 2024 11:57:38 -0400 Subject: [PATCH 07/11] Update tolerance for checking rad. equil. (single precision was failing) --- tests/rte_lw_solver_unit_tests.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/rte_lw_solver_unit_tests.F90 b/tests/rte_lw_solver_unit_tests.F90 index b04b6a9d4..a5897edf9 100644 --- a/tests/rte_lw_solver_unit_tests.F90 +++ b/tests/rte_lw_solver_unit_tests.F90 @@ -316,7 +316,7 @@ function check_gray_rad_equil(sfc_T, lw_tau, top_at_1, up_flux, net_flux) ! if(.not. allclose(net_flux(:,:), & spread(net_flux(:,1), dim=2, ncopies=size(net_flux,2)), & - tol = 70._wp)) then + tol = 100._wp)) then call report_err("Net flux not constant with tau in gray radiative equilibrium") check_gray_rad_equil = .false. end if From 88801cac4a6b3b1debb529f5cd7f279fba562c4c Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Fri, 19 Apr 2024 12:19:11 -0400 Subject: [PATCH 08/11] Update tolerance again for checking rad. equil. (nvhpc was failing) --- tests/rte_lw_solver_unit_tests.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/rte_lw_solver_unit_tests.F90 b/tests/rte_lw_solver_unit_tests.F90 index a5897edf9..e662af4ed 100644 --- a/tests/rte_lw_solver_unit_tests.F90 +++ b/tests/rte_lw_solver_unit_tests.F90 @@ -316,7 +316,7 @@ function check_gray_rad_equil(sfc_T, lw_tau, top_at_1, up_flux, net_flux) ! if(.not. allclose(net_flux(:,:), & spread(net_flux(:,1), dim=2, ncopies=size(net_flux,2)), & - tol = 100._wp)) then + tol = 128._wp)) then call report_err("Net flux not constant with tau in gray radiative equilibrium") check_gray_rad_equil = .false. end if From f27cf85a000b30fa6d19f1eb4cb7173ed23553e4 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Fri, 19 Apr 2024 12:33:23 -0400 Subject: [PATCH 09/11] OK so what is the tolerance we need for RE? --- tests/rte_lw_solver_unit_tests.F90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/rte_lw_solver_unit_tests.F90 b/tests/rte_lw_solver_unit_tests.F90 index e662af4ed..31c30e442 100644 --- a/tests/rte_lw_solver_unit_tests.F90 +++ b/tests/rte_lw_solver_unit_tests.F90 @@ -317,6 +317,8 @@ function check_gray_rad_equil(sfc_T, lw_tau, top_at_1, up_flux, net_flux) if(.not. allclose(net_flux(:,:), & spread(net_flux(:,1), dim=2, ncopies=size(net_flux,2)), & tol = 128._wp)) then + print *, maxval((net_flux(:,:) - & + spread(net_flux(:,1), dim=2, ncopies=size(net_flux,2)))/spacing(net_flux(:,:))) call report_err("Net flux not constant with tau in gray radiative equilibrium") check_gray_rad_equil = .false. end if From 059b1e6fb5057375513ea1bded6abe0c7d4778a3 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Fri, 19 Apr 2024 12:42:24 -0400 Subject: [PATCH 10/11] OK so what is the tolerance we need for RE (2/n)? --- tests/rte_lw_solver_unit_tests.F90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/tests/rte_lw_solver_unit_tests.F90 b/tests/rte_lw_solver_unit_tests.F90 index 31c30e442..f2a834a2f 100644 --- a/tests/rte_lw_solver_unit_tests.F90 +++ b/tests/rte_lw_solver_unit_tests.F90 @@ -307,6 +307,8 @@ function check_gray_rad_equil(sfc_T, lw_tau, top_at_1, up_flux, net_flux) if(.not. allclose(up_flux(:,toa), & gray_rad_equil_olr(sfc_t, lw_tau), tol=4._wp)) then call report_err("OLR is not consistent with gray radiative equilibrium") + print *, maxval((up_flux(:,toa) - & + gray_rad_equil_olr(sfc_t, lw_tau))/spacing(gray_rad_equil_olr(sfc_t, lw_tau))) check_gray_rad_equil = .false. end if ! @@ -316,7 +318,7 @@ function check_gray_rad_equil(sfc_T, lw_tau, top_at_1, up_flux, net_flux) ! if(.not. allclose(net_flux(:,:), & spread(net_flux(:,1), dim=2, ncopies=size(net_flux,2)), & - tol = 128._wp)) then + tol = 100._wp)) then print *, maxval((net_flux(:,:) - & spread(net_flux(:,1), dim=2, ncopies=size(net_flux,2)))/spacing(net_flux(:,:))) call report_err("Net flux not constant with tau in gray radiative equilibrium") From 98cf0ee9f527057fa284acc0b2712b021e3ab845 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Fri, 19 Apr 2024 12:59:01 -0400 Subject: [PATCH 11/11] Revised thresholds for OLR (single precision nvfortran) --- tests/rte_lw_solver_unit_tests.F90 | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/tests/rte_lw_solver_unit_tests.F90 b/tests/rte_lw_solver_unit_tests.F90 index f2a834a2f..3f77aa606 100644 --- a/tests/rte_lw_solver_unit_tests.F90 +++ b/tests/rte_lw_solver_unit_tests.F90 @@ -305,10 +305,8 @@ function check_gray_rad_equil(sfc_T, lw_tau, top_at_1, up_flux, net_flux) ! Check top-of-atmosphere energy balance ! if(.not. allclose(up_flux(:,toa), & - gray_rad_equil_olr(sfc_t, lw_tau), tol=4._wp)) then + gray_rad_equil_olr(sfc_t, lw_tau), tol=8._wp)) then call report_err("OLR is not consistent with gray radiative equilibrium") - print *, maxval((up_flux(:,toa) - & - gray_rad_equil_olr(sfc_t, lw_tau))/spacing(gray_rad_equil_olr(sfc_t, lw_tau))) check_gray_rad_equil = .false. end if ! @@ -319,8 +317,6 @@ function check_gray_rad_equil(sfc_T, lw_tau, top_at_1, up_flux, net_flux) if(.not. allclose(net_flux(:,:), & spread(net_flux(:,1), dim=2, ncopies=size(net_flux,2)), & tol = 100._wp)) then - print *, maxval((net_flux(:,:) - & - spread(net_flux(:,1), dim=2, ncopies=size(net_flux,2)))/spacing(net_flux(:,:))) call report_err("Net flux not constant with tau in gray radiative equilibrium") check_gray_rad_equil = .false. end if