Refactored shortwave two-stream calculation; accumulate broadband flu…

…xes within solvers; reduced resolution k-distributions (#138) Efficiency related changes: Refactored short-wave computes direct-beam source function alongside two-stream direct and diffuse reflectivity and transmissivity. Updates to reduced-resolution k-distributions; these are now ready for scientific use. Added option to solver to accumulate broadband fluxes as they are computed, rather than in a separate reduce() operation. Other corrections and cleanups. Co-authored-by: Dmitry Alexeev <dalexeev@nvidia.com>
earth-system-radiation · Sep 22, 2021 · 3b5d157 · 3b5d157
1 parent 124eb32
commit 3b5d157
Show file tree

Hide file tree

Showing 29 changed files with 1,687 additions and 1,252 deletions.
diff --git a/.github/workflows/containerized-ci.yml b/.github/workflows/containerized-ci.yml
@@ -38,6 +38,7 @@ jobs:
         cd ${RRTMGP_ROOT}
         ${FC} --version
         make libs
+        make -C build separate-libs
     ############################################################################
     - name: Cache RFMIP files
       id: cache-rfmip-files

diff --git a/.github/workflows/continuous-integration.yml b/.github/workflows/continuous-integration.yml
@@ -101,6 +101,7 @@ jobs:
         cd ${RRTMGP_ROOT}
         ${FC} --version
         make libs
+        make -C build separate-libs
 
     ############################################################################
     # Set up Python and packages

diff --git a/.github/workflows/self-hosted-ci.yml b/.github/workflows/self-hosted-ci.yml
@@ -32,7 +32,7 @@ jobs:
     env:
       FCFLAGS: ${{ matrix.FCFLAGS }}
       RTE_KERNELS: ${{ matrix.RTE_KERNELS }}
-      RUN_CMD: "srun -C gpu -A pr55 -p cscsci -t 15:00"
+      RUN_CMD: "srun -C gpu -A d56 -p cscsci -t 15:00"
     steps:
     - name: Check out code
       uses: actions/checkout@v2
@@ -64,11 +64,16 @@ jobs:
         source compiler_modules
         export RRTMGP_ROOT=$PWD
         export FC=ftn
+        # Compiler needs more temporary space than normally available 
+        mkdir -p $PWD/tmp
+        export TMPDIR=$PWD/tmp
         make clean
         make libs
+        make -C build separate-libs
     - name: Run
       run: |
         set -e
+        export TMPDIR=$PWD/tmp
         source compiler_modules
         module load cray-python
         export RRTMGP_ROOT=$PWD

diff --git a/README.md b/README.md
@@ -24,7 +24,7 @@ Relative to commit `69d36c9` to `master` on Apr 20, 2020, the required arguments
 2. Set environment variables `RRTMGP_ROOT` to the top-level RTE+RRTMGP directory and `RTE_KERNELS` to `openacc` if you want the OpenACC/OpenMP kernels rather than the default.
 3. `make libs` in the top-level directory will make the RTE and RRTMGP libraries.
 4. The examples and testing codes use netCDF. Set the variables `NCHOME` and `NFHOME` to the roots of the C and Fortran netCDF installations, then `make tests` to build and run these. (A few files need to be downloaded for `examples/rfmaip-clear-sky`. The default is to download these with `wget` but a Python script is also available.)
-5. Evaluating the results of the tests requires `Python` with the `xarray` package and its depdencies installed. Comparisons can be made with `make check` in the top level directory.
+5. Evaluating the results of the tests requires `Python` with the `xarray` package and its dependencies installed. Comparisons can be made with `make check` in the top level directory.
 6. `make` invoked without a target in the top level attempts all three steps.
 
 ## Examples

diff --git a/build/Makefile b/build/Makefile
@@ -14,8 +14,10 @@ RRTMGP_KERNEL_DIR += ../rrtmgp/kernels
 VPATH = $(RTE_DIR):$(RTE_KERNEL_DIR):$(RRTMGP_DIR):$(RRTMGP_KERNEL_DIR)
 #
 # Compiler variables FC, FCFLAGS must be set in the environment
-
-all: librte.a librrtmgp.a
+#
+# Make all the libraries though we'll only use the interface + kernels
+all:           librte.a                           librrtmgp.a
+separate-libs: librtekernels.a librteinterfacef.a librrtmgpkernels.a librrtmgpinterfacef.a
 
 COMPILE = $(FC) $(FCFLAGS) $(FCINCLUDE) -c
 %.o: %.F90
@@ -24,11 +26,23 @@ COMPILE = $(FC) $(FCFLAGS) $(FCINCLUDE) -c
 include $(RTE_DIR)/Make.depends
 include $(RRTMGP_DIR)/Make.depends
 
-librte.a: $(RTE_SRC)
-	ar -rvs librte.a $(RTE_SRC)
+librte.a: $(RTE_FORTRAN_KERNELS) $(RTE_FORTRAN_INTERFACE)
+	ar -rvs librte.a $(RTE_FORTRAN_KERNELS) $(RTE_FORTRAN_INTERFACE)
+
+librtekernels.a: $(RTE_FORTRAN_KERNELS)
+	ar -rvs librtekernels.a $(RTE_FORTRAN_KERNELS)
+
+librteinterfacef.a: $(RTE_FORTRAN_INTERFACE)
+	ar -rvs librteinterfacef.a $(RTE_FORTRAN_INTERFACE)
+
+librrtmgp.a: $(RRTMGP_FORTRAN_KERNELS)  $(RRTMGP_FORTRAN_INTERFACE)
+	ar -rvs librrtmgp.a $(RRTMGP_FORTRAN_KERNELS) $(RRTMGP_FORTRAN_INTERFACE)
+
+librrtmgpkernels.a: $(RRTMGP_FORTRAN_KERNELS)
+	ar -rvs librtekernels.a $(RRTMGP_FORTRAN_KERNELS)
 
-librrtmgp.a: $(RRTMGP_SRC)
-	ar -rvs librrtmgp.a $(RRTMGP_SRC)
+librrtmgpinterfacef.a: $(RRTMGP_FORTRAN_INTERFACE)
+	ar -rvs librteinterfacef.a $(RRTMGP_FORTRAN_INTERFACE)
 
 clean:
-	rm -f *.optrpt *.mod *.o librrtmgp.a librte.a
+	rm -f *.optrpt *.mod *.o lib*.a
diff --git a/environment.yml b/environment.yml
@@ -0,0 +1,13 @@
+#
+# Python modules below are needed to run tests and check results 
+#
+name: rte_rrtmgp
+dependencies:
+  - urllib3
+  - netcdf4
+  - xarray
+  - dask
+  - scipy
+  - matplotlib
+  - seaborn
+  - colorcet
diff --git a/examples/all-sky/README.md b/examples/all-sky/README.md
@@ -4,8 +4,4 @@ This example provides a modestly more realistic setting the clear-sky problem do
 
 The example uses the first of the Garand atmosphere used for developing RRTMGP, as described in the [paper](https://doi.org/10.1029/2019MS001621) documenting the code, repeats the column a user-specified number of times, computes the optical properties of an arbitrary cloud in each column, and computes the broadband fluxes. Fractional cloudiness is not considered, and the clouds are extensive but quite boring, with uniform condensate amount and particle size everywhere (though with different values for liquid and ice).
 
-1. Build the RTE+RRTMGP libraries in `../../build/`. This will require setting environmental variables `FC` for the Fortran compiler and `FCFLAGS`.
-2. Build the executables in this directory, which will require providing the locations of the netCDF C and Fortran libraries and module files as environmental
-variables `NCHOME` and `NFHOME`, as well a variable `RRTMGP_ROOT` pointing to the root of the installation (the absolute path to `../../`).
-4. Use Python script `run-rfmip-examples.py` to run the examples. The script takes some optional arguments, see `run-rfmip-examples.py -h`
-5. Python script `compare-to-reference.py` will compare the results to reference answers for 128 columns, produced on a Mac with Intel 19 Fortran compiler.
+Note that this example is run, and the results checked automatically, when `make` is invoked in the root directory. 
diff --git a/examples/all-sky/rrtmgp_allsky.F90 b/examples/all-sky/rrtmgp_allsky.F90
@@ -35,7 +35,7 @@ subroutine vmr_2d_to_1d(gas_concs, gas_concs_garand, name, sz1, sz2)
 end subroutine vmr_2d_to_1d
 ! ----------------------------------------------------------------------------------
 program rte_rrtmgp_clouds
-  use mo_rte_kind,           only: wp, i8
+  use mo_rte_kind,           only: wp, i8, wl
   use mo_optical_props,      only: ty_optical_props, &
                                    ty_optical_props_arry, ty_optical_props_1scl, ty_optical_props_2str
   use mo_gas_optics_rrtmgp,  only: ty_gas_optics_rrtmgp
@@ -49,6 +49,7 @@ program rte_rrtmgp_clouds
   use mo_load_cloud_coefficients, &
                              only: load_cld_lutcoeff, load_cld_padecoeff
   use mo_garand_atmos_io,    only: read_atmos, write_lw_fluxes, write_sw_fluxes
+  use mo_rte_config,         only: rte_config_checks
   implicit none
   ! ----------------------------------------------------------------------------------
   ! Variables
@@ -338,6 +339,9 @@ program rte_rrtmgp_clouds
   !
   !!$omp parallel do firstprivate(fluxes)
   do iloop = 1, nloops
+    ! Omit the checks starting with the second iteration
+    if (iloop == 2) call rte_config_checks(logical(.false., wl))
+
     call system_clock(start)
     call stop_on_err(                                      &
       cloud_optics%cloud_optics(lwp, iwp, rel, rei, clouds))
@@ -347,8 +351,10 @@ program rte_rrtmgp_clouds
     fluxes%flux_up => flux_up(:,:)
     fluxes%flux_dn => flux_dn(:,:)
     if(is_lw) then
-      !$acc enter data create(lw_sources, lw_sources%lay_source, lw_sources%lev_source_inc, lw_sources%lev_source_dec, lw_sources%sfc_source)
-      !$omp target enter data map(alloc:lw_sources%lay_source, lw_sources%lev_source_inc, lw_sources%lev_source_dec, lw_sources%sfc_source)
+      !$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)
       call stop_on_err(k_dist%gas_optics(p_lay, p_lev, &
                                          t_lay, t_sfc, &
                                          gas_concs,    &
@@ -360,10 +366,11 @@ program rte_rrtmgp_clouds
                               lw_sources,      &
                               emis_sfc,        &
                               fluxes))
-      !$acc exit data delete(lw_sources%lay_source, lw_sources%lev_source_inc, lw_sources%lev_source_dec, lw_sources%sfc_source, lw_sources)
-      !$omp target exit data map(release:lw_sources%lay_source, lw_sources%lev_source_inc, lw_sources%lev_source_dec, lw_sources%sfc_source)
+      !$acc        end data
+      !$omp end target data
+
     else
-      !$acc enter data create(toa_flux)
+      !$acc        enter data create(   toa_flux)
       !$omp target enter data map(alloc:toa_flux)
       fluxes%flux_dn_dir => flux_dir(:,:)
 
@@ -378,7 +385,7 @@ program rte_rrtmgp_clouds
                               mu0,   toa_flux, &
                               sfc_alb_dir, sfc_alb_dif, &
                               fluxes))
-      !$acc exit data delete(toa_flux)
+      !$acc        exit data delete(     toa_flux)
       !$omp target exit data map(release:toa_flux)
     end if
     !print *, "******************************************************************"

diff --git a/examples/mo_simple_netcdf.F90 b/examples/mo_simple_netcdf.F90
@@ -170,14 +170,17 @@ function write_2d_field(ncid, varName, var) result(err_msg)
     character(len=128)                   :: err_msg
 
     integer :: varid
+    integer :: stat
 
     err_msg = ""
     if(nf90_inq_varid(ncid, trim(varName), varid) /= NF90_NOERR) then
       err_msg = "write_field: can't find variable " // trim(varName)
       return
     end if
-    if(nf90_put_var(ncid, varid, var)  /= NF90_NOERR) &
-      err_msg = "write_field: can't write variable " // trim(varName)
+    stat = nf90_put_var(ncid, varid, var)
+    if(stat /= NF90_NOERR) &
+      err_msg = "write_field: can't write variable " // trim(varName) // &
+      " netcdf err: " // nf90_strerror(stat)
 
   end function write_2d_field
   !--------------------------------------------------------------------------------------------------------------------
@@ -188,14 +191,17 @@ function write_3d_field(ncid, varName, var) result(err_msg)
     character(len=128)                     :: err_msg
 
     integer :: varid
+    integer :: stat
 
     err_msg = ""
     if(nf90_inq_varid(ncid, trim(varName), varid) /= NF90_NOERR) then
       err_msg = "write_field: can't find variable " // trim(varName)
       return
     end if
-    if(nf90_put_var(ncid, varid, var)  /= NF90_NOERR) &
-      err_msg = "write_field: can't write variable " // trim(varName)
+    stat = nf90_put_var(ncid, varid, var)
+    if(stat /= NF90_NOERR) &
+      err_msg = "write_field: can't write variable " // trim(varName) // &
+      " netcdf err: " // nf90_strerror(stat)
 
   end function write_3d_field
   !--------------------------------------------------------------------------------------------------------------------
@@ -206,14 +212,17 @@ function write_4d_field(ncid, varName, var) result(err_msg)
     character(len=128)                     :: err_msg
 
     integer :: varid
+    integer :: stat
 
     err_msg = ""
     if(nf90_inq_varid(ncid, trim(varName), varid) /= NF90_NOERR) then
       err_msg = "write_field: can't find variable " // trim(varName)
       return
     end if
-    if(nf90_put_var(ncid, varid, var)  /= NF90_NOERR) &
-      err_msg = "write_field: can't write variable " // trim(varName)
+    stat = nf90_put_var(ncid, varid, var)
+    if(stat /= NF90_NOERR) &
+      err_msg = "write_field: can't write variable " // trim(varName) // &
+      " netcdf err: " // nf90_strerror(stat)
 
   end function write_4d_field
   !--------------------------------------------------------------------------------------------------------------------

diff --git a/examples/rfmip-clear-sky/Makefile b/examples/rfmip-clear-sky/Makefile
@@ -49,7 +49,7 @@ ADDITIONS = mo_simple_netcdf.o mo_rfmip_io.o mo_load_coefficients.o
 
 all: rrtmgp_rfmip_lw rrtmgp_rfmip_sw
 
-rrtmgp_rfmip_lw:   rrtmgp_rfmip_lw.o   $(ADDITIONS) $(RRTMGP_BUILD)/librte.a $(RRTMGP_BUILD)/librrtmgp.a
+rrtmgp_rfmip_lw:   rrtmgp_rfmip_lw.o   $(ADDITIONS) $(RRTMGP_BUILD)/librrtmgp.a $(RRTMGP_BUILD)/librte.a
 
 rrtmgp_rfmip_lw.o: rrtmgp_rfmip_lw.F90 $(ADDITIONS)
 

diff --git a/examples/rfmip-clear-sky/Readme.md b/examples/rfmip-clear-sky/Readme.md
@@ -3,17 +3,15 @@ This directory contains programs and support infrastructure for running
 the [RTE+RRTMGP](https://github.com/RobertPincus/rte-rrtmgp) radiation parameterization for the
 [RFMIP](https://www.earthsystemcog.org/projects/rfmip/) cases.
 
-1. Build the RTE+RRTMGP libraries in `../../build/`. This will require setting
-environmental variables `FC` for the Fortran compiler and `FCFLAGS`.
-2. Build the executables in this directory, which will require providing the
-locations of the netCDF C and Fortran libraries and module files as environmental
-variables `NCHOME` and `NFHOME`, as well a variable `RRTMGP_ROOT` pointing to the root of the installation
-(the absolute path to `../../`).
-3. Use Python script `stage_files.py` to download relevant files from the
+Note that this example is run, and the results checked automatically, when `make` is invoked in the root directory.
+
+A Python script `stage_files.py` is used to download relevant files from the
 [RFMIP web site](https://www.earthsystemcog.org/projects/rfmip/resources/).This script invokes another Python script to create empty output files.
-4. Use Python script `run-rfmip-examples.py` to run the examples. The script takes
+
+Use Python script `run-rfmip-examples.py` to run the examples. The script takes
 some optional arguments, see `run-rfmip-examples.py -h`
-5. Python script `compare-to-reference.py` will compare the results to reference
+
+Python script `compare-to-reference.py` will compare the results to reference
 answers produced on a Mac with Intel 19 Fortran compiler. Differences are normally
 within 10<sup>-6</sup> W/m<sup>2</sup>.
 

diff --git a/examples/rfmip-clear-sky/compare-to-reference.py b/examples/rfmip-clear-sky/compare-to-reference.py
@@ -15,7 +15,7 @@
 def construct_esgf_remote_name(var):
     #
     # For a given variable name, provide the OpenDAP URL for the RTE+RRTMGP RFMIP results
-    #   This doesn't seem to work on CSCS Piz Daint within the netcdf-python module 
+    #   This doesn't seem to work on CSCS Piz Daint within the netcdf-python module
     #
     esgf_url_base = "http://esgf3.dkrz.de/thredds/dodsC/cmip6/RFMIP/RTE-RRTMGP-Consortium/RTE-RRTMGP-181204/rad-irf/r1i1p1f1/Efx/"
     # DKRZ server has been unstable - better to try the other if one fails
@@ -61,9 +61,9 @@ def construct_esgf_file(var):
     failed = False
     for v in vars:
       if np.all(np.isnan(tst.variables[v].values)):
-        raise Exception(var + ": all test values are missing. Were the tests run?")
+        raise Exception(v + ": all test values are missing. Were the tests run?")
       if np.any(np.isnan(tst.variables[v].values)):
-        raise Exception(var + ": some test values are missing. Now that is strange.")
+        raise Exception(v + ": some test values are missing. Now that is strange.")
 
       diff = abs((tst-ref).variables[v].values)
       avg  = 0.5*(tst+ref).variables[v].values

diff --git a/examples/rfmip-clear-sky/mo_rfmip_io.F90 b/examples/rfmip-clear-sky/mo_rfmip_io.F90
@@ -205,25 +205,25 @@ subroutine determine_gas_names(concentrationFile, kdistFile, forcing_index, name
     character(len=32), dimension(11) :: &
       chem_name = ['co   ', &
                    'ch4  ', &
-        				   'o2   ', &
-        				   'n2o  ', &
-        				   'n2   ', &
-        				   'co2  ', &
-        				   'CCl4 ', &
-        				   'ch4  ', &
-        				   'CH3Br', &
-   			           'CH3Cl', &
+                   'o2   ', &
+                   'n2o  ', &
+                   'n2   ', &
+                   'co2  ', &
+                   'CCl4 ', &
+                   'ch4  ', &
+                   'CH3Br', &
+                   'CH3Cl', &
                    'cfc22'], &
       conc_name = ['carbon_monoxide     ', &
                    'methane             ', &
                    'oxygen              ', &
-          			   'nitrous_oxide       ', &
-          			   'nitrogen            ', &
-        				   'carbon_dioxide      ', &
-        				   'carbon_tetrachloride', &
-        				   'methane             ', &
-        				   'methyl_bromide      ', &
-        				   'methyl_chloride     ', &
+                   'nitrous_oxide       ', &
+                   'nitrogen            ', &
+                   'carbon_dioxide      ', &
+                   'carbon_tetrachloride', &
+                   'methane             ', &
+                   'methyl_bromide      ', &
+                   'methyl_chloride     ', &
                    'hcfc22              ']
     ! ----------------
     select case (forcing_index)

diff --git a/extensions/cloud_optics/mo_cloud_optics.F90 b/extensions/cloud_optics/mo_cloud_optics.F90
@@ -374,10 +374,10 @@ function cloud_optics(this, &
                         optical_props) result(error_msg)
     class(ty_cloud_optics), &
               intent(in   ) :: this
-    real(wp), intent(in   ) :: clwp  (:,:), &   ! cloud ice water path    (units?)
-                               ciwp  (:,:), &   ! cloud liquid water path (units?)
-                               reliq (:,:), &   ! cloud ice particle effective size (microns)
-                               reice (:,:)      ! cloud liquid particle effective radius (microns)
+    real(wp), intent(in   ) :: clwp  (:,:), &   ! cloud liquid water path (g/m2)
+                               ciwp  (:,:), &   ! cloud ice water path    (g/m2)
+                               reliq (:,:), &   ! cloud liquid particle effective size (microns)
+                               reice (:,:)      ! cloud ice particle effective radius  (microns)
     class(ty_optical_props_arry), &
               intent(inout) :: optical_props
                                                ! Dimensions: (ncol,nlay,nbnd)
@@ -558,7 +558,7 @@ function cloud_optics(this, &
       type is (ty_optical_props_nstr)
         error_msg = "cloud optics: n-stream calculations not yet supported"
       end select
-    end if 
+    end if
     !$acc end data
     !$omp end target data
   end function cloud_optics
@@ -643,7 +643,7 @@ subroutine compute_all_from_table(ncol, nlay, nbnd, mask, lwp, re, &
     integer  :: icol, ilay, ibnd
     integer  :: index
     real(wp) :: fint
-    real(wp) :: t, ts, tsg  ! tau, tau*ssa, tau*ssa*g
+    real(wp) :: t, ts  ! tau, tau*ssa, tau*ssa*g
     ! ---------------------------
     !$acc parallel loop gang vector default(present) collapse(3)
     !$omp target teams distribute parallel do simd collapse(3)
@@ -697,7 +697,7 @@ subroutine compute_all_from_pade(ncol, nlay, nbnd, nsizes, &
                                     intent(in) :: coeffs_asy
     real(wp), dimension(ncol,nlay,nbnd)        :: tau, taussa, taussag
     ! ---------------------------
-    integer  :: icol, ilay, ibnd, irad, count
+    integer  :: icol, ilay, ibnd, irad
     real(wp) :: t, ts
 
     !$acc parallel loop gang vector default(present) collapse(3)

diff --git a/rrtmgp/Make.depends b/rrtmgp/Make.depends
@@ -1,10 +1,12 @@
-RRTMGP_SRC = \
+RRTMGP_FORTRAN_KERNELS = \
+  mo_rrtmgp_util_reorder_kernels.o \
+  mo_gas_optics_kernels.o
+
+RRTMGP_FORTRAN_INTERFACE = \
   mo_rrtmgp_util_string.o \
 	mo_rrtmgp_constants.o \
-  mo_rrtmgp_util_reorder_kernels.o \
   mo_rrtmgp_util_reorder.o \
   mo_gas_concentrations.o \
-  mo_gas_optics_kernels.o \
   mo_gas_optics.o \
   mo_gas_optics_rrtmgp.o