Simplify LW source functions (#250)

Single level source for LW no-scattering solver; spectral mapping as the sqrt of the product of the mappings in each adjacent layer

.github/workflows/containerized-ci.yml

@@ -59,7 +59,7 @@ jobs:
       RUN_CMD:
       # https://github.com/earth-system-radiation/rte-rrtmgp/issues/194
       OMP_TARGET_OFFLOAD: DISABLED
-      FAILURE_THRESHOLD: 7.e-4
+      FAILURE_THRESHOLD: 5.8e-2 # 7.e-4
 
     steps:
     #

.github/workflows/continuous-integration.yml

@@ -29,7 +29,7 @@ jobs:
       RRTMGP_ROOT: ${{ github.workspace }}
       RRTMGP_DATA: ${{ github.workspace }}/rrtmgp-data
       RUN_CMD:
-      FAILURE_THRESHOLD: 7.e-4
+      FAILURE_THRESHOLD: 5.8e-2 # 7.e-4
     steps:
     #
     # Relax failure thresholds for single precision

.github/workflows/self-hosted-ci.yml

@@ -50,7 +50,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: 7.e-4
+      FAILURE_THRESHOLD: 5.8e-2 # 7.e-4
     steps:
     #
     # Checks-out repository under $GITHUB_WORKSPACE

examples/all-sky/rrtmgp_allsky.F90

@@ -310,10 +310,10 @@ 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_inc) &
-      !$acc             create(               lw_sources%lev_source_dec, lw_sources%sfc_source,     lw_sources%sfc_source_Jac)
-      !$omp target data map(alloc:            lw_sources%lay_source,     lw_sources%lev_source_inc) &
-      !$omp             map(alloc:            lw_sources%lev_source_dec, lw_sources%sfc_source,     lw_sources%sfc_source_Jac)
+      !$acc        data create(   lw_sources, lw_sources%lay_source,     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             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, &
                                          gas_concs,    &

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_inc, source%lev_source_dec, source%sfc_source)
-  !$omp target enter data map(alloc:source%lay_source, source%lev_source_inc, source%lev_source_dec, source%sfc_source)
+  !$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)
   ! --------------------------------------------------
   !
   ! 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_inc, source%lev_source_dec, source%sfc_source)
-  !$omp target exit data map(release:source%lay_source, source%lev_source_inc, source%lev_source_dec, source%sfc_source)
+  !$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)
   ! --------------------------------------------------m
   call unblock_and_write(trim(flxup_file), 'rlu', flux_up)

rrtmgp-frontend/mo_gas_optics_rrtmgp.F90

@@ -890,9 +890,9 @@ function source(this,                               &
     !-------------------------------------------------------------------
     ! Compute internal (Planck) source functions at layers and levels,
     !  which depend on mapping from spectral space that creates k-distribution.
-    !$acc        data copyin(sources) copyout( sources%lay_source, sources%lev_source_inc, sources%lev_source_dec) &
+    !$acc        data copyin(sources) copyout( sources%lay_source, sources%lev_source) &
     !$acc                             copyout( sources%sfc_source, sources%sfc_source_Jac)
-    !$omp target data                 map(from:sources%lay_source, sources%lev_source_inc, sources%lev_source_dec) &
+    !$omp target data                 map(from:sources%lay_source, sources%lev_source) &
     !$omp                             map(from:sources%sfc_source, sources%sfc_source_Jac)
 
     !$acc kernels copyout(top_at_1)
@@ -907,7 +907,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_inc, sources%lev_source_dec, &
+                sources%sfc_source, sources%lay_source, sources%lev_source, &
                 sources%sfc_source_Jac)
     !$acc end        data
     !$omp end target data

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_inc, lev_src_dec, sfc_source_Jac) bind(C, name="rrtmgp_compute_Planck_source")
+                    sfc_src, lay_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
@@ -593,25 +593,25 @@ subroutine compute_Planck_source(                        &
     real(wp), dimension(nPlanckTemp,nbnd),        intent(in) :: totplnk
     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,ngpt), intent(out) :: lev_src_inc, lev_src_dec
-    real(wp), dimension(ncol,     ngpt), intent(out) :: sfc_source_Jac
+    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
     ! -----------------
     ! local
     real(wp), parameter                             :: delta_Tsurf = 1.0_wp
 
     integer  :: ilay, icol, igpt, ibnd, itropo, iflav
     integer  :: gptS, gptE
     real(wp), dimension(2), parameter :: one = [1._wp, 1._wp]
-    real(wp) :: pfrac
+    real(wp) :: pfrac, pfrac_m1 ! Planck fraction in this layer and the one below 
     real(wp) :: planck_function_1, planck_function_2
     ! -----------------
 
     !$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_inc,lev_src_dec,sfc_source_Jac)
+    !$acc             copyout(  sfc_src,lay_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_inc,lev_src_dec,sfc_source_Jac)
+    !$omp             map(from: sfc_src,lay_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)
@@ -625,30 +625,38 @@ subroutine compute_Planck_source(                        &
           ! itropo = 1 lower atmosphere; itropo = 2 upper atmosphere
           itropo = merge(1,2,tropo(icol,ilay))  !WS moved itropo inside loop for GPU
           iflav = gpoint_flavor(itropo, igpt) !eta interpolation depends on band's flavor
+          ! interpolation in temperature, pressure, and eta
           pfrac = &
-            ! interpolation in temperature, pressure, and eta
             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 layer source irradiance for g-point, equals band irradiance x fraction for g-point
-          planck_function_1 = interpolate1D(tlev(icol,ilay),   temp_ref_min, totplnk_delta, totplnk(:,ibnd))
-          planck_function_2 = interpolate1D(tlev(icol,ilay+1), temp_ref_min, totplnk_delta, totplnk(:,ibnd))
-          lev_src_dec(icol,ilay,igpt) = pfrac * planck_function_1
-          lev_src_inc(icol,ilay,igpt) = pfrac * planck_function_2
+          lay_src  (icol,ilay,igpt) = pfrac * planck_function_1
 
+          ! Compute layer source irradiance for g-point, equals band irradiance x fraction for g-point
+          planck_function_1 = interpolate1D(tlev(icol,ilay), temp_ref_min, totplnk_delta, totplnk(:,ibnd))
+          if (ilay == 1) then 
+            lev_src(icol,ilay,  igpt) = pfrac * planck_function_1
+          else if (ilay == nlay) then 
+            lev_src(icol,ilay,  igpt) = pfrac * planck_function_1
+            planck_function_2 = interpolate1D(tlev(icol,nlay+1), temp_ref_min, totplnk_delta, totplnk(:,ibnd))
+            lev_src(icol,nlay+1,igpt) = pfrac * planck_function_2
+          else
+            pfrac_m1 = &
+              interpolate3D(one, fmajor(:,:,:,icol,ilay-1,iflav), pfracin, &
+                            igpt, jeta(:,icol,ilay-1,iflav), jtemp(icol,ilay-1),jpress(icol,ilay-1)+itropo)
+            lev_src(icol,ilay,  igpt) = sqrt(pfrac * pfrac_m1) * planck_function_1
+          end if 
           if (ilay == sfc_lay) then
             planck_function_1 = interpolate1D(tsfc(icol)              , temp_ref_min, totplnk_delta, totplnk(:,ibnd))
             planck_function_2 = interpolate1D(tsfc(icol) + delta_Tsurf, temp_ref_min, totplnk_delta, totplnk(:,ibnd))
             sfc_src       (icol,igpt) = pfrac * planck_function_1
             sfc_source_Jac(icol,igpt) = pfrac * (planck_function_2 - planck_function_1)
           end if
         end do ! igpt
-
       end do ! icol
     end do ! ilay
-
     !$acc end        data
     !$omp end target data
   end subroutine compute_Planck_source

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_inc, lev_src_dec, sfc_source_Jac) bind(C, name="rrtmgp_compute_Planck_source")
+                    sfc_src, lay_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
@@ -597,11 +597,10 @@ subroutine compute_Planck_source(                        &
     real(wp), dimension(nPlanckTemp,nbnd),        intent(in) :: totplnk       !! Total Planck function by band at each temperature 
     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,ngpt), intent(out) :: lev_src_inc, lev_src_dec
-      !! Planck emission at layer boundaries, using spectral mapping in the direction of propagation 
-    real(wp), dimension(ncol,     ngpt), intent(out) :: sfc_source_Jac 
+    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 
     ! -----------------
     ! local
@@ -675,11 +674,13 @@ subroutine compute_Planck_source(                        &
       end do
     end do
 
-    ! compute level source irradiances for each g-point, one each for upward and downward paths
+    ! 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)
+    end do
     do ilay = 1, nlay
       do icol = 1, ncol
-      planck_function(icol,     1,1:nbnd) = interpolate1D(tlev(icol,     1),temp_ref_min, totplnk_delta, totplnk)
-      planck_function(icol,ilay+1,1:nbnd) = interpolate1D(tlev(icol,ilay+1),temp_ref_min, totplnk_delta, totplnk)
+        planck_function(icol,ilay+1,1:nbnd) = interpolate1D(tlev(icol,ilay+1),temp_ref_min, totplnk_delta, totplnk)
       end do
     end do
 
@@ -690,12 +691,19 @@ subroutine compute_Planck_source(                        &
       gptS = band_lims_gpt(1, ibnd)
       gptE = band_lims_gpt(2, ibnd)
       do igpt = gptS, gptE
-        do ilay = 1, nlay
+        do icol = 1, ncol
+          lev_src(icol,     1,igpt) = pfrac(icol,   1,igpt) * planck_function(icol,     1,ibnd)
+        end do
+        do ilay = 2, nlay
           do icol = 1, ncol
-            lev_src_inc(icol,ilay,igpt) = pfrac(icol,ilay,igpt) *planck_function(icol,ilay+1,ibnd)
-            lev_src_dec(icol,ilay,igpt) = pfrac(icol,ilay,igpt) *planck_function(icol,ilay  ,ibnd)
+            lev_src(icol,ilay,igpt) = sqrt(pfrac(icol,ilay-1, igpt) *  & 
+                                           pfrac(icol,ilay,   igpt)) & 
+                                                            * planck_function(icol,ilay,  ibnd)
           end do
         end do
+        do icol = 1, ncol
+          lev_src(icol,nlay+1,igpt) = pfrac(icol,nlay,igpt) * planck_function(icol,nlay+1,ibnd)
+        end do
       end do
     end do
 

rte-frontend/mo_rte_lw.F90

@@ -353,8 +353,8 @@ 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_inc, &
-                                sources%lev_source_dec,            &
+                                sources%lay_source,                &
+                                sources%lev_source,                &
                                 sfc_emis_gpt, sources%sfc_source,  &
                                 inc_flux_diffuse,                  &
                                 gpt_flux_up, gpt_flux_dn,          &
@@ -371,8 +371,8 @@ function rte_lw(optical_props, top_at_1, &
             ! two-stream calculation with scattering
             !
             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, &
+                                   optical_props%tau, optical_props%ssa, optical_props%g, &
+                                   sources%lay_source, sources%lev_source,                &
                                    sfc_emis_gpt, sources%sfc_source,       &
                                    inc_flux_diffuse,                       &
                                    gpt_flux_up, gpt_flux_dn)
@@ -396,8 +396,8 @@ 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_inc, &
-                                  sources%lev_source_dec,            &
+                                  sources%lay_source,                &
+                                  sources%lev_source,                &
                                   sfc_emis_gpt, sources%sfc_source,  &
                                   inc_flux_diffuse,                  &
                                   gpt_flux_up, gpt_flux_dn,          &