From 2516f688229b1e8fce7193cc9e0f54b11bd8b311 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Mon, 20 May 2024 17:53:13 -0400 Subject: [PATCH] Update LW quadrature angles (#282) Updates weights and secants for LW quadrature to the "Gauss-Jacobi-5" described in R. J. Hogan 2023, doi:10.1002/qj.4598. Changes answers for LW calculations. --- .github/workflows/containerized-ci.yml | 3 +- .github/workflows/continuous-integration.yml | 3 +- .github/workflows/self-hosted-ci.yml | 3 +- .gitlab/common.yml | 4 +- .gitlab/lumi.yml | 3 ++ examples/rfmip-clear-sky/stage_files.py | 45 -------------------- examples/rfmip-clear-sky/stage_files.sh | 5 --- rte-frontend/mo_rte_lw.F90 | 25 ++++++----- rte-kernels/accel/mo_rte_solver_kernels.F90 | 12 +++--- rte-kernels/mo_rte_solver_kernels.F90 | 12 +++--- tests/rte_lw_solver_unit_tests.F90 | 8 ++-- 11 files changed, 41 insertions(+), 82 deletions(-) delete mode 100755 examples/rfmip-clear-sky/stage_files.py delete mode 100644 examples/rfmip-clear-sky/stage_files.sh diff --git a/.github/workflows/containerized-ci.yml b/.github/workflows/containerized-ci.yml index 2548e2940..022fe2217 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: develop # # Build libraries, examples and tests (expect success) # diff --git a/.github/workflows/continuous-integration.yml b/.github/workflows/continuous-integration.yml index 138994fd4..9921917d6 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: develop # # Synchronize the package index # diff --git a/.github/workflows/self-hosted-ci.yml b/.github/workflows/self-hosted-ci.yml index 5ef3cd98b..3f344cd93 100644 --- a/.github/workflows/self-hosted-ci.yml +++ b/.github/workflows/self-hosted-ci.yml @@ -48,7 +48,7 @@ jobs: RRTMGP_DATA: ${{ github.workspace }}/rrtmgp-data RTE_KERNELS: ${{ matrix.rte-kernels }} RUN_CMD: "srun -C gpu -A d56 -p cscsci -t 15:00" - FAILURE_THRESHOLD: 5.8e-2 # 7.e-4 + FAILURE_THRESHOLD: 7.e-4 steps: # # Checks-out repository under $GITHUB_WORKSPACE @@ -62,6 +62,7 @@ jobs: with: repository: earth-system-radiation/rrtmgp-data path: rrtmgp-data + ref: develop # # Finalize build environment # diff --git a/.gitlab/common.yml b/.gitlab/common.yml index 196c641f4..914f3c381 100644 --- a/.gitlab/common.yml +++ b/.gitlab/common.yml @@ -1,7 +1,7 @@ .dp: variables: FPMODEL: DP - FAILURE_THRESHOLD: "5.8e-2" + FAILURE_THRESHOLD: "7.e-4" .sp: variables: @@ -16,7 +16,7 @@ RRTMGP_DATA: ${CI_PROJECT_DIR}/rrtmgp-data # Convenience variables: RRTMGP_DATA_REPO: https://github.com/earth-system-radiation/rrtmgp-data.git - RRTMGP_DATA_TAG: + RRTMGP_DATA_TAG: develop script: # # Build libraries, examples and tests diff --git a/.gitlab/lumi.yml b/.gitlab/lumi.yml index 1a642d992..9e2c1ad4c 100644 --- a/.gitlab/lumi.yml +++ b/.gitlab/lumi.yml @@ -8,6 +8,9 @@ include: .default: tags: - lumi + id_tokens: + CI_JOB_JWT: + aud: https://gitlab.com variables: SCHEDULER_PARAMETERS: >- diff --git a/examples/rfmip-clear-sky/stage_files.py b/examples/rfmip-clear-sky/stage_files.py deleted file mode 100755 index b23a6986e..000000000 --- a/examples/rfmip-clear-sky/stage_files.py +++ /dev/null @@ -1,45 +0,0 @@ -#! /usr/bin/env python -# -# This script downloads and/or creates files needed for the RFMIP off-line test -# cases -# -import sys - -import subprocess -import urllib.request -from pathlib import Path - -# -# Download and/or create input files and output template files -# -rte_rrtmgp_dir = Path("..").joinpath("..") -rfmip_dir = rte_rrtmgp_dir.joinpath("examples", "rfmip-clear-sky") -conds_file = "multiple_input4MIPs_radiation_RFMIP_UColorado-RFMIP-1-2_none.nc" -conds_url = ("http://aims3.llnl.gov/thredds/fileServer/user_pub_work/" - "input4MIPs/CMIP6/RFMIP/UColorado/UColorado-RFMIP-1-2/" + - "atmos/fx/multiple/none/v20190401/" + conds_file) -# -# The official RFMIP conditions are available from the ESFG, as above, but this -# fails from time to time, so we use the copy at Lamont-Doherty Earth -# Observatory, which we have to access via ftp(!) -# -conds_url = ("ftp://ftp.ldeo.columbia.edu/pub/robertp/rte-rrtmgp/" - "continuous-integration/" + conds_file) -output_files = [f"r{wl}{d}_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc" - for wl in ['l', 's'] for d in ['u', 'd']] -# -# Remove previous versions of files -# -for f in output_files: - Path(f).unlink(missing_ok=True) -Path(conds_file).unlink(missing_ok=True) - -# -# Download the profiles for RFMIP; download or make the empty output files -# -print("Downloading RFMIP input files") -urllib.request.urlretrieve(conds_url, conds_file) - -print("Downloading output templates") -for f in output_files: - urllib.request.urlretrieve(conds_url.replace(conds_file, f), f) diff --git a/examples/rfmip-clear-sky/stage_files.sh b/examples/rfmip-clear-sky/stage_files.sh deleted file mode 100644 index 7fa81db9f..000000000 --- a/examples/rfmip-clear-sky/stage_files.sh +++ /dev/null @@ -1,5 +0,0 @@ -wget ftp://ftp.ldeo.columbia.edu/pub/robertp/rte-rrtmgp/continuous-integration/multiple_input4MIPs_radiation_RFMIP_UColorado-RFMIP-1-2_none.nc -wget ftp://ftp.ldeo.columbia.edu/pub/robertp/rte-rrtmgp/continuous-integration/rld_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc -wget ftp://ftp.ldeo.columbia.edu/pub/robertp/rte-rrtmgp/continuous-integration/rlu_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc -wget ftp://ftp.ldeo.columbia.edu/pub/robertp/rte-rrtmgp/continuous-integration/rsd_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc -wget ftp://ftp.ldeo.columbia.edu/pub/robertp/rte-rrtmgp/continuous-integration/rsu_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc 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 ! ------------------------------------------------------------------------------------------------- diff --git a/tests/rte_lw_solver_unit_tests.F90 b/tests/rte_lw_solver_unit_tests.F90 index fc18afc44..baf08fd3f 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 ! Diffusivity angle, from single-angle RRTMGP solver 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, & @@ -311,7 +311,7 @@ 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") check_gray_rad_equil = .false. end if @@ -322,7 +322,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