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.
earth-system-radiation · May 20, 2024 · 2516f68 · 2516f68
1 parent da9fc8b
commit 2516f68
Show file tree

Hide file tree

Showing 11 changed files with 41 additions and 82 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
diff --git a/examples/rfmip-clear-sky/stage_files.sh b/examples/rfmip-clear-sky/stage_files.sh
diff --git 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
@@ -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,24 +219,24 @@ 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
     !
     ! 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
@@ -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
@@ -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