Reduced-resolution k-distributions, refactored gas optics, misc impro…

…vements (#149)

Reduced-resolution k-distributions files halve computational cost with modest decreases in accuracy. Gas optics is refactored to avoid the transposition of large arrays. Small changes to documentation, CI, consistency in CPU and GPU implementations. 

OpenMP GPU implementation now fails on Cray CCE compiler; merging anyway.

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

@@ -6,6 +6,7 @@ jobs:
   CI:
     runs-on: daint
     strategy:
+      fail-fast: false
       matrix:
         include:
         - config_name: pgi_default_gpu
@@ -64,7 +65,7 @@ jobs:
         source compiler_modules
         export RRTMGP_ROOT=$PWD
         export FC=ftn
-        # Compiler needs more temporary space than normally available 
+        # Compiler needs more temporary space than normally available
         mkdir -p $PWD/tmp
         export TMPDIR=$PWD/tmp
         make clean

Compiler-flags.md

@@ -27,4 +27,4 @@ root of the RTE+RRTMGP installation.
 ### Debugging flags
 `FCFLAGS: "-g -Minfo -Mbounds -Mchkptr -Mstandard -Kieee -Mchkstk -Mallocatable=03  -Mpreprocess"`
 ### Optimization flags:  
-`FCFLAGS: "-m64 -O3 -g -traceback -heap-arrays -assume realloc_lhs -extend-source 132"`
+`FCFLAGS: "-g -O3 -fast -Minfo -Mallocatable=03 -Mpreprocess"`

examples/mo_load_coefficients.F90

@@ -35,6 +35,7 @@ subroutine stop_on_err(msg)
     use iso_fortran_env, only : error_unit
     character(len=*), intent(in) :: msg
 
+
     if(msg /= "") then
       write(error_unit, *) msg
       error stop 1

extensions/mo_rrtmgp_clr_all_sky.F90

@@ -176,12 +176,13 @@ function rte_lw(k_dist, gas_concs, p_lay, t_lay, p_lev,    &
                             inc_flux, n_gauss_angles)
 
     call sources%finalize()
+    call optical_props%finalize()
   end function rte_lw
   ! --------------------------------------------------
   function rte_sw(k_dist, gas_concs, p_lay, t_lay, p_lev, &
                                  mu0, sfc_alb_dir, sfc_alb_dif, cloud_props, &
                                  allsky_fluxes, clrsky_fluxes,           &
-                                 aer_props, col_dry, inc_flux, tsi_scaling) result(error_msg)
+                                 aer_props, col_dry, inc_flux) result(error_msg)
     class(ty_gas_optics),              intent(in   ) :: k_dist       !< derived type with spectral information
     type(ty_gas_concs),                intent(in   ) :: gas_concs    !< derived type encapsulating gas concentrations
     real(wp), dimension(:,:),          intent(in   ) :: p_lay, t_lay !< pressure [Pa], temperature [K] at layer centers (ncol,nlay)
@@ -198,7 +199,6 @@ function rte_sw(k_dist, gas_concs, p_lay, t_lay, p_lev, &
     real(wp), dimension(:,:), &
               optional,       intent(in ) :: col_dry, &  !< Molecular number density (ncol, nlay)
                                              inc_flux    !< incident flux at domain top [W/m2] (ncol, ngpts)
-    real(wp), optional,       intent(in ) :: tsi_scaling !< Optional scaling for total solar irradiance
 
     character(len=128)                    :: error_msg
     ! --------------------------------
@@ -227,9 +227,6 @@ function rte_sw(k_dist, gas_concs, p_lay, t_lay, p_lev, &
     ! ------------------------------------------------------------------------------------
     !  Error checking
     !
-    if(present(inc_flux) .and. present(tsi_scaling)) &
-      error_msg = "rrtmpg_sw: only one of inc_flux, tsi_scaling may be supplied."
-
     if(present(aer_props)) then
       if(any([aer_props%get_ncol(), &
               aer_props%get_nlay()] /= [ncol, nlay])) &
@@ -238,11 +235,6 @@ function rte_sw(k_dist, gas_concs, p_lay, t_lay, p_lev, &
         error_msg = "rrtmpg_sw: aerosol properties inconsistently sized"
     end if
 
-    if(present(tsi_scaling)) then
-      if(tsi_scaling <= 0._wp) &
-        error_msg = "rrtmpg_sw: tsi_scaling is < 0"
-    end if
-
     if(present(inc_flux)) then
       if(any([size(inc_flux, 1), &
               size(inc_flux, 2)] /= [ncol, ngpt])) &
@@ -291,7 +283,6 @@ function rte_sw(k_dist, gas_concs, p_lay, t_lay, p_lev, &
     ! If users have supplied an incident flux, use that
     !
     if(present(inc_flux))    toa_flux(:,:) = inc_flux(:,:)
-    if(present(tsi_scaling)) toa_flux(:,:) = toa_flux(:,:) * tsi_scaling
     ! ----------------------------------------------------
     ! Clear sky is gases + aerosols (if they're supplied)
     !
@@ -315,6 +306,7 @@ function rte_sw(k_dist, gas_concs, p_lay, t_lay, p_lev, &
                                sfc_alb_dir, sfc_alb_dif,        &
                                allsky_fluxes)
 
+    call optical_props%finalize()
   end function rte_sw
 
 end module mo_rrtmgp_clr_all_sky

tests/clear_sky_regression.F90

@@ -141,7 +141,7 @@ program rte_clear_sky_regression
   is_lw = .not. is_sw
   print *, "k-distribution is for the " // merge("longwave ", "shortwave", is_lw)
   print *, "  pressure    limits (Pa):", k_dist%get_press_min(), k_dist%get_press_max()
-  print *, "  temperature limits (Pa):", k_dist%get_temp_min(),  k_dist%get_temp_max()
+  print *, "  temperature limits (K):", k_dist%get_temp_min(),  k_dist%get_temp_max()
   !
   ! Problem sizes
   !
@@ -218,6 +218,7 @@ program rte_clear_sky_regression
       call load_and_init(k_dist_2, k_dist_file_2, gas_concs)
       print *, "Alternate k-distribution is for the " // merge("longwave ", "shortwave", .not. k_dist_2%source_is_external())
       print *, "  Resolution :", k_dist_2%get_nband(), k_dist_2%get_ngpt()
+      ngpt = k_dist_2%get_ngpt()
       call atmos%finalize()
       call make_optical_props_1scl(k_dist_2)
       call stop_on_err(lw_sources%alloc(ncol, nlay, k_dist_2))
@@ -566,6 +567,7 @@ end subroutine lw_clear_sky_incr
   subroutine lw_clear_sky_alt
     real(wp), dimension(ncol, nlay+1), target :: flux_net
     real(wp), dimension(ncol, nlay)           :: heating_rate
+    real(wp), dimension(ncol, ngpt)           :: lw_Ds
 
     fluxes%flux_net => flux_net
     call stop_on_err(k_dist_2%gas_optics(p_lay, p_lev, &
@@ -581,10 +583,21 @@ subroutine lw_clear_sky_alt
     call write_broadband_field(input_file, flux_up,  "lw_flux_up_alt",  "LW flux up, fewer g-points")
     call write_broadband_field(input_file, flux_dn,  "lw_flux_dn_alt",  "LW flux dn, fewer g-points")
     call write_broadband_field(input_file, flux_net, "lw_flux_net_alt", "LW flux ne, fewer g-pointst")
-
     call stop_on_err(compute_heating_rate(flux_up, flux_dn, p_lev, heating_rate))
     call write_broadband_field(input_file, heating_rate,  &
                                                      "lw_flux_hr_alt",  "LW heating rate, fewer g-points", vert_dim_name = "layer")
+
+    call stop_on_err(k_dist_2%compute_optimal_angles(atmos, lw_Ds))
+    call stop_on_err(rte_lw(atmos, top_at_1, &
+                            lw_sources,      &
+                            sfc_emis,        &
+                            fluxes, lw_Ds=lw_Ds))
+    call write_broadband_field(input_file, flux_up,  "lw_flux_up_alt_oa",  "LW flux up, fewer g-points, opt. angle")
+    call write_broadband_field(input_file, flux_dn,  "lw_flux_dn_alt_oa",  "LW flux dn, fewer g-points, opt. angle")
+    call write_broadband_field(input_file, flux_net, "lw_flux_net_alt_oa", "LW flux ne, fewer g-points, opt. angle")
+    call stop_on_err(compute_heating_rate(flux_up, flux_dn, p_lev, heating_rate))
+    call write_broadband_field(input_file, heating_rate,  &
+                                                     "lw_flux_hr_alt_oa",  "LW heating rate, fewer g-points, opt. angle", vert_dim_name = "layer")
     call k_dist_2%finalize()
   end subroutine lw_clear_sky_alt
   ! ----------------------------------------------------------------------------

tests/validation-plots.py

@@ -81,18 +81,20 @@ def construct_lbl_esgf_name(var, esgf_node="llnl"):
         #
         # Accuracy - 3-angle and single-angle
         #
-        variants = [[gpi.lw_flux_dn, gpi.lw_flux_dn_alt, gpi.lw_flux_dn_optang, gpi.lw_flux_dn_3ang, gpi.lw_flux_dn_2str],
-                    [gpi.lw_flux_up, gpi.lw_flux_up_alt, gpi.lw_flux_up_optang, gpi.lw_flux_up_3ang, gpi.lw_flux_up_2str],
+        variants = [[gpi.lw_flux_dn, gpi.lw_flux_dn_alt, gpi.lw_flux_dn_optang, gpi.lw_flux_dn_alt_oa, gpi.lw_flux_dn_3ang, gpi.lw_flux_dn_2str, gpi.lw_flux_dn_1rescl],
+                    [gpi.lw_flux_up, gpi.lw_flux_up_alt, gpi.lw_flux_up_optang, gpi.lw_flux_up_alt_oa, gpi.lw_flux_up_3ang, gpi.lw_flux_up_2str, gpi.lw_flux_up_1rescl],
                     [gpi.lw_flux_net,
                     gpi.lw_flux_net_alt,
                     gpi.lw_flux_dn_optang - gpi.lw_flux_up_optang,
+                    gpi.lw_flux_dn_alt_oa - gpi.lw_flux_up_alt_oa,
                     gpi.lw_flux_dn_3ang   - gpi.lw_flux_up_3ang,
-                    gpi.lw_flux_dn_2str   - gpi.lw_flux_up_2str]]
+                    gpi.lw_flux_dn_2str   - gpi.lw_flux_up_2str,
+                    gpi.lw_flux_dn_1rescl - gpi.lw_flux_up_1rescl]]
         refs = [lbli.rld,  lbli.rlu,  lbli.rld -  lbli.rlu]
         titles = ["Accuracy wrt LBLRTM: LW down", "Accuracy wrt LBLRTM: LW up", "Accuracy: LW net"]
         for v, r, t in zip(variants, refs, titles):
             make_comparison_plot(v, \
-                                 labels = ["default","fewer-g-points", "optimal-angle", "3-angle", "2-stream"], \
+                                 labels = ["default","fewer-g-points", "optimal-angle", "fewer points + optimal-angle", "3-angle", "2-stream", "rescaled"], \
                                  reference = r,            \
                                  vscale = plev/100.)
             plt.ylabel("Pressure (Pa)")