diff --git a/.gitignore b/.gitignore index 687e3f623..8e033b713 100644 --- a/.gitignore +++ b/.gitignore @@ -14,4 +14,11 @@ Makefile.libs Makefile.conf -src/ +# Intel code coverage +*.dyn +*.dpi +*.spi +*.spl +*.HTML +*.html +*.gif diff --git a/.travis.yml b/.travis.yml index eaba33de1..f2cc6ad31 100644 --- a/.travis.yml +++ b/.travis.yml @@ -7,7 +7,8 @@ language: cpp # - RTE_KERNELS="openacc" matrix: include: - - compiler: gcc # Using + - name: "PGI community edition " + compiler: gcc env: - FC=pgfortran - FCFLAGS="-Mallocatable=03 -Mstandard -Mbounds -Mchkptr -Kieee -Mchkstk" @@ -20,10 +21,11 @@ matrix: - gcc-8 - gfortran-8 - libnetcdf-dev - - compiler: gcc + - name: "gcc-9" + compiler: gcc env: - FC=gfortran-9 - - FCFLAGS="-ffree-line-length-none -m64 -std=f2003 -march=native -fbounds-check -finit-real=nan -DUSE_CBOOL" + - FCFLAGS="-ffree-line-length-none -m64 -std=f2008 -march=native -fbounds-check -finit-real=nan -DUSE_CBOOL" - TEST_PGI=no - NCHOME=/usr/ addons: @@ -34,10 +36,11 @@ matrix: - gcc-9 - gfortran-9 - libnetcdf-dev - - compiler: gcc + - name: "gcc-8" + compiler: gcc env: - FC=gfortran-8 - - FCFLAGS="-ffree-line-length-none -m64 -std=f2003 -march=native -fbounds-check -finit-real=nan -DUSE_CBOOL" + - FCFLAGS="-ffree-line-length-none -m64 -std=f2008 -march=native -fbounds-check -finit-real=nan -DUSE_CBOOL" - TEST_PGI=no - NCHOME=/usr/ addons: @@ -64,6 +67,7 @@ before_install: fi # Install netcdf-fortran for PGI +# git clone https://github.com/Unidata/netcdf-fortran.git - git clone https://github.com/Unidata/netcdf-fortran.git --branch v4.4.4 - cd netcdf-fortran - FC=$FC ./configure --prefix=${HOME}/netcdff @@ -82,57 +86,91 @@ install: # conda create --yes -n test python=$TRAVIS_PYTHON_VERSION - conda create --yes -n test python=3.7 - source activate test -- conda install --yes urllib3 netcdf4 xarray dask +- conda install --yes urllib3 netcdf4 xarray dask scipy before_script: - export RRTMGP_ROOT=$PWD -- export RRTMGP_DIR=$PWD/build +- export RRTMGP_BUILD=$PWD/build - export NFHOME=${HOME}/netcdff - export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$NFHOME/lib - cd examples/rfmip-clear-sky - python3 ./stage_files.py - cd ../.. script: -# Build library - default kernels -- cd build || exit 1 -- make || exit 1 -- cd .. || exit 1 +# +# Default kernels +# +- cd ${RRTMGP_ROOT} +# +# Build library +# +- make -C build clean || exit 1 +- make -C build || exit 1 +# # Build and run RFMIP examples -- cd examples/rfmip-clear-sky || exit 1 +# +- cd ${RRTMGP_ROOT}/examples/rfmip-clear-sky || exit 1 - make clean - make || exit 1 - python3 ./run-rfmip-examples.py -- python3 ./compare-to-reference.py --fail=1.e-4 -- cd ../.. +- python3 ./compare-to-reference.py --fail=7.e-4 +- make clean +# # Build and run all-sky example -- cd examples/all-sky || exit 1 +# +- cd ${RRTMGP_ROOT}/examples/all-sky || exit 1 - make clean - make || exit 1 - python3 ./run-allsky-example.py - python3 ./compare-to-reference.py -- cd ../.. +- make clean # -# Build library - Open-acc kernels +# Build and regression tests # -- export RTE_KERNELS="openacc" -- cd build || exit 1 +- cd ${RRTMGP_ROOT}/tests|| exit 1 - make clean - make || exit 1 -- cd .. || exit 1 +- cp ${RRTMGP_ROOT}/examples/rfmip-clear-sky/multiple_input4MIPs_radiation_RFMIP_UColorado-RFMIP-1-2_none.nc test_atmospheres.nc +- ./clear_sky_regression test_atmospheres.nc ${RRTMGP_ROOT}/rrtmgp/data/rrtmgp-data-lw-g256-2018-12-04.nc +- ./clear_sky_regression test_atmospheres.nc ${RRTMGP_ROOT}/rrtmgp/data/rrtmgp-data-sw-g224-2018-12-04.nc +- python3 verification.py +# +# Open-acc kernels +# +- cd ${RRTMGP_ROOT} +# +# Build library +# +- make -C build clean || exit 1 +- make -C build || exit 1 +# # Build and run RFMIP examples -- cd examples/rfmip-clear-sky || exit 1 +# +- cd ${RRTMGP_ROOT}/examples/rfmip-clear-sky || exit 1 - make clean - make || exit 1 - python3 ./run-rfmip-examples.py -- python3 ./compare-to-reference.py --fail=1.e-4 -- cd ../.. +- python3 ./compare-to-reference.py --fail=7.e-4 +- make clean +# # Build and run all-sky example -- cd examples/all-sky || exit 1 +# +- cd ${RRTMGP_ROOT}/examples/all-sky || exit 1 - make clean - make || exit 1 - python3 ./run-allsky-example.py - python3 ./compare-to-reference.py -- cd ../.. +- make clean +# +# Build and regression tests +# +- cd ${RRTMGP_ROOT}/tests|| exit 1 +- make clean +- make || exit 1 +- cp ${RRTMGP_ROOT}/examples/rfmip-clear-sky/multiple_input4MIPs_radiation_RFMIP_UColorado-RFMIP-1-2_none.nc test_atmospheres.nc +- ./clear_sky_regression test_atmospheres.nc ${RRTMGP_ROOT}/rrtmgp/data/rrtmgp-data-lw-g256-2018-12-04.nc +- ./clear_sky_regression test_atmospheres.nc ${RRTMGP_ROOT}/rrtmgp/data/rrtmgp-data-sw-g224-2018-12-04.nc +- python3 verification.py notifications: email: on_success: never diff --git a/Contributing.md b/Contributing.md new file mode 100644 index 000000000..e5a81aaa0 --- /dev/null +++ b/Contributing.md @@ -0,0 +1,19 @@ +### Contributing to RTE+RRTMGP + +Thanks for considering making a contribution to RTE+RRTMGP. + +The code in this repository is intended to work with compilers supporting the Fortran 2008 standard. It is also expected to run end-to-end on GPUs when compiled with OpenACC. Commits are tested with [Travis](https://travis-ci.com) against gfortran versions 8 and 9 and against various versions > 19.9 of the PGI compiler using resources provided by the [Swiss Supercomputing Center](https://cscs.ch). The testing uses two general codes in `examples/`for which results are compared against existing implemetations, and custom codes in tests/ intended to excercise all code options. + +##### Did you find a bug? + +Please file an issue on the [Github page](https://github.com/RobertPincus/rte-rrtmgp/issues). + +##### Did you write a patch that fixes a bug? + +Please fork this repository, branch from `develop`, make your changes, and open a Github [pull request](https://github.com/RobertPincus/rte-rrtmgp/pulls) against branch `develop`. + +##### Did you add functionality? + +Please fork this repository, branch from `develop`, make your changes, and open a Github [pull request](https://github.com/RobertPincus/rte-rrtmgp/pulls) against branch `develop`, adding a new regression test or comparison against the reference in `tests/verification.py` or `tests/validation-plots.py` as appropriate. + +RTE+RRTMGP is intended to be a core that users can extend with custom code to suit their own needs. \ No newline at end of file diff --git a/README.md b/README.md index 15128baa1..486c9a182 100644 --- a/README.md +++ b/README.md @@ -6,9 +6,16 @@ RRTMGP uses a k-distribution to provide an optical description (absorption and p RTE computes fluxes given spectrally-resolved optical descriptions and source functions. The fluxes are normally summarized or reduced via a user extensible class. -Example programs and documenation are evolving - please see examples/ in the repo and Wiki on the project's Github page. Suggestions are welcome. Meanwhile for questions please contact Robert Pincus and Eli Mlawer at rrtmgp@aer.com. +Example programs and documentation are evolving - please see examples/ in the repo and Wiki on the project's Github page. Suggestions are welcome. Meanwhile for questions please contact Robert Pincus and Eli Mlawer at rrtmgp@aer.com. + +## Recent changes + +1. The default method for solution for longwave problems that include scattering has been changed from 2-stream methods to a re-scaled and refined no-scattering calculation following [Tang et al. 2018](https://doi.org/10.1175/JAS-D-18-0014.1). +2. In RRTMGP gas optics, the spectrally-resolved solar source function in can be adjusted by specifying the total solar irradiance (`gas_optics%set_tsi(tsi)`) and/or the facular and sunspot indicies (`gas_optics%set_solar_variability(mg_index, sb_index, tsi)`)from the [NRLSSI2 model of solar variability](http://doi.org/10.1175/BAMS-D-14-00265.1). +3. `rte_lw()` now includes optional arguments for computing the Jacobian (derivative) of broadband flux with respect to changes in surface temperature. In calculations neglecting scattering only the Jacobian of upwelling flux is computed. When using re-scaling to account for scattering the Jacobians of both up- and downwelling flux are computed. + +Relative to commit `69d36c9` to `master` on Apr 20, 2020, the required arguments to both the longwave and shortwave versions of `ty_gas_optics_rrtmgp%load()`have changed. -In the most recent revision, the default method for solution for longwave problems that include scattering has been changed from 2-stream methods to a re-scaled and refined no-scattering calculation following [Tang et al. 2018](https://doi.org/10.1175/JAS-D-18-0014.1). ## Building the libraries. diff --git a/azure-pipelines.yml b/azure-pipelines.yml index c775826e0..48bb0b79b 100644 --- a/azure-pipelines.yml +++ b/azure-pipelines.yml @@ -4,23 +4,35 @@ jobs: strategy: matrix: pgi_19_9_gpu: - compiler_module: PGI/19.9.0 + compiler_module: PGI/19.9 + PrgEnv: PrgEnv-pgi accel_module: cudatoolkit - FCFLAGS: "-O3 -ta=tesla:cc60,cuda10.1 -Minfo -Mallocatable=03 -gopt -Minline,reshape,maxsize:40" + FCFLAGS: "-O3 -ta=tesla:cc60,cuda10.1 -Mallocatable=03 -gopt -Minline,reshape,maxsize:40" RTE_KERNELS: openacc pgi_default_gpu: compiler_module: pgi + PrgEnv: PrgEnv-pgi accel_module: craype-accel-nvidia60 # Generic accelerator flag - FCFLAGS: "-O3 -acc -Minfo -Mallocatable=03 -gopt" + FCFLAGS: "-O3 -acc -Mallocatable=03 -gopt" RTE_KERNELS: openacc + pgi_19_10_cpu: + compiler_module: PGI/19.10 + PrgEnv: PrgEnv-pgi + accel_module: + # Error checking flags + FCFLAGS: "-Mallocatable=03 -Mstandard -Mbounds -Mchkptr -Kieee -Mchkstk" + RUN_CMD: pgi_19_9_cpu: - compiler_module: PGI/19.9.0 + compiler_module: PGI/19.9 + PrgEnv: PrgEnv-pgi accel_module: # Error checking flags FCFLAGS: "-Mallocatable=03 -Mstandard -Mbounds -Mchkptr -Kieee -Mchkstk" + RUN_CMD: pgi_default_cpu: compiler_module: pgi + PrgEnv: PrgEnv-pgi accel_module: # Error checking flags FCFLAGS: "-Mallocatable=03 -Mstandard -Mbounds -Mchkptr -Kieee -Mchkstk" @@ -32,21 +44,21 @@ jobs: steps: - script: | set -e - + echo " export PATH=$CRAY_BINUTILS_BIN:$PATH module load daint-gpu - module swap PrgEnv-cray PrgEnv-pgi + module swap PrgEnv-cray $(PrgEnv) module load cray-netcdf cray-hdf5 module swap pgi $(compiler_module) module load $(accel_module) module load cray-python/3.6.5.7 export LD_LIBRARY_PATH=$CRAY_LD_LIBRARY_PATH:$LD_LIBRARY_PATH - + echo Environment: module list " > modules - + displayName: 'Create module environment' - script: | set -e @@ -60,29 +72,41 @@ jobs: set -e source modules export RRTMGP_ROOT=$PWD - export RRTMGP_DIR=$PWD/build + export RRTMGP_BUILD=$PWD/build export FC=ftn + make -C build/ clean + make -C build/ -j 8 + make -C tests clean + make -C tests -j 1 + make -C examples/all-sky clean make -C examples/all-sky -j 8 + make -C examples/rfmip-clear-sky clean make -C examples/rfmip-clear-sky -j 8 displayName: 'Make' - script: | set -e source modules - cd examples/rfmip-clear-sky + export RRTMGP_ROOT=$PWD + cd ${RRTMGP_ROOT}/examples/rfmip-clear-sky srun -C gpu -A c15 -p cscsci python ./run-rfmip-examples.py --block_size 1800 - cd ../.. - cd examples/all-sky + cd ${RRTMGP_ROOT}/examples/all-sky srun -C gpu -A c15 -p cscsci python ./run-allsky-example.py + cd ${RRTMGP_ROOT}/tests + cp ${RRTMGP_ROOT}/examples/rfmip-clear-sky/multiple_input4MIPs_radiation_RFMIP_UColorado-RFMIP-1-2_none.nc test_atmospheres.nc + srun -C gpu -A c15 -p cscsci ./clear_sky_regression test_atmospheres.nc ${RRTMGP_ROOT}/rrtmgp/data/rrtmgp-data-lw-g256-2018-12-04.nc + srun -C gpu -A c15 -p cscsci ./clear_sky_regression test_atmospheres.nc ${RRTMGP_ROOT}/rrtmgp/data/rrtmgp-data-sw-g224-2018-12-04.nc displayName: 'Run' - script: | set -e source modules + export RRTMGP_ROOT=$PWD # This module will unload some of the build modules, so do the checks separately module load netcdf-python - cd examples/rfmip-clear-sky - python ./compare-to-reference.py --fail=1.e-4 - cd ../.. - cd examples/all-sky + cd ${RRTMGP_ROOT}/examples/rfmip-clear-sky + python ./compare-to-reference.py --fail=7.e-4 + cd ${RRTMGP_ROOT}/examples/all-sky python ./compare-to-reference.py + cd ${RRTMGP_ROOT}/tests + python verification.py displayName: 'Check results' diff --git a/build/Makefile.conf.ifort b/build/Makefile.conf.ifort index b2b47bfb3..325892b63 100644 --- a/build/Makefile.conf.ifort +++ b/build/Makefile.conf.ifort @@ -8,11 +8,11 @@ export FC = ifort # # Optimized # -export FCFLAGS += -m64 -O3 -g -traceback -assume realloc_lhs -extend-source 132 +export FCFLAGS += -m64 -O3 -g -traceback -heap-arrays -assume realloc_lhs -extend-source 132 export F77FLAGS += -m64 -O3 -g -traceback # can add -qopt-report-phase=vec # # Debugging # -# export FCFLAGS = -m64 -g -traceback -assume realloc_lhs -extend-source 132 -check bounds,uninit,pointers,stack -stand f03 +# export FCFLAGS = -m64 -g -traceback -heap-arrays -assume realloc_lhs -extend-source 132 -check bounds,uninit,pointers,stack -stand f08 # export F77FLAGS = -m64 -g -traceback -check bounds,uninit,pointers,stack diff --git a/examples/all-sky/compare-to-reference.py b/examples/all-sky/compare-to-reference.py index fb0813b6d..554dddb12 100644 --- a/examples/all-sky/compare-to-reference.py +++ b/examples/all-sky/compare-to-reference.py @@ -31,7 +31,7 @@ # If a version of the file exists in the reference directory, no need to download (can be over-ridden) if (args.download_reference or not os.path.exists(ref_file)): os.makedirs(args.ref_dir, exist_ok=True) - urllib.request.urlretrieve("https://owncloud.gwdg.de/index.php/s/wgCL7imbA0QRCEf/download", ref_file) + urllib.request.urlretrieve("https://owncloud.gwdg.de/index.php/s/OjbNzRTlXUk0G5w/download", ref_file) tst = xr.open_dataset(tst_file) ref = xr.open_dataset(ref_file) @@ -42,8 +42,10 @@ if np.any(np.isnan(tst.variables[v].values)): raise Exception("Some test values are missing. Now that is strange.") - diff = abs((tst-ref).variables[v].values) - avg = 0.5*(tst+ref).variables[v].values + # express as (tst-ref).variables[v].values when replacing reference file + # to have same number of columns + diff = abs((tst.variables[v]-ref.variables[v]).values) + avg = 0.5*(tst.variables[v]+ref.variables[v]).values # Division raises a runtime warning when we divide by zero even if the # values in those locations will be ignored. with np.errstate(divide='ignore', invalid='ignore'): diff --git a/examples/all-sky/garand-atmos-1.nc b/examples/all-sky/garand-atmos-1.nc new file mode 100644 index 000000000..fa92c2ea0 Binary files /dev/null and b/examples/all-sky/garand-atmos-1.nc differ diff --git a/examples/all-sky/rrtmgp-allsky.nc b/examples/all-sky/rrtmgp-allsky.nc deleted file mode 100644 index a9bf2b0ae..000000000 Binary files a/examples/all-sky/rrtmgp-allsky.nc and /dev/null differ diff --git a/examples/all-sky/run-allsky-example.py b/examples/all-sky/run-allsky-example.py index 3f2c08de0..8160a3e50 100644 --- a/examples/all-sky/run-allsky-example.py +++ b/examples/all-sky/run-allsky-example.py @@ -2,11 +2,14 @@ # # This script runs runs the RTE+RRTMGP all-sky examples # -import os, subprocess, glob +import os +import shutil +import subprocess +import glob import argparse rte_rrtmgp_dir = os.path.join("..", "..") -# This files lives in $RRTMGP_DIR/examples/all-sky/ +# This files lives in $RRTMGP_ROOT/examples/all-sky/ all_sky_dir = "." # Code should be run in the all_sky_dir directory @@ -18,8 +21,10 @@ # In the local directory all_sky_exe_name = os.path.join(all_sky_dir, "rrtmgp_allsky") +input_file = os.path.join(all_sky_dir, "garand-atmos-1.nc") atmos_file = os.path.join(all_sky_dir, "rrtmgp-allsky.nc") + if __name__ == '__main__': parser = argparse.ArgumentParser(description="Runs all-sky examples, resetting output.") parser.add_argument("--run_command", type=str, default="", @@ -38,7 +43,7 @@ os.chdir(all_sky_dir) # Remove cloudy-sky fluxes from the file containing the atmospheric profiles - subprocess.run(["git", "checkout", atmos_file]) + shutil.copyfile(input_file, atmos_file) subprocess.run([all_sky_exe_name, atmos_file, lw_gas_coeffs_file, lw_clouds_coeff_file, ncol_str, nloops_str]) subprocess.run([all_sky_exe_name, atmos_file, sw_gas_coeffs_file, sw_clouds_coeff_file, ncol_str, nloops_str]) diff --git a/examples/mo_load_coefficients.F90 b/examples/mo_load_coefficients.F90 index edf2ea4a5..6bc0c4386 100644 --- a/examples/mo_load_coefficients.F90 +++ b/examples/mo_load_coefficients.F90 @@ -70,9 +70,12 @@ subroutine load_and_init(kdist, filename, available_gases) real(wp), dimension(:,:,:), allocatable :: kminor_lower, kminor_upper real(wp), dimension(:,:,: ), allocatable :: rayl_lower, rayl_upper - real(wp), dimension(: ), allocatable :: solar_src + real(wp), dimension(: ), allocatable :: solar_quiet, solar_facular, solar_sunspot + real(wp) :: tsi_default, mg_default, sb_default real(wp), dimension(:,: ), allocatable :: totplnk real(wp), dimension(:,:,:,:), allocatable :: planck_frac + real(wp), dimension(:,:) , allocatable :: optimal_angle_fit + ! ----------------- ! ! Book-keeping variables @@ -92,7 +95,8 @@ subroutine load_and_init(kdist, filename, available_gases) nminor_absorber_intervals_upper, & ncontributors_lower, & ncontributors_upper, & - ninternalSourcetemps + ninternalSourcetemps, & + nfit_coeffs ! -------------------------------------------------- ! ! How big are the various arrays? @@ -117,6 +121,8 @@ subroutine load_and_init(kdist, filename, available_gases) = get_dim_size(ncid,'temperature_Planck') ncontributors_lower = get_dim_size(ncid,'contributors_lower') ncontributors_upper = get_dim_size(ncid,'contributors_upper') + nfit_coeffs = get_dim_size(ncid,'fit_coeffs') ! Will be 0 for SW + ! ----------------- ! ! Read the many arrays @@ -177,6 +183,7 @@ subroutine load_and_init(kdist, filename, available_gases) ! totplnk = read_field(ncid, 'totplnk', ninternalSourcetemps, nbnds) planck_frac = read_field(ncid, 'plank_fraction', ngpts, nmixingfracs, npress+1, ntemps) + optimal_angle_fit = read_field(ncid, 'optimal_angle_fit', nfit_coeffs, nbnds) call stop_on_err(kdist%load(available_gases, & gas_names, & key_species, & @@ -200,12 +207,18 @@ subroutine load_and_init(kdist, filename, available_gases) kminor_start_lower, & kminor_start_upper, & totplnk, planck_frac, & - rayl_lower, rayl_upper)) + rayl_lower, rayl_upper, & + optimal_angle_fit)) else ! ! Solar source doesn't have an dependencies yet ! - solar_src = read_field(ncid, 'solar_source', ngpts) + solar_quiet = read_field(ncid, 'solar_source_quiet', ngpts) + solar_facular = read_field(ncid, 'solar_source_facular', ngpts) + solar_sunspot = read_field(ncid, 'solar_source_sunspot', ngpts) + tsi_default = read_field(ncid, 'tsi_default') + mg_default = read_field(ncid, 'mg_default') + sb_default = read_field(ncid, 'sb_default') call stop_on_err(kdist%load(available_gases, & gas_names, & key_species, & @@ -228,7 +241,8 @@ subroutine load_and_init(kdist, filename, available_gases) scale_by_complement_upper, & kminor_start_lower, & kminor_start_upper, & - solar_src, & + solar_quiet, solar_facular, solar_sunspot, & + tsi_default, mg_default, sb_default, & rayl_lower, rayl_upper)) end if ! -------------------------------------------------- diff --git a/examples/rfmip-clear-sky/Makefile b/examples/rfmip-clear-sky/Makefile index 47f106971..1a47d561d 100644 --- a/examples/rfmip-clear-sky/Makefile +++ b/examples/rfmip-clear-sky/Makefile @@ -1,15 +1,15 @@ # -# Here set variables RRTMGP_DIR, NCHOME, NFHOME, TIME_DIR (for GPTL) +# Here set variables RRTMGP_BUILD, NCHOME, NFHOME, TIME_DIR (for GPTL) # or have those variables set in the environment # -include Makefile.libs --include $(RRTMGP_DIR)/Makefile.conf +-include $(RRTMGP_BUILD)/Makefile.conf # # RRTMGP library, module files # -LDFLAGS += -L$(RRTMGP_DIR) +LDFLAGS += -L$(RRTMGP_BUILD) LIBS += -lrrtmgp -lrte -FCINCLUDE += -I$(RRTMGP_DIR) +FCINCLUDE += -I$(RRTMGP_BUILD) # # netcdf library, module files @@ -53,11 +53,11 @@ 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_DIR)/librte.a $(RRTMGP_DIR)/librrtmgp.a +rrtmgp_rfmip_lw: rrtmgp_rfmip_lw.o $(ADDITIONS) $(RRTMGP_BUILD)/librte.a $(RRTMGP_BUILD)/librrtmgp.a rrtmgp_rfmip_lw.o: rrtmgp_rfmip_lw.F90 $(ADDITIONS) -rrtmgp_rfmip_sw: rrtmgp_rfmip_sw.o $(ADDITIONS) $(RRTMGP_DIR)/librte.a $(RRTMGP_DIR)/librrtmgp.a +rrtmgp_rfmip_sw: rrtmgp_rfmip_sw.o $(ADDITIONS) $(RRTMGP_BUILD)/librte.a $(RRTMGP_BUILD)/librrtmgp.a rrtmgp_rfmip_sw.o: rrtmgp_rfmip_sw.F90 $(ADDITIONS) diff --git a/examples/rfmip-clear-sky/Makefile.libs.macos b/examples/rfmip-clear-sky/Makefile.libs.macos index 135fea006..465cdd3aa 100644 --- a/examples/rfmip-clear-sky/Makefile.libs.macos +++ b/examples/rfmip-clear-sky/Makefile.libs.macos @@ -1,5 +1,5 @@ # Location of RTE+RRTMGP libraries, module files. -export RRTMGP_DIR = ../../build +export RRTMGP_BUILD = $(RRTMGP_ROOT)/build # Sets macros FC, FCFLAGS consistent with RTE+RRTMGP # NetCDF C and Fortran libraries, module files diff --git a/examples/rfmip-clear-sky/Makefile.libs.olcf b/examples/rfmip-clear-sky/Makefile.libs.olcf index 9754e1fbb..38ebdc75c 100644 --- a/examples/rfmip-clear-sky/Makefile.libs.olcf +++ b/examples/rfmip-clear-sky/Makefile.libs.olcf @@ -1,5 +1,5 @@ # Location of RTE+RRTMGP libraries, module files. -export RRTMGP_DIR = ../../build +export RRTMGP_BUILD = ../../build # Sets macros FC, FCFLAGS consistent with RTE+RRTMGP # NetCDF C and Fortran libraries, module files @@ -7,4 +7,4 @@ NCHOME = $(OLCF_NETCDF_ROOT) NFHOME = $(OLCF_NETCDF_FORTRAN_ROOT) # GPTL libraries and module files if desired (https://jmrosinski.github.io/GPTL/) -# export TIME_DIR +# export TIME_DIR diff --git a/examples/rfmip-clear-sky/stage_files.py b/examples/rfmip-clear-sky/stage_files.py index a1556e2ab..0009bc48e 100755 --- a/examples/rfmip-clear-sky/stage_files.py +++ b/examples/rfmip-clear-sky/stage_files.py @@ -2,13 +2,11 @@ # # This script downloads and creates files needed for the RFMIP off-line test cases # -import os, subprocess, glob -from shutil import copy2 +import os +import subprocess +import glob import urllib.request -# Will be needed by scripts to generate output file templates -from netCDF4 import Dataset -import time, uuid, argparse -import json +# This script invokes templ_scr_url, which uses Python modules time, uuid, argparse, netCDF4 # # Download and/or create input files and output template files diff --git a/extensions/solar_variability/mo_solar_variability.F90 b/extensions/solar_variability/mo_solar_variability.F90 new file mode 100644 index 000000000..93a126bcf --- /dev/null +++ b/extensions/solar_variability/mo_solar_variability.F90 @@ -0,0 +1,185 @@ +! This code is part of RRTM for GCM Applications - Parallel (RRTMGP) +! +! Contacts: Robert Pincus and Eli Mlawer +! email: rrtmgp@aer.com +! +! Copyright 2015-2020, Atmospheric and Environmental Research and +! Regents of the University of Colorado. All right reserved. +! +! Use and duplication is permitted under the terms of the +! BSD 3-clause license, see http://opensource.org/licenses/BSD-3-Clause +! ------------------------------------------------------------------------------------------------- +! Description: Optional calculation of solar variability facular and +! sunspot indices + +module mo_solar_variability + use mo_rte_kind, only: wp + + implicit none + private + type, public :: ty_solar_var + ! + ! Data + ! + real(wp), dimension(:,:), allocatable :: avgcyc_ind ! solar variabilty index lookup table + ! time-averaged over solar cycles 13-24. + ! (NRLSSI2 facular "Bremen" index and + ! sunspot "SPOT67" index) + ! (nsolarterms, nsolarfrac) -> (2,134) + contains + ! + ! Public procedures + ! + procedure, public :: solar_var_ind_interp + procedure, public :: load + procedure, public :: finalize + ! + end type ty_solar_var + +contains + ! ------------------------------------------------------------------------------ + ! + ! Routine to load mean facular and sunspot index tables + ! + ! ------------------------------------------------------------------------------ + function load(this, avgcyc_ind) result(error_msg) + class(ty_solar_var), intent(inout) :: this + ! Lookup table of mean solar cycle facular brightening and sunspot dimming indices + real(wp), dimension(:,:), intent(in ) :: avgcyc_ind + character(len=128) :: error_msg + ! ------- + ! + ! Local variables + ! + integer :: nsolarterms, nsolarfrac + + error_msg = "" + ! + ! LUT index dimensions + ! + nsolarterms= size(avgcyc_ind,dim=1) + nsolarfrac = size(avgcyc_ind,dim=2) + ! + ! Allocate LUT index array + allocate(this%avgcyc_ind(nsolarterms, nsolarfrac)) + + ! Load LUT index array + this%avgcyc_ind = avgcyc_ind + + end function load + !-------------------------------------------------------------------------------------------------------------------- + ! + ! Finalize + ! + !-------------------------------------------------------------------------------------------------------------------- + subroutine finalize(this) + class(ty_solar_var), intent(inout) :: this + + ! Lookup table solar variability indices + if(allocated(this%avgcyc_ind)) then + deallocate(this%avgcyc_ind) + end if + + end subroutine finalize + ! ------------------------------------------------------------------------------ + ! + ! Facular brightening and sunspot dimming indices are derived from the + ! averaged solar cycle, which is the mean of Solar Cycles 13-24. The user specifices + ! the solar cycle fraction (0 to 1) and the indices are interpolated to the + ! requested fractional position within the cycle, where 0 is close to solar minimum. + ! + function solar_var_ind_interp(this, & + solcycfrac, & + mg_index, sb_index) & + result(error_msg) + + class(ty_solar_var), intent(in ) :: this ! Solar variability + real(wp), intent(in ) :: solcycfrac ! solar cycle fraction + real(wp), intent(out ) :: mg_index ! Facular brightening NRLSSI2 index + ! interpolated from the mean solar cycle + ! to the provided solar cycle fraction + real(wp), intent(out ) :: sb_index ! Sunspot dimmng NRLSSI2 index + ! interpolated from the mean solar cycle + ! to the provided solar cycle fraction + character(len=128) :: error_msg + + ! ---------------------------------------------------------- + ! Local variables + ! + integer :: nsolfrac ! Number of solar fraction points in facular + ! and sunspot tables + integer :: sfid ! Solar variability solar cycle fraction index + + real(wp) :: intrvl_len ! Fractional interval length of mgavgcyc and sbavgcyc + real(wp) :: intrvl_len_hf ! Fractional half interval length of mgavgcyc and sbavgcyc + real(wp) :: fraclo, frachi, intfrac ! Solar variability interpolation factors + + ! ---------------------------------------------------------- + ! Error checking + error_msg = "" + ! + ! Check input data sizes and values + ! + if (solcycfrac .lt. 0._wp .or. solcycfrac .gt. 1._wp) & + error_msg = 'solar_var_ind_interp: solcycfrac out of range' + if(error_msg /= '') return + ! + ! Interpolate solar variability indices to requested solar cycle fraction, + ! and derive final facular and sunspot indices + ! + ! nsolfrac is the length of the time dimension of the interpolation tables + ! of facular and sunspot indices over the mean solar cycle (this%avgcyc_ind). + ! The end-points of avgcyc_ind represent the indices at solcycfrac values of + ! 0 (first day of the first year) and 1 (last day of the 11th year), while + ! the intervening values of avgcyc_ind represent the indices at the center + ! of each month over the mean 11-year solar cycle. + if (allocated (this%avgcyc_ind)) then + nsolfrac = size(this%avgcyc_ind,2) + ! Define indices for the lowest allowable value of solcycfrac + if (solcycfrac .eq. 0._wp) then + mg_index = this%avgcyc_ind(1,1) + sb_index = this%avgcyc_ind(2,1) + ! Define indices for the highest allowable value of solcycfrac + elseif (solcycfrac .eq. 1._wp) then + mg_index = this%avgcyc_ind(1,nsolfrac) + sb_index = this%avgcyc_ind(2,nsolfrac) + ! Define indices for intervening values of solcycfrac + else + intrvl_len = 1._wp / (nsolfrac-2) + intrvl_len_hf = 0.5_wp * intrvl_len + ! Define interpolation fractions for the first interval, which represents + ! the first half of the first month of the first year of the mean 11-year + ! solar cycle + if (solcycfrac .le. intrvl_len_hf) then + sfid = 1 + fraclo = 0._wp + frachi = intrvl_len_hf + endif + ! Define interpolation fractions for the intervening intervals, which represent + ! the center point of each month in each year of the mean 11-year solar cycle + if (solcycfrac .gt. intrvl_len_hf .and. solcycfrac .lt. 1._wp-intrvl_len_hf) then + sfid = floor((solcycfrac-intrvl_len_hf) * (nsolfrac-2)) + 2 + fraclo = (sfid-2) * intrvl_len + intrvl_len_hf + frachi = fraclo + intrvl_len + endif + ! Define interpolation fractions for the last interval, which represents + ! the last half of the last month of the last year of the mean 11-year + ! solar cycle + if (solcycfrac .ge. 1._wp-intrvl_len_hf) then + sfid = (nsolfrac-2) + 1 + fraclo = 1._wp - intrvl_len_hf + frachi = 1._wp + endif + ! Interpolate the facular (mg_index) and sunspot (sb_index) indices for the + ! requested value of solcycfrac + intfrac = (solcycfrac - fraclo) / (frachi - fraclo) + mg_index = this%avgcyc_ind(1,sfid) + & + intfrac * (this%avgcyc_ind(1,sfid+1) - this%avgcyc_ind(1,sfid)) + sb_index = this%avgcyc_ind(2,sfid) + & + intfrac * (this%avgcyc_ind(2,sfid+1) - this%avgcyc_ind(2,sfid)) + endif + endif + + end function solar_var_ind_interp + ! -------------------------------------------------------------------------------------- +end module mo_solar_variability diff --git a/extensions/solar_variability/rrtmgp-solar-var-tables.nc b/extensions/solar_variability/rrtmgp-solar-var-tables.nc new file mode 100644 index 000000000..3545eb2d0 Binary files /dev/null and b/extensions/solar_variability/rrtmgp-solar-var-tables.nc differ diff --git a/rrtmgp/data/rrtmgp-data-lw-g256-2018-12-04.nc b/rrtmgp/data/rrtmgp-data-lw-g256-2018-12-04.nc index 86f7a258d..991ce908a 100644 Binary files a/rrtmgp/data/rrtmgp-data-lw-g256-2018-12-04.nc and b/rrtmgp/data/rrtmgp-data-lw-g256-2018-12-04.nc differ diff --git a/rrtmgp/data/rrtmgp-data-sw-g224-2018-12-04.nc b/rrtmgp/data/rrtmgp-data-sw-g224-2018-12-04.nc index 63bc19ec5..03cb071b6 100644 Binary files a/rrtmgp/data/rrtmgp-data-sw-g224-2018-12-04.nc and b/rrtmgp/data/rrtmgp-data-sw-g224-2018-12-04.nc differ diff --git a/rrtmgp/kernels-openacc/mo_gas_optics_kernels.F90 b/rrtmgp/kernels-openacc/mo_gas_optics_kernels.F90 index 3e6f32a9d..e6e7996ed 100644 --- a/rrtmgp/kernels-openacc/mo_gas_optics_kernels.F90 +++ b/rrtmgp/kernels-openacc/mo_gas_optics_kernels.F90 @@ -538,7 +538,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) bind(C, name="compute_Planck_source") + sfc_src, lay_src, lev_src_inc, lev_src_dec, sfc_source_Jac) bind(C, name="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 @@ -561,8 +561,12 @@ subroutine compute_Planck_source( & real(wp), dimension(ngpt, ncol), intent(out) :: sfc_src real(wp), dimension(ngpt,nlay,ncol), intent(out) :: lay_src real(wp), dimension(ngpt,nlay,ncol), intent(out) :: lev_src_inc, lev_src_dec + + real(wp), dimension(ngpt, ncol), 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] @@ -573,6 +577,7 @@ subroutine compute_Planck_source( & !$acc enter data copyin(tlay,tlev,tsfc,fmajor,jeta,tropo,jtemp,jpress,gpoint_bands,pfracin,totplnk,gpoint_flavor) !$acc enter data create(sfc_src,lay_src,lev_src_inc,lev_src_dec) !$acc enter data create(pfrac,planck_function) + !$acc enter data create(sfc_source_Jac) ! Calculation of fraction of band's Planck irradiance associated with each g-point !$acc parallel loop collapse(3) @@ -596,7 +601,8 @@ subroutine compute_Planck_source( & ! !$acc parallel loop do icol = 1, ncol - call interpolate1D(tsfc(icol), temp_ref_min, totplnk_delta, totplnk, planck_function(1:nbnd,1,icol)) + call interpolate1D(tsfc(icol) , temp_ref_min, totplnk_delta, totplnk, planck_function(1:nbnd,1,icol)) + call interpolate1D(tsfc(icol) + delta_Tsurf, temp_ref_min, totplnk_delta, totplnk, planck_function(1:nbnd,2,icol)) end do ! ! Map to g-points @@ -604,7 +610,9 @@ subroutine compute_Planck_source( & !$acc parallel loop collapse(2) do igpt = 1, ngpt do icol = 1, ncol - sfc_src(igpt,icol) = pfrac(igpt,sfc_lay,icol) * planck_function(gpoint_bands(igpt), 1, icol) + sfc_src (igpt,icol) = pfrac(igpt,sfc_lay,icol) * planck_function(gpoint_bands(igpt),1,icol) + sfc_source_Jac(igpt,icol) = pfrac(igpt,sfc_lay,icol) * & + (planck_function(gpoint_bands(igpt),2,icol) - planck_function(gpoint_bands(igpt),1,icol)) end do end do ! icol @@ -667,6 +675,7 @@ subroutine compute_Planck_source( & !$acc exit data delete(tlay,tlev,tsfc,fmajor,jeta,tropo,jtemp,jpress,gpoint_bands,pfracin,totplnk,gpoint_flavor) !$acc exit data delete(pfrac,planck_function) !$acc exit data copyout(sfc_src,lay_src,lev_src_inc,lev_src_dec) + !$acc exit data copyout(sfc_source_Jac) end subroutine compute_Planck_source ! ---------------------------------------------------------- diff --git a/rrtmgp/kernels/mo_gas_optics_kernels.F90 b/rrtmgp/kernels/mo_gas_optics_kernels.F90 index 9b9f99e7e..1c11ab36c 100644 --- a/rrtmgp/kernels/mo_gas_optics_kernels.F90 +++ b/rrtmgp/kernels/mo_gas_optics_kernels.F90 @@ -481,7 +481,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) bind(C, name="compute_Planck_source") + sfc_src, lay_src, lev_src_inc, lev_src_dec, sfc_source_Jac) bind(C, name="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 @@ -504,8 +504,12 @@ subroutine compute_Planck_source( & real(wp), dimension(ngpt, ncol), intent(out) :: sfc_src real(wp), dimension(ngpt,nlay,ncol), intent(out) :: lay_src real(wp), dimension(ngpt,nlay,ncol), intent(out) :: lev_src_inc, lev_src_dec + + real(wp), dimension(ngpt, ncol), 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] @@ -536,7 +540,8 @@ subroutine compute_Planck_source( & ! Compute surface source irradiance for g-point, equals band irradiance x fraction for g-point ! do icol = 1, ncol - planck_function(1:nbnd,1,icol) = interpolate1D(tsfc(icol), temp_ref_min, totplnk_delta, totplnk) + planck_function(1:nbnd,1,icol) = interpolate1D(tsfc(icol) , temp_ref_min, totplnk_delta, totplnk) + planck_function(1:nbnd,2,icol) = interpolate1D(tsfc(icol) + delta_Tsurf, temp_ref_min, totplnk_delta, totplnk) ! ! Map to g-points ! @@ -544,7 +549,9 @@ subroutine compute_Planck_source( & gptS = band_lims_gpt(1, ibnd) gptE = band_lims_gpt(2, ibnd) do igpt = gptS, gptE - sfc_src(igpt, icol) = pfrac(igpt,sfc_lay,icol) * planck_function(ibnd, 1, icol) + sfc_src (igpt, icol) = pfrac(igpt,sfc_lay,icol) * planck_function(ibnd,1,icol) + sfc_source_Jac(igpt, icol) = pfrac(igpt,sfc_lay,icol) * & + (planck_function(ibnd, 2, icol) - planck_function(ibnd,1,icol)) end do end do end do ! icol diff --git a/rrtmgp/mo_gas_optics_rrtmgp.F90 b/rrtmgp/mo_gas_optics_rrtmgp.F90 index e89183579..f2e8f2d1c 100644 --- a/rrtmgp/mo_gas_optics_rrtmgp.F90 +++ b/rrtmgp/mo_gas_optics_rrtmgp.F90 @@ -127,11 +127,22 @@ module mo_gas_optics_rrtmgp ! planck_frac(g-point, eta, pressure, temperature) real(wp), dimension(:,:), allocatable :: totplnk ! integrated Planck irradiance by band; (Planck temperatures,band) real(wp) :: totplnk_delta ! temperature steps in totplnk + real(wp), dimension(:,:), allocatable :: optimal_angle_fit ! coefficients of linear function + ! of vertical path clear-sky transmittance that is used to + ! determine the secant of single angle used for the + ! no-scattering calculation, + ! optimal_angle_fit(coefficient, band) ! ----------------------------------------------------------------------------------- - ! Solar source function spectral mapping - ! Allocated only when gas optics object is external-source + ! Solar source function spectral mapping with solar variability capability + ! Allocated when gas optics object is external-source + ! n-solar-terms: quiet sun, facular brightening and sunspot dimming components + ! following the NRLSSI2 model of Coddington et al. 2016, doi:10.1175/BAMS-D-14-00265.1. ! - real(wp), dimension(:), allocatable :: solar_src ! incoming solar irradiance(g-point) + real(wp), dimension(:), allocatable :: solar_source ! incoming solar irradiance, computed from other three terms (g-point) + real(wp), dimension(:), allocatable :: solar_source_quiet ! incoming solar irradiance, quiet sun term (g-point) + real(wp), dimension(:), allocatable :: solar_source_facular ! incoming solar irradiance, facular term (g-point) + real(wp), dimension(:), allocatable :: solar_source_sunspot ! incoming solar irradiance, sunspot term (g-point) + ! ! ----------------------------------------------------------------------------------- ! Ancillary @@ -153,13 +164,15 @@ module mo_gas_optics_rrtmgp procedure, public :: get_press_max procedure, public :: get_temp_min procedure, public :: get_temp_max + procedure, public :: compute_optimal_angles + procedure, public :: set_solar_variability + procedure, public :: set_tsi ! Internal procedures procedure, private :: load_int procedure, private :: load_ext procedure, public :: gas_optics_int procedure, public :: gas_optics_ext procedure, private :: check_key_species_present - procedure, private :: get_minor_list ! Interpolation table dimensions procedure, private :: get_nflav procedure, private :: get_neta @@ -185,7 +198,7 @@ module mo_gas_optics_rrtmgp pure function get_ngas(this) ! return the number of gases registered in the spectral configuration class(ty_gas_optics_rrtmgp), intent(in) :: this - integer :: get_ngas + integer :: get_ngas get_ngas = size(this%gas_names) end function get_ngas @@ -257,7 +270,7 @@ function gas_optics_int(this, & ! External source -- check arrays sizes and values ! input data sizes and values ! - !$acc enter data copyin(tsfc,tlev) + !$acc enter data copyin(tsfc) if(.not. extents_are(tsfc, ncol)) & error_msg = "gas_optics(): array tsfc has wrong size" if(any_vals_outside(tsfc, this%temp_ref_min, this%temp_ref_max)) & @@ -265,6 +278,7 @@ function gas_optics_int(this, & if(error_msg /= '') return if(present(tlev)) then + !$acc enter data copyin(tlev) if(.not. extents_are(tlev, ncol, nlay+1)) & error_msg = "gas_optics(): array tlev has wrong size" if(any_vals_outside(tlev, this%temp_ref_min, this%temp_ref_max)) & @@ -282,13 +296,26 @@ function gas_optics_int(this, & ! ! Interpolate source function ! - error_msg = source(this, & - ncol, nlay, nband, ngpt, & - play, plev, tlay, tsfc, & - jtemp, jpress, jeta, tropo, fmajor, & - sources, & - tlev) - !$acc exit data delete(tsfc,tlev) + if(present(tlev)) then + ! + ! present status of optional argument should be passed to source() + ! but isn't with PGI 19.10 + ! + error_msg = source(this, & + ncol, nlay, nband, ngpt, & + play, plev, tlay, tsfc, & + jtemp, jpress, jeta, tropo, fmajor, & + sources, & + tlev) + !$acc exit data delete(tlev) + else + error_msg = source(this, & + ncol, nlay, nband, ngpt, & + play, plev, tlay, tsfc, & + jtemp, jpress, jeta, tropo, fmajor, & + sources) + end if + !$acc exit data delete(tsfc) !$acc exit data delete(jtemp, jpress, tropo, fmajor, jeta) end function gas_optics_int !------------------------------------------------------------------------------------------ @@ -356,7 +383,7 @@ function gas_optics_ext(this, & !$acc parallel loop collapse(2) do igpt = 1,ngpt do icol = 1,ncol - toa_src(icol,igpt) = this%solar_src(igpt) + toa_src(icol,igpt) = this%solar_source(igpt) end do end do !$acc exit data copyout(toa_src) @@ -442,6 +469,10 @@ function compute_gas_taus(this, & error_msg = "gas_optics(): array tlay has wrong size" if(.not. extents_are(plev, ncol, nlay+1)) & error_msg = "gas_optics(): array plev has wrong size" + if(optical_props%get_ncol() /= ncol .or. & + optical_props%get_nlay() /= nlay .or. & + optical_props%get_ngpt() /= ngpt) & + error_msg = "gas_optics(): optical properties have the wrong extents" if(error_msg /= '') return if(any_vals_outside(play, this%press_ref_min,this%press_ref_max)) & @@ -489,6 +520,7 @@ function compute_gas_taus(this, & ! Compute dry air column amounts (number of molecule per cm^2) if user hasn't provided them ! idx_h2o = string_loc_in_array('h2o', this%gas_names) + col_dry_wk => col_dry_arr !$acc enter data create(col_dry_wk, col_dry_arr, col_gas) if (present(col_dry)) then col_dry_wk => col_dry @@ -596,6 +628,83 @@ function compute_gas_taus(this, & end function compute_gas_taus !------------------------------------------------------------------------------------------ ! + ! Compute the spectral solar source function adjusted to account for solar variability + ! following the NRLSSI2 model of Coddington et al. 2016, doi:10.1175/BAMS-D-14-00265.1. + ! as specified by the facular brightening (mg_index) and sunspot dimming (sb_index) + ! indices provided as input. + ! + ! Users provide the NRLSSI2 facular ("Bremen") index and sunspot ("SPOT67") index. + ! Changing either of these indicies will change the total solar irradiance (TSI) + ! Code in extensions/mo_solar_variability may be used to compute the value of these + ! indices through an average solar cycle + ! Users may also specify the TSI, either alone or in conjunction with the facular and sunspot indices + ! + !------------------------------------------------------------------------------------------ + function set_solar_variability(this, & + mg_index, sb_index, tsi) & + result(error_msg) + ! + ! Updates the spectral distribution and, optionally, + ! the integrated value of the solar source function + ! Modifying either index will change the total solar irradiance + ! + class(ty_gas_optics_rrtmgp), intent(inout) :: this + ! + real(wp), intent(in) :: mg_index, & ! facular brightening index (NRLSSI2 facular "Bremen" index) + sb_index ! sunspot dimming index (NRLSSI2 sunspot "SPOT67" index) + real(wp), optional, intent(in) :: tsi ! total solar irradiance + character(len=128) :: error_msg + ! ---------------------------------------------------------- + integer :: igpt + real(wp), parameter :: a_offset = 0.1495954_wp + real(wp), parameter :: b_offset = 0.00066696_wp + ! ---------------------------------------------------------- + error_msg = "" + if(mg_index < 0._wp) error_msg = 'mg_index out of range' + if(sb_index < 0._wp) error_msg = 'sb_index out of range' + if(error_msg /= "") return + ! + ! Calculate solar source function for provided facular and sunspot indices + ! + !$acc parallel loop + do igpt = 1, size(this%solar_source_quiet) + this%solar_source(igpt) = this%solar_source_quiet(igpt) + & + (mg_index - a_offset) * this%solar_source_facular(igpt) + & + (sb_index - b_offset) * this%solar_source_sunspot(igpt) + end do + ! + ! Scale solar source to input TSI value + ! + if (present(tsi)) error_msg = this%set_tsi(tsi) + + end function set_solar_variability + !------------------------------------------------------------------------------------------ + function set_tsi(this, tsi) result(error_msg) + ! + ! Scale the solar source function without changing the spectral distribution + ! + class(ty_gas_optics_rrtmgp), intent(inout) :: this + real(wp), intent(in ) :: tsi ! user-specified total solar irradiance; + character(len=128) :: error_msg + + real(wp) :: norm + ! ---------------------------------------------------------- + error_msg = "" + if(tsi < 0._wp) then + error_msg = 'tsi out of range' + else + ! + ! Scale the solar source function to the input tsi + ! + !$acc kernels + norm = 1._wp/sum(this%solar_source(:)) + this%solar_source(:) = this%solar_source(:) * tsi * norm + !$acc end kernels + end if + + end function set_tsi + !------------------------------------------------------------------------------------------ + ! ! Compute Planck source functions at layer centers and levels ! function source(this, & @@ -603,7 +712,7 @@ function source(this, & play, plev, tlay, tsfc, & jtemp, jpress, jeta, tropo, fmajor, & sources, & ! Planck sources - tlev) & ! optional input + tlev) & ! optional input result(error_msg) ! inputs class(ty_gas_optics_rrtmgp), intent(in ) :: this @@ -627,6 +736,8 @@ function source(this, & integer :: icol, ilay, igpt real(wp), dimension(ngpt,nlay,ncol) :: lay_source_t, lev_source_inc_t, lev_source_dec_t real(wp), dimension(ngpt, ncol) :: sfc_source_t + real(wp), dimension(ngpt, ncol) :: sfc_source_Jac + ! Variables for temperature at layer edges [K] (ncol, nlay+1) real(wp), dimension( ncol,nlay+1), target :: tlev_arr real(wp), dimension(:,:), pointer :: tlev_wk @@ -669,23 +780,32 @@ function source(this, & !$acc enter data copyin(sources) !$acc enter data create(sources%lay_source, sources%lev_source_inc, sources%lev_source_dec, sources%sfc_source) !$acc enter data create(sfc_source_t, lay_source_t, lev_source_inc_t, lev_source_dec_t) attach(tlev_wk) + !$acc enter data create(sfc_source_Jac) + !$acc enter data create(sources%sfc_source_Jac) call compute_Planck_source(ncol, nlay, nbnd, ngpt, & get_nflav(this), this%get_neta(), this%get_npres(), this%get_ntemp(), this%get_nPlanckTemp(), & tlay, tlev_wk, tsfc, merge(1,nlay,play(1,1) > play(1,nlay)), & 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, & - sfc_source_t, lay_source_t, lev_source_inc_t, lev_source_dec_t) + sfc_source_t, lay_source_t, lev_source_inc_t, lev_source_dec_t, & + sfc_source_Jac) !$acc parallel loop collapse(2) do igpt = 1, ngpt do icol = 1, ncol - sources%sfc_source(icol,igpt) = sfc_source_t(igpt,icol) + sources%sfc_source (icol,igpt) = sfc_source_t (igpt,icol) + sources%sfc_source_Jac(icol,igpt) = sfc_source_Jac(igpt,icol) end do end do call reorder123x321(lay_source_t, sources%lay_source) call reorder123x321(lev_source_inc_t, sources%lev_source_inc) call reorder123x321(lev_source_dec_t, sources%lev_source_dec) + ! + ! Transposition of a 2D array, for which we don't have a routine in mo_rrtmgp_util_reorder. + ! + !$acc exit data delete(sfc_source_Jac) !$acc exit data delete(sfc_source_t, lay_source_t, lev_source_inc_t, lev_source_dec_t) detach(tlev_wk) + !$acc exit data copyout(sources%sfc_source_Jac) !$acc exit data copyout(sources%lay_source, sources%lev_source_inc, sources%lev_source_dec, sources%sfc_source) !$acc exit data copyout(sources) end function source @@ -713,7 +833,9 @@ function load_int(this, available_gases, gas_names, key_species, & scale_by_complement_upper, & kminor_start_lower, & kminor_start_upper, & - totplnk, planck_frac, rayl_lower, rayl_upper) result(err_message) + totplnk, planck_frac, & + rayl_lower, rayl_upper, & + optimal_angle_fit) result(err_message) class(ty_gas_optics_rrtmgp), intent(inout) :: this class(ty_gas_concs), intent(in ) :: available_gases ! Which gases does the host model have available? character(len=*), dimension(:), intent(in ) :: gas_names @@ -729,6 +851,7 @@ function load_int(this, available_gases, gas_names, key_species, & real(wp), dimension(:,:,:,:), intent(in ) :: planck_frac real(wp), dimension(:,:,:), intent(in ), & allocatable :: rayl_lower, rayl_upper + real(wp), dimension(:,:), intent(in ) :: optimal_angle_fit character(len=*), dimension(:), intent(in ) :: gas_minor,identifier_minor character(len=*), dimension(:), intent(in ) :: minor_gases_lower, & minor_gases_upper @@ -767,13 +890,16 @@ function load_int(this, available_gases, gas_names, key_species, & rayl_lower, rayl_upper) ! Planck function tables ! - allocate(this%totplnk (size(totplnk, 1), size(totplnk, 2)), & - this%planck_frac(size(planck_frac,1), size(planck_frac,2), size(planck_frac,3), size(planck_frac,4)) ) - !$acc enter data create(this%totplnk, this%planck_frac) + allocate(this%totplnk (size(totplnk, 1), size(totplnk, 2)), & + this%planck_frac (size(planck_frac,1), size(planck_frac,2), size(planck_frac,3), size(planck_frac,4)), & + this%optimal_angle_fit(size(optimal_angle_fit, 1), size(optimal_angle_fit, 2))) + !$acc enter data create(this%totplnk, this%planck_frac, this%optimal_angle_fit) !$acc kernels this%totplnk = totplnk this%planck_frac = planck_frac + this%optimal_angle_fit = optimal_angle_fit !$acc end kernels + ! Temperature steps for Planck function interpolation ! Assumes that temperature minimum and max are the same for the absorption coefficient grid and the ! Planck grid and the Planck grid is equally spaced @@ -801,9 +927,11 @@ function load_ext(this, available_gases, gas_names, key_species, & scale_by_complement_upper, & kminor_start_lower, & kminor_start_upper, & - solar_src, rayl_lower, rayl_upper) result(err_message) + solar_quiet, solar_facular, solar_sunspot, & + tsi_default, mg_default, sb_default, & + rayl_lower, rayl_upper) result(err_message) class(ty_gas_optics_rrtmgp), intent(inout) :: this - class(ty_gas_concs), intent(in ) :: available_gases ! Which gases does the host model have available? + class(ty_gas_concs), intent(in ) :: available_gases ! Which gases does the host model have available? character(len=*), & dimension(:), intent(in) :: gas_names integer, dimension(:,:,:), intent(in) :: key_species @@ -823,22 +951,28 @@ function load_ext(this, available_gases, gas_names, key_species, & integer, dimension(:,:), intent(in) :: & minor_limits_gpt_lower, & minor_limits_gpt_upper - logical(wl), dimension(:), intent(in) :: & + logical(wl), dimension(:), intent(in) :: & minor_scales_with_density_lower, & minor_scales_with_density_upper - character(len=*), dimension(:),intent(in) :: & + character(len=*),dimension(:),intent(in) :: & scaling_gas_lower, & scaling_gas_upper - logical(wl), dimension(:), intent(in) :: & + logical(wl), dimension(:), intent(in) :: & scale_by_complement_lower, & scale_by_complement_upper - integer, dimension(:), intent(in) :: & + integer, dimension(:), intent(in) :: & kminor_start_lower, & kminor_start_upper - real(wp), dimension(:), intent(in), allocatable :: solar_src - ! allocatable status to change when solar source is present in file - real(wp), dimension(:,:,:), intent(in), allocatable :: rayl_lower, rayl_upper + real(wp), dimension(:), intent(in) :: solar_quiet, & + solar_facular, & + solar_sunspot + real(wp), intent(in) :: tsi_default, & + mg_default, sb_default + real(wp), dimension(:,:,:), intent(in), & + allocatable :: rayl_lower, rayl_upper character(len = 128) err_message + + integer :: ngpt ! ---- !$acc enter data create(this) err_message = init_abs_coeffs(this, & @@ -861,14 +995,20 @@ function load_ext(this, available_gases, gas_names, key_species, & kminor_start_lower, & kminor_start_upper, & rayl_lower, rayl_upper) + if(err_message /= "") return ! - ! Solar source table init + ! Spectral solar irradiance terms init ! - allocate(this%solar_src(size(solar_src))) - !$acc enter data create(this%solar_src) + ngpt = size(solar_quiet) + allocate(this%solar_source_quiet(ngpt), this%solar_source_facular(ngpt), & + this%solar_source_sunspot(ngpt), this%solar_source(ngpt)) + !$acc enter data create(this%solar_source_quiet, this%solar_source_facular, this%solar_source_sunspot, this%solar_source) !$acc kernels - this%solar_src = solar_src + this%solar_source_quiet = solar_quiet + this%solar_source_facular = solar_facular + this%solar_source_sunspot = solar_sunspot !$acc end kernels + err_message = this%set_solar_variability(mg_default, sb_default) end function load_ext !-------------------------------------------------------------------------------------------------------------------- ! @@ -1126,32 +1266,6 @@ function check_key_species_present(this, gas_desc) result(error_msg) end function check_key_species_present !-------------------------------------------------------------------------------------------------------------------- ! - ! Function to define names of key and minor gases to be used by gas_optics(). - ! The final list gases includes those that are defined in gas_optics_specification - ! and are provided in ty_gas_concs. - ! - function get_minor_list(this, gas_desc, ngas, names_spec) - class(ty_gas_optics_rrtmgp), intent(in) :: this - class(ty_gas_concs), intent(in) :: gas_desc - integer, intent(in) :: ngas - character(32), dimension(ngas), intent(in) :: names_spec - - ! List of minor gases to be used in gas_optics() - character(len=32), dimension(:), allocatable :: get_minor_list - ! Logical flag for minor species in specification (T = minor; F = not minor) - logical, dimension(size(names_spec)) :: gas_is_present - integer :: igas, icnt - - if (allocated(get_minor_list)) deallocate(get_minor_list) - do igas = 1, this%get_ngas() - gas_is_present(igas) = string_in_array(names_spec(igas), gas_desc%gas_name) - end do - icnt = count(gas_is_present) - allocate(get_minor_list(icnt)) - get_minor_list(:) = pack(this%gas_names, mask=gas_is_present) - end function get_minor_list - !-------------------------------------------------------------------------------------------------------------------- - ! ! Inquiry functions ! !-------------------------------------------------------------------------------------------------------------------- @@ -1170,7 +1284,7 @@ end function source_is_internal pure function source_is_external(this) class(ty_gas_optics_rrtmgp), intent(in) :: this logical :: source_is_external - source_is_external = allocated(this%solar_src) + source_is_external = allocated(this%solar_source) end function source_is_external !-------------------------------------------------------------------------------------------------------------------- @@ -1279,6 +1393,57 @@ function get_col_dry(vmr_h2o, plev, latitude) result(col_dry) end function get_col_dry !-------------------------------------------------------------------------------------------------------------------- ! + ! Compute a transport angle that minimizes flux errors at surface and TOA based on empirical fits + ! + function compute_optimal_angles(this, optical_props, optimal_angles) result(err_msg) + ! input + class(ty_gas_optics_rrtmgp), intent(in ) :: this + class(ty_optical_props_arry), intent(in ) :: optical_props + real(wp), dimension(:,:), intent( out) :: optimal_angles + character(len=128) :: err_msg + !---------------------------- + integer :: ncol, nlay, ngpt + integer :: icol, ilay, igpt, bnd + real(wp) :: t, trans_total + !---------------------------- + ncol = optical_props%get_ncol() + nlay = optical_props%get_nlay() + ngpt = optical_props%get_ngpt() + + err_msg="" + if(.not. this%gpoints_are_equal(optical_props)) & + err_msg = "gas_optics%compute_optimal_angles: optical_props has different spectral discretization than gas_optics" + if(.not. extents_are(optimal_angles, ncol, ngpt)) & + err_msg = "gas_optics%compute_optimal_angles: optimal_angles different dimension (ncol)" + if (err_msg /= "") return + + ! + ! column transmissivity + ! + !$acc parallel loop gang vector collapse(2) copyin(optical_props, optical_props%tau, optical_props%gpt2band) copyout(optimal_angles) + do icol = 1, ncol + do igpt = 1, ngpt + ! + ! Column transmissivity + ! + t = 0._wp + trans_total = 0._wp + do ilay = 1, nlay + t = t + optical_props%tau(icol,ilay,igpt) + end do + trans_total = exp(-t) + ! + ! Optimal transport angle is a linear fit to column transmissivity + ! + bnd = optical_props%gpt2band(igpt) + optimal_angles(icol,igpt) = this%optimal_angle_fit(1,bnd)*trans_total + & + this%optimal_angle_fit(2,bnd) + end do + end do + + end function compute_optimal_angles + !-------------------------------------------------------------------------------------------------------------------- + ! ! Internal procedures ! !-------------------------------------------------------------------------------------------------------------------- diff --git a/rte/kernels-openacc/mo_rte_solver_kernels.F90 b/rte/kernels-openacc/mo_rte_solver_kernels.F90 index 19384bf6a..f234e1500 100644 --- a/rte/kernels-openacc/mo_rte_solver_kernels.F90 +++ b/rte/kernels-openacc/mo_rte_solver_kernels.F90 @@ -61,7 +61,8 @@ module mo_rte_solver_kernels ! --------------------------------------------------------------- subroutine lw_solver_noscat(ncol, nlay, ngpt, top_at_1, D, weight, & tau, lay_source, lev_source_inc, lev_source_dec, sfc_emis, sfc_src, & - radn_up, radn_dn) bind(C, name="lw_solver_noscat") + radn_up, radn_dn, & + sfc_srcJac, radn_upJac) bind(C, name="lw_solver_noscat") integer, intent( in) :: ncol, nlay, ngpt ! Number of columns, layers, g-points logical(wl), intent( in) :: top_at_1 real(wp), dimension(ncol, ngpt), intent( in) :: D ! secant of propagation angle [] @@ -79,11 +80,13 @@ subroutine lw_solver_noscat(ncol, nlay, ngpt, top_at_1, D, weight, real(wp), dimension(ncol,nlay+1,ngpt), intent( out) :: radn_up ! Radiances [W/m2-str] ! Top level must contain incident flux boundary condition + real(wp), dimension(ncol ,ngpt), intent(in ) :: sfc_srcJac ! surface temperature Jacobian of surface source function [W/m2/K] + real(wp), dimension(ncol,nlay+1,ngpt), intent(out) :: radn_upJac ! surface temperature Jacobian of Radiances [W/m2-str / K] ! Local variables, no g-point dependency real(wp), dimension(ncol,nlay,ngpt) :: tau_loc, & ! path length (tau/mu) trans ! transmissivity = exp(-tau) real(wp), dimension(ncol,nlay,ngpt) :: source_dn, source_up - real(wp), dimension(ncol, ngpt) :: source_sfc, sfc_albedo + real(wp), dimension(ncol, ngpt) :: source_sfc, sfc_albedo, source_sfcJac real(wp), dimension(:,:,:), pointer :: lev_source_up, lev_source_dn ! Mapping increasing/decreasing indicies to up/down @@ -130,6 +133,14 @@ subroutine lw_solver_noscat(ncol, nlay, ngpt, top_at_1, D, weight, end do end do + !$acc enter data copyin(sfc_srcJac) + !$acc enter data create(source_sfcJac, radn_upJac) + !$acc parallel loop collapse(2) + do igpt = 1, ngpt + do icol = 1, ncol + source_sfcJac(icol,igpt) = sfc_emis(icol,igpt) * sfc_srcJac(icol,igpt) + end do + end do ! NOTE: This kernel produces small differences between GPU and CPU ! implementations on Ascent with PGI, we assume due to floating point ! differences in the exp() function. These differences are small in the @@ -157,7 +168,7 @@ subroutine lw_solver_noscat(ncol, nlay, ngpt, top_at_1, D, weight, ! call lw_transport_noscat(ncol, nlay, ngpt, top_at_1, & tau_loc, trans, sfc_albedo, source_dn, source_up, source_sfc, & - radn_up, radn_dn) + radn_up, radn_dn, source_sfcJac, radn_upJac) ! ! Convert intensity to flux assuming azimuthal isotropy and quadrature weight @@ -166,11 +177,14 @@ subroutine lw_solver_noscat(ncol, nlay, ngpt, top_at_1, D, weight, do igpt = 1, ngpt do ilev = 1, nlay+1 do icol = 1, ncol - radn_dn(icol,ilev,igpt) = 2._wp * pi * weight * radn_dn(icol,ilev,igpt) - radn_up(icol,ilev,igpt) = 2._wp * pi * weight * radn_up(icol,ilev,igpt) + radn_dn (icol,ilev,igpt) = 2._wp * pi * weight * radn_dn (icol,ilev,igpt) + radn_up (icol,ilev,igpt) = 2._wp * pi * weight * radn_up (icol,ilev,igpt) + radn_upJac(icol,ilev,igpt) = 2._wp * pi * weight * radn_upJac(icol,ilev,igpt) end do end do end do + !$acc exit data delete(sfc_srcJac, source_sfcJac) + !$acc exit data copyout(radn_upJac) !$acc exit data copyout(radn_dn,radn_up) !$acc exit data delete(d,tau,sfc_src,sfc_emis,lev_source_dec,lev_source_inc,lay_source,tau_loc,trans,source_dn,source_up,source_sfc,sfc_albedo) @@ -185,7 +199,8 @@ end subroutine lw_solver_noscat ! ! --------------------------------------------------------------- subroutine lw_solver_noscat_GaussQuad(ncol, nlay, ngpt, top_at_1, nmus, Ds, weights, & - tau, lay_source, lev_source_inc, lev_source_dec, sfc_emis, sfc_src, flux_up, flux_dn) & + tau, lay_source, lev_source_inc, lev_source_dec, sfc_emis, sfc_src, flux_up, flux_dn,& + sfc_srcJac, flux_upJac) & bind (C, name="lw_solver_noscat_GaussQuad") integer, intent(in ) :: ncol, nlay, ngpt ! Number of columns, layers, g-points logical(wl), intent(in ) :: top_at_1 @@ -203,8 +218,10 @@ subroutine lw_solver_noscat_GaussQuad(ncol, nlay, ngpt, top_at_1, nmus, Ds, weig real(wp), dimension(ncol,nlay+1,ngpt), intent(inout) :: flux_dn ! Radiances [W/m2-str] real(wp), dimension(ncol,nlay+1,ngpt), intent( out) :: flux_up ! Radiances [W/m2-str] ! Top level must contain incident flux boundary condition + real(wp), dimension(ncol ,ngpt), intent(in ) :: sfc_srcJac ! surface temperature Jacobian of surface source function [W/m2/K] + real(wp), dimension(ncol,nlay+1,ngpt), intent(out) :: flux_upJac ! surface temperature Jacobian of Radiances [W/m2-str / K] ! Local variables - real(wp), dimension(ncol,nlay+1,ngpt) :: radn_dn, radn_up ! Fluxes per quad angle + real(wp), dimension(ncol,nlay+1,ngpt) :: radn_dn, radn_up, radn_upJac ! Fluxes per quad angle real(wp), dimension(ncol, ngpt) :: Ds_ncol real(wp), dimension(ncol, ngpt) :: flux_top @@ -213,6 +230,8 @@ subroutine lw_solver_noscat_GaussQuad(ncol, nlay, ngpt, top_at_1, nmus, Ds, weig !$acc enter data copyin(Ds,weights,tau,lay_source,lev_source_inc,lev_source_dec,sfc_emis,sfc_src,flux_dn) !$acc enter data create(flux_up,radn_dn,radn_up,Ds_ncol,flux_top) + !$acc enter data copyin(sfc_srcJac) + !$acc enter data create(flux_upJac, radn_upJac) ! ------------------------------------ ! ------------------------------------ @@ -226,7 +245,7 @@ subroutine lw_solver_noscat_GaussQuad(ncol, nlay, ngpt, top_at_1, nmus, Ds, weig call lw_solver_noscat(ncol, nlay, ngpt, & top_at_1, Ds_ncol, weights(1), tau, & lay_source, lev_source_inc, lev_source_dec, sfc_emis, sfc_src, & - flux_up, flux_dn) + flux_up, flux_dn, sfc_srcJac, flux_upJac) ! ! For more than one angle use local arrays ! @@ -249,18 +268,21 @@ subroutine lw_solver_noscat_GaussQuad(ncol, nlay, ngpt, top_at_1, nmus, Ds, weig call lw_solver_noscat(ncol, nlay, ngpt, & top_at_1, Ds_ncol, weights(imu), tau, & lay_source, lev_source_inc, lev_source_dec, sfc_emis, sfc_src, & - radn_up, radn_dn) + radn_up, radn_dn, sfc_srcJac, radn_upJac) !$acc parallel loop collapse(3) do igpt = 1, ngpt do ilev = 1, nlay+1 do icol = 1, ncol - flux_up(icol,ilev,ngpt) = flux_up(icol,ilev,ngpt) + radn_up(icol,ilev,ngpt) - flux_dn(icol,ilev,ngpt) = flux_dn(icol,ilev,ngpt) + radn_dn(icol,ilev,ngpt) + flux_up (icol,ilev,igpt) = flux_up (icol,ilev,igpt) + radn_up (icol,ilev,igpt) + flux_dn (icol,ilev,igpt) = flux_dn (icol,ilev,igpt) + radn_dn (icol,ilev,igpt) + flux_upJac(icol,ilev,igpt) = flux_upJac(icol,ilev,igpt) + radn_upJac(icol,ilev,igpt) end do end do end do end do ! imu loop + !$acc exit data delete(sfc_srcJac, radn_upJac) + !$acc exit data copyout(flux_upJac) !$acc exit data copyout(flux_up,flux_dn) !$acc exit data delete(Ds,weights,tau,lay_source,lev_source_inc,lev_source_dec,sfc_emis,sfc_src,radn_dn,radn_up,Ds_ncol,flux_top) @@ -554,7 +576,7 @@ end subroutine lw_source_noscat ! --------------------------------------------------------------- subroutine lw_transport_noscat(ncol, nlay, ngpt, top_at_1, & tau, trans, sfc_albedo, source_dn, source_up, source_sfc, & - radn_up, radn_dn) bind(C, name="lw_transport_noscat") + radn_up, radn_dn, source_sfcJac, radn_upJac) bind(C, name="lw_transport_noscat") integer, intent(in ) :: ncol, nlay, ngpt ! Number of columns, layers, g-points logical(wl), intent(in ) :: top_at_1 ! real(wp), dimension(ncol,nlay ,ngpt), intent(in ) :: tau, & ! Absorption optical thickness, pre-divided by mu [] @@ -566,6 +588,8 @@ subroutine lw_transport_noscat(ncol, nlay, ngpt, top_at_1, & real(wp), dimension(ncol,nlay+1,ngpt), intent(inout) :: radn_dn ! Radiances [W/m2-str] real(wp), dimension(ncol,nlay+1,ngpt), intent( out) :: radn_up ! Radiances [W/m2-str] ! Top level must contain incident flux boundary condition + real(wp), dimension(ncol ,ngpt), intent(in ) :: source_sfcJac ! surface temperature Jacobian of surface source function [W/m2/K] + real(wp), dimension(ncol,nlay+1,ngpt), intent(out) :: radn_upJac ! surface temperature Jacobian of Radiances [W/m2-str / K] ! Local variables integer :: igpt, ilev, icol ! --------------------------------------------------- @@ -583,11 +607,13 @@ subroutine lw_transport_noscat(ncol, nlay, ngpt, top_at_1, & end do ! Surface reflection and emission - radn_up(icol,nlay+1,igpt) = radn_dn(icol,nlay+1,igpt)*sfc_albedo(icol,igpt) + source_sfc(icol,igpt) + radn_up (icol,nlay+1,igpt) = radn_dn(icol,nlay+1,igpt)*sfc_albedo(icol,igpt) + source_sfc (icol,igpt) + radn_upJac(icol,nlay+1,igpt) = source_sfcJac(icol,igpt) ! Upward propagation do ilev = nlay, 1, -1 - radn_up(icol,ilev,igpt) = trans(icol,ilev ,igpt)*radn_up(icol,ilev+1,igpt) + source_up(icol,ilev,igpt) + radn_up (icol,ilev,igpt) = trans(icol,ilev,igpt)*radn_up (icol,ilev+1,igpt) + source_up(icol,ilev,igpt) + radn_upJac(icol,ilev,igpt) = trans(icol,ilev,igpt)*radn_upJac(icol,ilev+1,igpt) end do end do end do @@ -604,11 +630,13 @@ subroutine lw_transport_noscat(ncol, nlay, ngpt, top_at_1, & end do ! Surface reflection and emission - radn_up(icol, 1,igpt) = radn_dn(icol, 1,igpt)*sfc_albedo(icol,igpt) + source_sfc(icol,igpt) + radn_up (icol,1,igpt) = radn_dn(icol,1,igpt)*sfc_albedo(icol,igpt) + source_sfc (icol,igpt) + radn_upJac(icol,1,igpt) = source_sfcJac(icol,igpt) ! Upward propagation do ilev = 2, nlay+1 - radn_up(icol,ilev,igpt) = trans(icol,ilev-1,igpt) * radn_up(icol,ilev-1,igpt) + source_up(icol,ilev-1,igpt) + radn_up (icol,ilev,igpt) = trans(icol,ilev-1,igpt) * radn_up (icol,ilev-1,igpt) + source_up(icol,ilev-1,igpt) + radn_upJac(icol,ilev,igpt) = trans(icol,ilev-1,igpt) * radn_upJac(icol,ilev-1,igpt) end do end do end do @@ -1216,7 +1244,8 @@ end subroutine apply_BC_0 ! ------------------------------------------------------------------------------------------------- subroutine lw_solver_1rescl(ncol, nlay, ngpt, top_at_1, D, & tau, scaling, lay_source, lev_source_inc, lev_source_dec, sfc_emis, sfc_src, & - radn_up, radn_dn) bind(C, name="lw_solver_1rescl") + radn_up, radn_dn, & + sfc_srcJac, rad_up_Jac, rad_dn_Jac) bind(C, name="lw_solver_1rescl") integer, intent(in ) :: ncol, nlay, ngpt ! Number of columns, layers, g-points logical(wl), intent(in ) :: top_at_1 real(wp), dimension(ncol, ngpt), intent(in ) :: D ! secant of propagation angle [] @@ -1233,11 +1262,15 @@ subroutine lw_solver_1rescl(ncol, nlay, ngpt, top_at_1, D, real(wp), dimension(ncol,nlay+1,ngpt), intent( out) :: radn_up ! Radiances [W/m2-str] real(wp), dimension(ncol,nlay+1,ngpt), intent(inout) :: radn_dn ! Top level must contain incident flux boundary condition + real(wp), dimension(ncol, ngpt), intent(in ) :: sfc_srcJac ! Surface Temperature Jacobian source function [W/m2/K] + real(wp), dimension(ncol,nlay+1,ngpt), intent( out) :: rad_up_Jac ! Surface Temperature Jacobians [W/m2-str/K] + real(wp), dimension(ncol,nlay+1,ngpt), intent( out) :: rad_dn_Jac ! Top level set to 0 ! Local variables, WITH g-point dependency real(wp), dimension(ncol,nlay,ngpt) :: tau_loc, & ! path length (tau/mu) trans ! transmissivity = exp(-tau) real(wp), dimension(ncol,nlay,ngpt) :: source_dn, source_up real(wp), dimension(ncol, ngpt) :: source_sfc, sfc_albedo + real(wp), dimension(ncol, ngpt) :: source_sfcJac real(wp), dimension(:,:,:), pointer :: lev_source_up, lev_source_dn ! Mapping increasing/decreasing indicies to up/down @@ -1271,6 +1304,9 @@ subroutine lw_solver_1rescl(ncol, nlay, ngpt, top_at_1, D, !$acc enter data create(sfcSource, An, Cn) !$acc enter data attach(lev_source_up,lev_source_dn) + !$acc enter data create(rad_up_Jac, rad_dn_Jac, source_sfcJac) + !$acc enter data copyin(sfc_srcJac) + ! NOTE: This kernel produces small differences between GPU and CPU ! implementations on Ascent with PGI, we assume due to floating point ! differences in the exp() function. These differences are small in the @@ -1288,6 +1324,9 @@ subroutine lw_solver_1rescl(ncol, nlay, ngpt, top_at_1, D, ! explanation of factor 0.4 note A of Table Cn(icol,ilev,igpt) = 0.4_wp*scaling(icol,ilev,igpt) An(icol,ilev,igpt) = (1._wp-trans(icol,ilev,igpt)*trans(icol,ilev,igpt)) + + ! initialize radn_dn_Jac + rad_dn_Jac(icol,ilev,igpt) = 0._wp end do end do end do @@ -1298,8 +1337,9 @@ subroutine lw_solver_1rescl(ncol, nlay, ngpt, top_at_1, D, ! ! Surface albedo, surface source function ! - sfc_albedo(icol,igpt) = 1._wp - sfc_emis(icol,igpt) - source_sfc(icol,igpt) = sfc_emis(icol,igpt) * sfc_src(icol,igpt) + sfc_albedo (icol,igpt) = 1._wp - sfc_emis(icol,igpt) + source_sfc (icol,igpt) = sfc_emis(icol,igpt) * sfc_src (icol,igpt) + source_sfcJac(icol,igpt) = sfc_emis(icol,igpt) * sfc_srcJac(icol,igpt) end do end do @@ -1317,12 +1357,17 @@ subroutine lw_solver_1rescl(ncol, nlay, ngpt, top_at_1, D, ! compute no-scattering fluxes call lw_transport_noscat(ncol, nlay, ngpt, top_at_1, & tau_loc, trans, sfc_albedo, source_dn, source_up, source_sfc, & - radn_up, radn_dn) + radn_up, radn_dn,& + source_sfcJac, rad_up_Jac) ! make adjustment call lw_transport_1rescl(ncol, nlay, ngpt, top_at_1, & tau_loc, trans, & - sfc_albedo, source_dn, source_up, source_sfc, & - radn_up, radn_dn, An, Cn) + sfc_albedo, source_dn, source_up, & + radn_up, radn_dn, An, Cn, rad_up_Jac, rad_dn_Jac) + + !$acc exit data copyout(rad_up_Jac, rad_dn_Jac) + !$acc exit data delete(source_sfcJac, sfc_srcJac) + !$acc exit data copyout(radn_dn,radn_up) !$acc exit data delete(sfcSource, An, Cn) @@ -1342,7 +1387,8 @@ end subroutine lw_solver_1rescl subroutine lw_solver_1rescl_GaussQuad(ncol, nlay, ngpt, top_at_1, nmus, Ds, weights, & tau, ssa, g, lay_source, lev_source_inc, lev_source_dec, & sfc_emis, sfc_src,& - flux_up, flux_dn) & + flux_up, flux_dn, & + sfc_src_Jac, flux_up_Jac, flux_dn_Jac) & bind(C, name="lw_solver_1rescl_GaussQuad") integer, intent(in ) :: ncol, nlay, ngpt ! Number of columns, layers, g-points logical(wl), intent(in ) :: top_at_1 @@ -1352,18 +1398,24 @@ subroutine lw_solver_1rescl_GaussQuad(ncol, nlay, ngpt, top_at_1, nmus, Ds, weig real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: ssa ! single-scattering albedo, real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: g ! asymmetry parameter [] real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: lay_source ! Planck source at layer average temperature [W/m2] - real(wp), dimension(ncol,nlay+1,ngpt), intent(in ) :: lev_source_inc + real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: lev_source_inc ! Planck source at layer edge for radiation in increasing ilay direction [W/m2] ! Includes spectral weighting that accounts for state-dependent frequency to g-space mapping - real(wp), dimension(ncol,nlay+1,ngpt), intent(in ) :: lev_source_dec + real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: lev_source_dec ! Planck source at layer edge for radiation in decreasing ilay direction [W/m2] real(wp), dimension(ncol, ngpt), intent(in ) :: sfc_emis ! Surface emissivity [] real(wp), dimension(ncol, ngpt), intent(in ) :: sfc_src ! Surface source function [W/m2] real(wp), dimension(ncol,nlay+1,ngpt), intent( out) :: flux_up ! Radiances [W/m2-str] real(wp), dimension(ncol,nlay+1,ngpt), intent(inout) :: flux_dn ! Top level must contain incident flux boundary condition + + real(wp), dimension(ncol ,ngpt), intent(in ) :: sfc_src_Jac ! surface temperature Jacobian of surface source function [W/m2/K] + real(wp), dimension(ncol,nlay+1,ngpt), intent(out) :: flux_up_Jac ! surface temperature Jacobian of Radiances [W/m2-str / K] + real(wp), dimension(ncol,nlay+1,ngpt), intent(out) :: flux_dn_Jac ! surface temperature Jacobian of Radiances [W/m2-str / K] + ! Local variables real(wp), dimension(ncol,nlay+1,ngpt) :: radn_dn, radn_up ! Fluxes per quad angle real(wp), dimension(ncol, ngpt) :: Ds_ncol + real(wp), dimension(ncol,nlay+1,ngpt) :: radn_dn_Jac, radn_up_Jac ! Surface temperature Jsacobians per quad angle integer :: imu, top_level,icol,ilev,igpt real :: weight @@ -1373,8 +1425,8 @@ subroutine lw_solver_1rescl_GaussQuad(ncol, nlay, ngpt, top_at_1, nmus, Ds, weig real(wp), dimension(ncol,nlay, ngpt) :: scaling ! scaling real(wp), parameter :: tresh=1.0_wp - 1e-6_wp - !$acc enter data copyin(Ds,weights,tau,ssa,g,lay_source,lev_source_inc,lev_source_dec,sfc_emis,sfc_src,flux_dn) - !$acc enter data create(flux_up,radn_dn,radn_up,Ds_ncol, scaling, tauLoc) + !$acc enter data copyin(Ds,weights,tau,ssa,g,lay_source,lev_source_inc,lev_source_dec,sfc_emis,sfc_src,flux_dn,sfc_src_Jac) + !$acc enter data create(flux_up,radn_dn,radn_up,Ds_ncol, scaling, tauLoc,flux_up_Jac,flux_dn_Jac,radn_dn_Jac, radn_up_Jac) ! Tang rescaling @@ -1397,18 +1449,20 @@ subroutine lw_solver_1rescl_GaussQuad(ncol, nlay, ngpt, top_at_1, nmus, Ds, weig ! Transport is for intensity ! convert flux at top of domain to intensity assuming azimuthal isotropy ! - radn_dn(1:ncol, top_level, 1:ngpt) = fluxTOA(1:ncol, 1:ngpt) / weight + radn_dn(1:ncol, top_level, 1:ngpt) = fluxTOA(1:ncol, 1:ngpt) / weight call lw_solver_1rescl(ncol, nlay, ngpt, & top_at_1, Ds_ncol, tauLoc, scaling, & lay_source, lev_source_inc, lev_source_dec, sfc_emis, sfc_src, & - flux_up, flux_dn) + flux_up, flux_dn,sfc_src_Jac,flux_up_Jac,flux_dn_Jac) !$acc parallel loop collapse(3) do igpt = 1, ngpt do ilev = 1, nlay+1 do icol = 1, ncol - flux_up(icol,ilev,igpt) = weight*flux_up(icol,ilev,igpt) - flux_dn(icol,ilev,igpt) = weight*flux_dn(icol,ilev,igpt) + flux_up (icol,ilev,igpt) = weight*flux_up (icol,ilev,igpt) + flux_dn (icol,ilev,igpt) = weight*flux_dn (icol,ilev,igpt) + flux_up_Jac(icol,ilev,igpt) = weight*flux_up_Jac(icol,ilev,igpt) + flux_dn_Jac(icol,ilev,igpt) = weight*flux_dn_Jac(icol,ilev,igpt) enddo enddo enddo @@ -1420,19 +1474,23 @@ subroutine lw_solver_1rescl_GaussQuad(ncol, nlay, ngpt, top_at_1, nmus, Ds, weig call lw_solver_1rescl(ncol, nlay, ngpt, & top_at_1, Ds_ncol, tauLoc, scaling, & lay_source, lev_source_inc, lev_source_dec, sfc_emis, sfc_src, & - radn_up, radn_dn) + radn_up, radn_dn,sfc_src_Jac,radn_up_Jac,radn_dn_Jac) !$acc parallel loop collapse(3) do igpt = 1, ngpt do ilev = 1, nlay+1 do icol = 1, ncol - flux_up(icol,ilev,igpt) = flux_up(icol,ilev,igpt) + weight*radn_up(icol,ilev,igpt) - flux_dn(icol,ilev,igpt) = flux_dn(icol,ilev,igpt) + weight*radn_dn(icol,ilev,igpt) + flux_up (icol,ilev,igpt) = flux_up (icol,ilev,igpt) + weight*radn_up (icol,ilev,igpt) + flux_dn (icol,ilev,igpt) = flux_dn (icol,ilev,igpt) + weight*radn_dn (icol,ilev,igpt) + flux_up_Jac(icol,ilev,igpt) = flux_up_Jac(icol,ilev,igpt) + weight*radn_up_Jac(icol,ilev,igpt) + flux_dn_Jac(icol,ilev,igpt) = flux_dn_Jac(icol,ilev,igpt) + weight*radn_dn_Jac(icol,ilev,igpt) enddo enddo enddo end do + !$acc exit data delete(sfc_src_Jac,radn_dn_Jac, radn_up_Jac) + !$acc exit data copyout(flux_up_Jac,flux_dn_Jac) !$acc exit data copyout(flux_up,flux_dn) !$acc exit data delete(Ds,weights,tau,ssa,g,tauLoc,scaling,lay_source,lev_source_inc,lev_source_dec,sfc_emis,sfc_src,radn_dn,radn_up,Ds_ncol) end subroutine lw_solver_1rescl_GaussQuad @@ -1532,19 +1590,21 @@ end subroutine scaling_1rescl_safe ! ! ------------------------------------------------------------------------------------------------- subroutine lw_transport_1rescl(ncol, nlay, ngpt, top_at_1, & - tau, trans, sfc_albedo, source_dn, source_up, source_sfc, & - radn_up, radn_dn, An, Cn) bind(C, name="lw_transport_1rescl") + tau, trans, sfc_albedo, source_dn, source_up, & + radn_up, radn_dn, An, Cn,& + rad_up_Jac, rad_dn_Jac) bind(C, name="lw_transport_1rescl") integer, intent(in ) :: ncol, nlay, ngpt ! Number of columns, layers, g-points logical(wl), intent(in ) :: top_at_1 ! real(wp), dimension(ncol,nlay ,ngpt), intent(in ) :: tau, & ! Absorption optical thickness, pre-divided by mu [] trans ! transmissivity = exp(-tau) real(wp), dimension(ncol ,ngpt), intent(in ) :: sfc_albedo ! Surface albedo real(wp), dimension(ncol,nlay ,ngpt), intent(in ) :: source_dn, & - source_up ! Diffuse radiation emitted by the layer - real(wp), dimension(ncol ,ngpt), intent(in ) :: source_sfc ! Surface source function [W/m2] + source_up ! Diffuse radiation emitted by the layer real(wp), dimension(ncol,nlay+1,ngpt), intent(inout) :: radn_up ! Radiances [W/m2-str] real(wp), dimension(ncol,nlay+1,ngpt), intent(inout) :: radn_dn !Top level must contain incident flux boundary condition - real(wp), dimension(ncol,nlay ,ngpt), intent(in ) :: An, Cn + real(wp), dimension(ncol,nlay ,ngpt), intent(in ) :: An, Cn + real(wp), dimension(ncol,nlay+1,ngpt), intent(inout) :: rad_up_Jac ! Radiances [W/m2-str] + real(wp), dimension(ncol,nlay+1,ngpt), intent(inout) :: rad_dn_Jac !Top level must contain incident flux boundary condition ! Local variables integer :: ilev, icol, igpt ! --------------------------------------------------- @@ -1559,7 +1619,9 @@ subroutine lw_transport_1rescl(ncol, nlay, ngpt, top_at_1, & do icol = 1, ncol ! 1st Upward propagation do ilev = nlay, 1, -1 - radn_up(icol,ilev,igpt) = trans(icol,ilev,igpt)*radn_up(icol,ilev+1,igpt) + source_up(icol,ilev,igpt) + radn_up (icol,ilev,igpt) = trans(icol,ilev,igpt)*radn_up (icol,ilev+1,igpt) + source_up(icol,ilev,igpt) + rad_up_Jac(icol,ilev,igpt) = trans(icol,ilev,igpt)*rad_up_Jac(icol,ilev+1,igpt) + adjustmentFactor = Cn(icol,ilev,igpt)*& ( An(icol,ilev,igpt)*radn_dn(icol,ilev,igpt) - & source_dn(icol,ilev,igpt) *trans(icol,ilev,igpt ) - & @@ -1568,12 +1630,16 @@ subroutine lw_transport_1rescl(ncol, nlay, ngpt, top_at_1, & enddo ! 2nd Downward propagation do ilev = 1, nlay - radn_dn(icol,ilev+1,igpt) = trans(icol,ilev,igpt)*radn_dn(icol,ilev,igpt) + source_dn(icol,ilev,igpt) + radn_dn (icol,ilev+1,igpt) = trans(icol,ilev,igpt)*radn_dn (icol,ilev,igpt) + source_dn(icol,ilev,igpt) + rad_dn_Jac(icol,ilev+1,igpt) = trans(icol,ilev,igpt)*rad_dn_Jac(icol,ilev,igpt) adjustmentFactor = Cn(icol,ilev,igpt)*( & An(icol,ilev,igpt)*radn_up(icol,ilev,igpt) - & source_up(icol,ilev,igpt)*trans(icol,ilev,igpt) - & source_dn(icol,ilev,igpt) ) - radn_dn(icol,ilev+1,igpt) = radn_dn(icol,ilev+1,igpt) + adjustmentFactor + radn_dn(icol,ilev+1,igpt) = radn_dn(icol,ilev+1,igpt) + adjustmentFactor + + adjustmentFactor = Cn(icol,ilev,igpt)*An(icol,ilev,igpt)*rad_up_Jac(icol,ilev,igpt) + rad_dn_Jac(icol,ilev+1,igpt) = rad_dn_Jac(icol,ilev+1,igpt) + adjustmentFactor enddo enddo enddo @@ -1583,7 +1649,8 @@ subroutine lw_transport_1rescl(ncol, nlay, ngpt, top_at_1, & do icol = 1, ncol ! Upward propagation do ilev = 1, nlay - radn_up(icol,ilev+1,igpt) = trans(icol,ilev,igpt) * radn_up(icol,ilev,igpt) + source_up(icol,ilev,igpt) + radn_up (icol,ilev+1,igpt) = trans(icol,ilev,igpt)*radn_up (icol,ilev,igpt) + source_up(icol,ilev,igpt) + rad_up_Jac(icol,ilev+1,igpt) = trans(icol,ilev,igpt)*rad_up_Jac(icol,ilev,igpt) adjustmentFactor = Cn(icol,ilev,igpt)*& ( An(icol,ilev,igpt)*radn_dn(icol,ilev+1,igpt) - & source_dn(icol,ilev,igpt) *trans(icol,ilev ,igpt) - & @@ -1592,12 +1659,16 @@ subroutine lw_transport_1rescl(ncol, nlay, ngpt, top_at_1, & end do ! 2st Downward propagation do ilev = nlay, 1, -1 - radn_dn(icol,ilev,igpt) = trans(icol,ilev,igpt)*radn_dn(icol,ilev+1,igpt) + source_dn(icol,ilev,igpt) + radn_dn (icol,ilev,igpt) = trans(icol,ilev,igpt)*radn_dn (icol,ilev+1,igpt) + source_dn(icol,ilev,igpt) + rad_dn_Jac(icol,ilev,igpt) = trans(icol,ilev,igpt)*rad_dn_Jac(icol,ilev+1,igpt) adjustmentFactor = Cn(icol,ilev,igpt)*( & An(icol,ilev,igpt)*radn_up(icol,ilev,igpt) - & source_up(icol,ilev,igpt)*trans(icol,ilev ,igpt ) - & source_dn(icol,ilev,igpt) ) - radn_dn(icol,ilev,igpt) = radn_dn(icol,ilev,igpt) + adjustmentFactor + radn_dn(icol,ilev,igpt) = radn_dn(icol,ilev,igpt) + adjustmentFactor + + adjustmentFactor = Cn(icol,ilev,igpt)*An(icol,ilev,igpt)*rad_up_Jac(icol,ilev,igpt) + rad_dn_Jac(icol,ilev,igpt) = rad_dn_Jac(icol,ilev,igpt) + adjustmentFactor end do enddo enddo diff --git a/rte/kernels/mo_rte_solver_kernels.F90 b/rte/kernels/mo_rte_solver_kernels.F90 index 156ff027a..7c8604472 100644 --- a/rte/kernels/mo_rte_solver_kernels.F90 +++ b/rte/kernels/mo_rte_solver_kernels.F90 @@ -62,7 +62,8 @@ module mo_rte_solver_kernels ! --------------------------------------------------------------- subroutine lw_solver_noscat(ncol, nlay, ngpt, top_at_1, D, weight, & tau, lay_source, lev_source_inc, lev_source_dec, sfc_emis, sfc_src, & - radn_up, radn_dn) bind(C, name="lw_solver_noscat") + radn_up, radn_dn, & + sfc_srcJac, radn_upJac) bind(C, name="lw_solver_noscat") integer, intent(in ) :: ncol, nlay, ngpt ! Number of columns, layers, g-points logical(wl), intent(in ) :: top_at_1 real(wp), dimension(ncol, ngpt), intent(in ) :: D ! secant of propagation angle [] @@ -79,11 +80,14 @@ subroutine lw_solver_noscat(ncol, nlay, ngpt, top_at_1, D, weight, real(wp), dimension(ncol,nlay+1,ngpt), intent( out) :: radn_up ! Radiances [W/m2-str] real(wp), dimension(ncol,nlay+1,ngpt), intent(inout) :: radn_dn ! Top level must contain incident flux boundary condition + real(wp), dimension(ncol ,ngpt), intent(in ) :: sfc_srcJac ! surface temperature Jacobian of surface source function [W/m2/K] + real(wp), dimension(ncol,nlay+1,ngpt), intent(out) :: radn_upJac ! surface temperature Jacobian of Radiances [W/m2-str / K] + ! Local variables, no g-point dependency real(wp), dimension(ncol,nlay) :: tau_loc, & ! path length (tau/mu) trans ! transmissivity = exp(-tau) real(wp), dimension(ncol,nlay) :: source_dn, source_up - real(wp), dimension(ncol ) :: source_sfc, sfc_albedo + real(wp), dimension(ncol ) :: source_sfc, sfc_albedo, source_sfcJac real(wp), dimension(:,:,:), pointer :: lev_source_up, lev_source_dn ! Mapping increasing/decreasing indicies to up/down @@ -129,17 +133,20 @@ subroutine lw_solver_noscat(ncol, nlay, ngpt, top_at_1, D, weight, ! sfc_albedo(:) = 1._wp - sfc_emis(:,igpt) source_sfc(:) = sfc_emis(:,igpt) * sfc_src(:,igpt) + source_sfcJac(:) = sfc_emis(:,igpt) * sfc_srcJac(:,igpt) ! ! Transport ! call lw_transport_noscat(ncol, nlay, top_at_1, & tau_loc, trans, sfc_albedo, source_dn, source_up, source_sfc, & - radn_up(:,:,igpt), radn_dn(:,:,igpt)) + radn_up(:,:,igpt), radn_dn(:,:,igpt), & + source_sfcJac, radn_upJac(:,:,igpt)) ! ! Convert intensity to flux assuming azimuthal isotropy and quadrature weight ! - radn_dn(:,:,igpt) = 2._wp * pi * weight * radn_dn(:,:,igpt) - radn_up(:,:,igpt) = 2._wp * pi * weight * radn_up(:,:,igpt) + radn_dn (:,:,igpt) = 2._wp * pi * weight * radn_dn (:,:,igpt) + radn_up (:,:,igpt) = 2._wp * pi * weight * radn_up (:,:,igpt) + radn_upJac(:,:,igpt) = 2._wp * pi * weight * radn_upJac(:,:,igpt) end do ! g point loop end subroutine lw_solver_noscat @@ -151,7 +158,8 @@ end subroutine lw_solver_noscat ! ! --------------------------------------------------------------- subroutine lw_solver_noscat_GaussQuad(ncol, nlay, ngpt, top_at_1, nmus, Ds, weights, & - tau, lay_source, lev_source_inc, lev_source_dec, sfc_emis, sfc_src, flux_up, flux_dn) & + tau, lay_source, lev_source_inc, lev_source_dec, sfc_emis, sfc_src, flux_up, flux_dn,& + sfc_srcJac, flux_upJac) & bind(C, name="lw_solver_noscat_GaussQuad") integer, intent(in ) :: ncol, nlay, ngpt ! Number of columns, layers, g-points logical(wl), intent(in ) :: top_at_1 @@ -168,9 +176,13 @@ subroutine lw_solver_noscat_GaussQuad(ncol, nlay, ngpt, top_at_1, nmus, Ds, weig real(wp), dimension(ncol, ngpt), intent(in ) :: sfc_src ! Surface source function [W/m2] real(wp), dimension(ncol,nlay+1,ngpt), intent( out) :: flux_up ! Radiances [W/m2-str] real(wp), dimension(ncol,nlay+1,ngpt), intent(inout) :: flux_dn ! Top level must contain incident flux boundary condition + + real(wp), dimension(ncol ,ngpt), intent(in ) :: sfc_srcJac ! surface temperature Jacobian of surface source function [W/m2/K] + real(wp), dimension(ncol,nlay+1,ngpt), intent(out) :: flux_upJac ! surface temperature Jacobian of Radiances [W/m2-str / K] ! Local variables real(wp), dimension(ncol,nlay+1,ngpt) :: radn_dn, radn_up ! Fluxes per quad angle real(wp), dimension(ncol, ngpt) :: Ds_ncol + real(wp), dimension(ncol,nlay+1,ngpt) :: radn_upJac ! perturbed Fluxes per quad angle integer :: imu, top_level ! ------------------------------------ @@ -181,7 +193,8 @@ subroutine lw_solver_noscat_GaussQuad(ncol, nlay, ngpt, top_at_1, nmus, Ds, weig call lw_solver_noscat(ncol, nlay, ngpt, & top_at_1, Ds_ncol, weights(1), tau, & lay_source, lev_source_inc, lev_source_dec, sfc_emis, sfc_src, & - flux_up, flux_dn) + flux_up, flux_dn, & + sfc_srcJac, flux_upJac) ! ! For more than one angle use local arrays ! @@ -193,9 +206,12 @@ subroutine lw_solver_noscat_GaussQuad(ncol, nlay, ngpt, top_at_1, nmus, Ds, weig call lw_solver_noscat(ncol, nlay, ngpt, & top_at_1, Ds_ncol, weights(imu), tau, & lay_source, lev_source_inc, lev_source_dec, sfc_emis, sfc_src, & - radn_up, radn_dn) - flux_up(:,:,:) = flux_up(:,:,:) + radn_up(:,:,:) - flux_dn(:,:,:) = flux_dn(:,:,:) + radn_dn(:,:,:) + radn_up, radn_dn, & + sfc_srcJac, radn_upJac) + flux_up (:,:,:) = flux_up (:,:,:) + radn_up (:,:,:) + flux_dn (:,:,:) = flux_dn (:,:,:) + radn_dn (:,:,:) + flux_upJac(:,:,:) = flux_upJac(:,:,:) + radn_upJac(:,:,:) + end do end subroutine lw_solver_noscat_GaussQuad ! ------------------------------------------------------------------------------------------------- @@ -427,7 +443,8 @@ end subroutine lw_source_noscat ! ------------------------------------------------------------------------------------------------- subroutine lw_transport_noscat(ncol, nlay, top_at_1, & tau, trans, sfc_albedo, source_dn, source_up, source_sfc, & - radn_up, radn_dn) bind(C, name="lw_transport_noscat") + radn_up, radn_dn,& + source_sfcJac, radn_upJac) bind(C, name="lw_transport_noscat") integer, intent(in ) :: ncol, nlay ! Number of columns, layers, g-points logical(wl), intent(in ) :: top_at_1 ! real(wp), dimension(ncol,nlay ), intent(in ) :: tau, & ! Absorption optical thickness, pre-divided by mu [] @@ -438,6 +455,9 @@ subroutine lw_transport_noscat(ncol, nlay, top_at_1, & real(wp), dimension(ncol ), intent(in ) :: source_sfc ! Surface source function [W/m2] real(wp), dimension(ncol,nlay+1), intent( out) :: radn_up ! Radiances [W/m2-str] real(wp), dimension(ncol,nlay+1), intent(inout) :: radn_dn !Top level must contain incident flux boundary condition + + real(wp), dimension(ncol ), intent(in ) :: source_sfcJac ! surface temperature Jacobian of surface source function [W/m2/K] + real(wp), dimension(ncol,nlay+1), intent(out) :: radn_upJac ! surface temperature Jacobian of Radiances [W/m2-str / K] ! Local variables integer :: ilev ! --------------------------------------------------- @@ -451,11 +471,13 @@ subroutine lw_transport_noscat(ncol, nlay, top_at_1, & end do ! Surface reflection and emission - radn_up(:,nlay+1) = radn_dn(:,nlay+1)*sfc_albedo(:) + source_sfc(:) + radn_up (:,nlay+1) = radn_dn(:,nlay+1)*sfc_albedo(:) + source_sfc (:) + radn_upJac(:,nlay+1) = source_sfcJac(:) ! Upward propagation do ilev = nlay, 1, -1 - radn_up(:,ilev) = trans(:,ilev )*radn_up(:,ilev+1) + source_up(:,ilev) + radn_up (:,ilev) = trans(:,ilev )*radn_up (:,ilev+1) + source_up(:,ilev) + radn_upJac(:,ilev) = trans(:,ilev )*radn_upJac(:,ilev+1) end do else ! @@ -467,11 +489,13 @@ subroutine lw_transport_noscat(ncol, nlay, top_at_1, & end do ! Surface reflection and emission - radn_up(:, 1) = radn_dn(:, 1)*sfc_albedo(:) + source_sfc(:) + radn_up (:, 1) = radn_dn(:,1)*sfc_albedo(:) + source_sfc (:) + radn_upJac(:, 1) = source_sfcJac(:) ! Upward propagation do ilev = 2, nlay+1 - radn_up(:,ilev) = trans(:,ilev-1) * radn_up(:,ilev-1) + source_up(:,ilev-1) + radn_up (:,ilev) = trans(:,ilev-1) * radn_up (:,ilev-1) + source_up(:,ilev-1) + radn_upJac(:,ilev) = trans(:,ilev-1) * radn_upJac(:,ilev-1) end do end if end subroutine lw_transport_noscat @@ -967,7 +991,8 @@ end subroutine apply_BC_0 ! ------------------------------------------------------------------------------------------------- subroutine lw_solver_1rescl(ncol, nlay, ngpt, top_at_1, D, & tau, scaling, lay_source, lev_source_inc, lev_source_dec, sfc_emis, sfc_src, & - radn_up, radn_dn) bind(C, name="lw_solver_1rescl") + radn_up, radn_dn, & + sfc_srcJac, rad_up_Jac, rad_dn_Jac) bind(C, name="lw_solver_1rescl") integer, intent(in ) :: ncol, nlay, ngpt ! Number of columns, layers, g-points logical(wl), intent(in ) :: top_at_1 real(wp), dimension(ncol, ngpt), intent(in ) :: D ! secant of propagation angle [] @@ -984,12 +1009,21 @@ subroutine lw_solver_1rescl(ncol, nlay, ngpt, top_at_1, D, real(wp), dimension(ncol,nlay+1,ngpt), intent( out) :: radn_up ! Radiances [W/m2-str] real(wp), dimension(ncol,nlay+1,ngpt), intent(inout) :: radn_dn ! Top level must contain incident flux boundary condition + real(wp), dimension(ncol, ngpt), intent(in ) :: sfc_srcJac ! Surface Temperature Jacobian source function [W/m2/K] + real(wp), dimension(ncol,nlay+1,ngpt), intent( out) :: rad_up_Jac ! Surface Temperature Jacobians [W/m2-str/K] + real(wp), dimension(ncol,nlay+1,ngpt), intent( out) :: rad_dn_Jac ! Top level set to 0 + ! Local variables, no g-point dependency real(wp), dimension(ncol,nlay) :: tau_loc, & ! path length (tau/mu) trans ! transmissivity = exp(-tau) real(wp), dimension(ncol,nlay) :: source_dn, source_up real(wp), dimension(ncol ) :: source_sfc, sfc_albedo real(wp), dimension(ncol,nlay) :: An, Cn + real(wp), dimension(ncol,nlay+1) :: dummyRadn_upJac + + real(wp), dimension(ncol,nlay+1) :: Radn_upJac + real(wp), dimension(ncol,nlay+1) :: Radn_dnJac + real(wp), dimension(ncol) :: source_sfcJac real(wp), dimension(:,:,:), pointer :: lev_source_up, lev_source_dn ! Mapping increasing/decreasing indicies to up/down @@ -1037,20 +1071,24 @@ subroutine lw_solver_1rescl(ncol, nlay, ngpt, top_at_1, D, ! ! Surface albedo, surface source function ! - sfc_albedo(:) = 1._wp - sfc_emis(:,igpt) - source_sfc(:) = sfc_emis(:,igpt) * sfc_src(:,igpt) + sfc_albedo(:) = 1._wp - sfc_emis(:,igpt) + source_sfc(:) = sfc_emis(:,igpt) * sfc_src (:,igpt) + source_sfcJac(:) = sfc_emis(:,igpt) * sfc_srcJac(:,igpt) ! ! Transport ! ! compute no-scattering fluxes call lw_transport_noscat(ncol, nlay, top_at_1, & tau_loc, trans, sfc_albedo, source_dn, source_up, source_sfc, & - radn_up(:,:,igpt), radn_dn(:,:,igpt)) + radn_up(:,:,igpt), radn_dn(:,:,igpt), & + source_sfcJac, rad_up_Jac(:,:,igpt)) + rad_dn_Jac(:,:,igpt) = 0._wp ! make adjustment call lw_transport_1rescl(ncol, nlay, top_at_1, trans, & - sfc_albedo, source_dn, source_up, source_sfc, & - radn_up(:,:,igpt), radn_dn(:,:,igpt), An, Cn) + source_dn, source_up, & + radn_up(:,:,igpt), radn_dn(:,:,igpt), An, Cn,& + rad_up_Jac(:,:,igpt), rad_dn_Jac(:,:,igpt)) end do ! g point loop end subroutine lw_solver_1rescl @@ -1064,8 +1102,9 @@ end subroutine lw_solver_1rescl ! --------------------------------------------------------------- subroutine lw_solver_1rescl_GaussQuad(ncol, nlay, ngpt, top_at_1, nmus, Ds, weights, & tau, ssa, g, lay_source, lev_source_inc, lev_source_dec, & - sfc_emis, sfc_src,& - flux_up, flux_dn) & + sfc_emis, sfc_src, & + flux_up, flux_dn, & + sfc_src_Jac, flux_up_Jac, flux_dn_Jac) & bind(C, name="lw_solver_1rescl_GaussQuad") integer, intent(in ) :: ncol, nlay, ngpt ! Number of columns, layers, g-points logical(wl), intent(in ) :: top_at_1 @@ -1075,18 +1114,23 @@ subroutine lw_solver_1rescl_GaussQuad(ncol, nlay, ngpt, top_at_1, nmus, Ds, weig real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: ssa ! single-scattering albedo, real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: g ! asymmetry parameter [] real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: lay_source ! Planck source at layer average temperature [W/m2] - real(wp), dimension(ncol,nlay+1,ngpt), intent(in ) :: lev_source_inc + real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: lev_source_inc ! Planck source at layer edge for radiation in increasing ilay direction [W/m2] ! Includes spectral weighting that accounts for state-dependent frequency to g-space mapping - real(wp), dimension(ncol,nlay+1,ngpt), intent(in ) :: lev_source_dec + real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: lev_source_dec ! Planck source at layer edge for radiation in decreasing ilay direction [W/m2] real(wp), dimension(ncol, ngpt), intent(in ) :: sfc_emis ! Surface emissivity [] real(wp), dimension(ncol, ngpt), intent(in ) :: sfc_src ! Surface source function [W/m2] real(wp), dimension(ncol,nlay+1,ngpt), intent( out) :: flux_up ! Radiances [W/m2-str] real(wp), dimension(ncol,nlay+1,ngpt), intent(inout) :: flux_dn ! Top level must contain incident flux boundary condition + + real(wp), dimension(ncol ,ngpt), intent(in ) :: sfc_src_Jac ! surface temperature Jacobian of surface source function [W/m2/K] + real(wp), dimension(ncol,nlay+1,ngpt), intent(out) :: flux_up_Jac ! surface temperature Jacobian of Radiances [W/m2-str / K] + real(wp), dimension(ncol,nlay+1,ngpt), intent(out) :: flux_dn_Jac ! surface temperature Jacobian of Radiances [W/m2-str / K] ! Local variables real(wp), dimension(ncol,nlay+1,ngpt) :: radn_dn, radn_up ! Fluxes per quad angle real(wp), dimension(ncol, ngpt) :: Ds_ncol + real(wp), dimension(ncol,nlay+1,ngpt) :: radn_dn_Jac, radn_up_Jac ! Fluxes per quad angle real(wp), dimension(ncol,nlay, ngpt) :: tauLoc ! rescaled Tau real(wp), dimension(ncol,nlay, ngpt) :: scaling ! scaling @@ -1114,15 +1158,15 @@ subroutine lw_solver_1rescl_GaussQuad(ncol, nlay, ngpt, top_at_1, nmus, Ds, weig ! convert flux at top of domain to intensity assuming azimuthal isotropy ! radn_dn(1:ncol, top_level, 1:ngpt) = fluxTOA(1:ncol, 1:ngpt) / weight - call lw_solver_1rescl(ncol, nlay, ngpt, & top_at_1, Ds_ncol, tauLoc, scaling, & lay_source, lev_source_inc, lev_source_dec, sfc_emis, sfc_src, & - flux_up, flux_dn) - - flux_up = flux_up * weight - flux_dn = flux_dn * weight + flux_up, flux_dn, sfc_src_Jac, flux_up_Jac, flux_dn_Jac) + flux_up = flux_up * weight + flux_dn = flux_dn * weight + flux_up_Jac = flux_up_Jac * weight + flux_dn_Jac = flux_dn_Jac * weight do imu = 2, nmus Ds_ncol(:,:) = Ds(imu) weight = 2._wp*pi*weights(imu) @@ -1133,10 +1177,12 @@ subroutine lw_solver_1rescl_GaussQuad(ncol, nlay, ngpt, top_at_1, nmus, Ds, weig call lw_solver_1rescl(ncol, nlay, ngpt, & top_at_1, Ds_ncol, tauLoc, scaling, & lay_source, lev_source_inc, lev_source_dec, sfc_emis, sfc_src, & - radn_up, radn_dn) + radn_up, radn_dn,sfc_src_Jac,radn_up_Jac,radn_dn_Jac) - flux_up(:,:,:) = flux_up(:,:,:) + weight*radn_up(:,:,:) - flux_dn(:,:,:) = flux_dn(:,:,:) + weight*radn_dn(:,:,:) + flux_up (:,:,:) = flux_up (:,:,:) + weight*radn_up (:,:,:) + flux_dn (:,:,:) = flux_dn (:,:,:) + weight*radn_dn (:,:,:) + flux_up_Jac(:,:,:) = flux_up_Jac(:,:,:) + weight*radn_up_Jac(:,:,:) + flux_dn_Jac(:,:,:) = flux_dn_Jac(:,:,:) + weight*radn_dn_Jac(:,:,:) end do end subroutine lw_solver_1rescl_GaussQuad ! ------------------------------------------------------------------------------------------------- @@ -1156,6 +1202,7 @@ pure subroutine scaling_1rescl(ncol, nlay, ngpt, tauLoc, scaling, tau, ssa, g) real(wp), dimension(ncol, nlay, ngpt), intent(inout) :: tauLoc real(wp), dimension(ncol, nlay, ngpt), intent(inout) :: scaling + integer :: icol, ilay, igpt real(wp) :: wb, ssal, scaleTau do igpt=1,ngpt @@ -1223,18 +1270,20 @@ end subroutine scaling_1rescl_safe ! ! ------------------------------------------------------------------------------------------------- subroutine lw_transport_1rescl(ncol, nlay, top_at_1, & - trans, sfc_albedo, source_dn, source_up, source_sfc, & - radn_up, radn_dn, An, Cn) bind(C, name="lw_transport_1rescl") + trans, source_dn, source_up, & + radn_up, radn_dn, An, Cn,& + radn_up_Jac, radn_dn_Jac) bind(C, name="lw_transport_1rescl") integer, intent(in ) :: ncol, nlay ! Number of columns, layers, g-points logical(wl), intent(in ) :: top_at_1 ! real(wp), dimension(ncol,nlay ), intent(in ) :: trans ! transmissivity = exp(-tau) - real(wp), dimension(ncol ), intent(in ) :: sfc_albedo ! Surface albedo real(wp), dimension(ncol,nlay ), intent(in ) :: source_dn, & source_up ! Diffuse radiation emitted by the layer - real(wp), dimension(ncol ), intent(in ) :: source_sfc ! Surface source function [W/m2] real(wp), dimension(ncol,nlay+1), intent(inout) :: radn_up ! Radiances [W/m2-str] real(wp), dimension(ncol,nlay+1), intent(inout) :: radn_dn !Top level must contain incident flux boundary condition real(wp), dimension(ncol,nlay), intent(in ) :: An, Cn + real(wp), dimension(ncol,nlay+1), intent(inout) :: radn_up_Jac ! Surface temperature Jacobians [W/m2-str/K] + real(wp), dimension(ncol,nlay+1), intent(inout) :: radn_dn_Jac !Top level must set to 0 + ! Local variables integer :: ilev, icol ! --------------------------------------------------- @@ -1245,21 +1294,26 @@ subroutine lw_transport_1rescl(ncol, nlay, top_at_1, & ! ! 1st Upward propagation do ilev = nlay, 1, -1 - radn_up(:,ilev) = trans(:,ilev )*radn_up(:,ilev+1) + source_up(:,ilev) + radn_up (:,ilev) = trans(:,ilev)*radn_up (:,ilev+1) + source_up(:,ilev) + radn_up_Jac(:,ilev) = trans(:,ilev)*radn_up_Jac(:,ilev+1) do icol=1,ncol adjustmentFactor = Cn(icol,ilev)*( An(icol,ilev)*radn_dn(icol,ilev) - & trans(icol,ilev)*source_dn(icol,ilev) - source_up(icol,ilev) ) - radn_up(icol,ilev) = radn_up(icol,ilev) + adjustmentFactor + radn_up (icol,ilev) = radn_up(icol,ilev) + adjustmentFactor enddo end do ! 2nd Downward propagation do ilev = 1, nlay - radn_dn(:,ilev+1) = trans(:,ilev)*radn_dn(:,ilev) + source_dn(:,ilev) + radn_dn (:,ilev+1) = trans(:,ilev)*radn_dn (:,ilev) + source_dn(:,ilev) + radn_dn_Jac(:,ilev+1) = trans(:,ilev)*radn_dn_Jac(:,ilev) do icol=1,ncol adjustmentFactor = Cn(icol,ilev)*( An(icol,ilev)*radn_up(icol,ilev) - & trans(icol,ilev)*source_up(icol,ilev) - source_dn(icol,ilev) ) - radn_dn(:,ilev+1) = radn_dn(:,ilev+1) + adjustmentFactor + radn_dn (icol,ilev+1) = radn_dn(icol,ilev+1) + adjustmentFactor + + adjustmentFactor = Cn(icol,ilev)*An(icol,ilev)*radn_up_Jac(icol,ilev) + radn_dn_Jac(icol,ilev+1) = radn_dn_Jac(icol,ilev+1) + adjustmentFactor enddo end do else @@ -1268,7 +1322,8 @@ subroutine lw_transport_1rescl(ncol, nlay, top_at_1, & ! ! Upward propagation do ilev = 1, nlay - radn_up(:,ilev+1) = trans(:,ilev) * radn_up(:,ilev) + source_up(:,ilev) + radn_up (:,ilev+1) = trans(:,ilev) * radn_up (:,ilev) + source_up(:,ilev) + radn_up_Jac(:,ilev+1) = trans(:,ilev) * radn_up_Jac(:,ilev) do icol=1,ncol adjustmentFactor = Cn(icol,ilev)*( An(icol,ilev)*radn_dn(icol,ilev+1) - & trans(icol,ilev)*source_dn(icol,ilev) - source_up(icol,ilev) ) @@ -1278,11 +1333,15 @@ subroutine lw_transport_1rescl(ncol, nlay, top_at_1, & ! 2st Downward propagation do ilev = nlay, 1, -1 - radn_dn(:,ilev) = trans(:,ilev )*radn_dn(:,ilev+1) + source_dn(:,ilev) + radn_dn (:,ilev) = trans(:,ilev)*radn_dn (:,ilev+1) + source_dn(:,ilev) + radn_dn_Jac(:,ilev) = trans(:,ilev)*radn_dn_Jac(:,ilev+1) do icol=1,ncol adjustmentFactor = Cn(icol,ilev)*( An(icol,ilev)*radn_up(icol,ilev) - & trans(icol,ilev)*source_up(icol,ilev) - source_dn(icol,ilev) ) - radn_dn(icol,ilev) = radn_dn(icol,ilev) + adjustmentFactor + radn_dn(icol,ilev) = radn_dn(icol,ilev) + adjustmentFactor + + adjustmentFactor = Cn(icol,ilev)*An(icol,ilev)*radn_up_Jac(icol,ilev) + radn_dn_Jac(icol,ilev) = radn_dn_Jac(icol,ilev) + adjustmentFactor enddo end do end if diff --git a/rte/mo_optical_props.F90 b/rte/mo_optical_props.F90 index b18751f3e..ffc8a3d7a 100644 --- a/rte/mo_optical_props.F90 +++ b/rte/mo_optical_props.F90 @@ -809,13 +809,15 @@ function increment(op_in, op_io) result(err_message) integer :: ncol, nlay, ngpt, nmom ! ----- err_message = "" - if(.not. op_in%bands_are_equal(op_io)) then - err_message = "ty_optical_props%increment: optical properties objects have different band structures" - return - end if ncol = op_io%get_ncol() nlay = op_io%get_nlay() ngpt = op_io%get_ngpt() + if(.not. op_in%bands_are_equal(op_io)) & + err_message = "ty_optical_props%increment: optical properties objects have different band structures" + if(.not. all([op_in%get_ncol(), op_in%get_nlay()] == [ncol, nlay])) & + err_message = "ty_optical_props%increment: optical properties objects have different ncol and/or nlay" + if(err_message /= "") return + if(op_in%gpoints_are_equal(op_io)) then ! ! Increment by gpoint diff --git a/rte/mo_rte_lw.F90 b/rte/mo_rte_lw.F90 index 25dde3cad..7ff091b2a 100644 --- a/rte/mo_rte_lw.F90 +++ b/rte/mo_rte_lw.F90 @@ -40,9 +40,9 @@ module mo_rte_lw ty_optical_props_arry, ty_optical_props_1scl, ty_optical_props_2str, ty_optical_props_nstr use mo_source_functions, & only: ty_source_func_lw - use mo_fluxes, only: ty_fluxes + use mo_fluxes, only: ty_fluxes, ty_fluxes_broadband use mo_rte_solver_kernels, & - only: apply_BC, lw_solver_noscat_GaussQuad, lw_solver_2stream,& + only: apply_BC, lw_solver_noscat, lw_solver_noscat_GaussQuad, lw_solver_2stream, & lw_solver_1rescl_GaussQuad implicit none private @@ -57,33 +57,43 @@ module mo_rte_lw function rte_lw(optical_props, top_at_1, & sources, sfc_emis, & fluxes, & - inc_flux, n_gauss_angles, use_2stream) result(error_msg) + inc_flux, n_gauss_angles, use_2stream, & + lw_Ds, flux_up_Jac, flux_dn_Jac) result(error_msg) class(ty_optical_props_arry), intent(in ) :: optical_props ! Array of ty_optical_props. This type is abstract ! and needs to be made concrete, either as an array ! (class ty_optical_props_arry) or in some user-defined way logical, intent(in ) :: top_at_1 ! Is the top of the domain at index 1? ! (if not, ordering is bottom-to-top) type(ty_source_func_lw), intent(in ) :: sources - real(wp), dimension(:,:), intent(in ) :: sfc_emis ! emissivity at surface [] (nband, ncol) - class(ty_fluxes), intent(inout) :: fluxes ! Array of ty_fluxes. Default computes broadband fluxes at all levels - ! if output arrays are defined. Can be extended per user desires. + real(wp), dimension(:,:), intent(in ) :: sfc_emis ! emissivity at surface [] (nband, ncol) + class(ty_fluxes), intent(inout) :: fluxes ! Array of ty_fluxes. Default computes broadband fluxes at all levels + ! if output arrays are defined. Can be extended per user desires. real(wp), dimension(:,:), & target, optional, intent(in ) :: inc_flux ! incident flux at domain top [W/m2] (ncol, ngpts) integer, optional, intent(in ) :: n_gauss_angles ! Number of angles used in Gaussian quadrature ! (no-scattering solution) logical, optional, intent(in ) :: use_2stream ! When 2-stream parameters (tau/ssa/g) are provided, use 2-stream methods ! Default is to use re-scaled longwave transport - character(len=128) :: error_msg ! If empty, calculation was successful + real(wp), dimension(:,:), & + optional, intent(in ) :: lw_Ds ! linear fit to column transmissivity (ncol,ngpt) + real(wp), dimension(:,:), & + target, optional, intent(inout) :: flux_up_Jac ! surface temperature flux Jacobian [W/m2/K] (ncol, ngpts) + real(wp), dimension(:,:), & + target, optional, intent(inout) :: flux_dn_Jac ! surface temperature flux Jacobian [W/m2/K] (ncol, ngpts) + character(len=128) :: error_msg ! If empty, calculation was successful ! -------------------------------- ! ! Local variables ! - integer :: ncol, nlay, ngpt, nband - integer :: n_quad_angs - integer :: icol, iband, igpt + integer :: ncol, nlay, ngpt, nband + integer :: n_quad_angs + integer :: icol, iband, igpt + real(wp) :: lw_Ds_wt + logical :: using_2stream real(wp), dimension(:,:,:), allocatable :: gpt_flux_up, gpt_flux_dn real(wp), dimension(:,:), allocatable :: sfc_emis_gpt - logical :: using_2stream + real(wp), dimension(:,:,:), allocatable :: gpt_flux_upJac, gpt_flux_dnJac + type(ty_fluxes_broadband) :: Jac_fluxes ! -------------------------------------------------- ! ! Weights and angle secants for first order (k=1) Gaussian quadrature. @@ -117,13 +127,16 @@ function rte_lw(optical_props, top_at_1, & ! ------------------------------------------------------------------------------------ ! - ! Error checking -- consistency of sizes and validity of values + ! Error checking -- input consistency of sizes and validity of values ! ! -------------------------------- - if(.not. fluxes%are_desired()) then + + if(.not. fluxes%are_desired()) & error_msg = "rte_lw: no space allocated for fluxes" - return - end if + if (present(flux_up_Jac)) then + if( .not. extents_are(flux_up_Jac, ncol, nlay+1)) & + error_msg = "rte_lw: flux Jacobian inconsistently sized" + endif ! ! Source functions @@ -158,18 +171,48 @@ function rte_lw(optical_props, top_at_1, & n_quad_angs = 1 if(present(n_gauss_angles)) then if(n_gauss_angles > max_gauss_pts) & - error_msg = "rte_lw: asking for too many quadrature points for no-scattering calculation" + error_msg = "rte_lw: asking for too many quadrature points for RT calculation" if(n_gauss_angles < 1) & - error_msg = "rte_lw: have to ask for at least one quadrature point for no-scattering calculation" + error_msg = "rte_lw: have to ask for at least one quadrature point for RT calculation" n_quad_angs = n_gauss_angles end if + if(len_trim(error_msg) > 0) return + ! ! Optionally - use 2-stream methods when low-order scattering properties are provided? ! using_2stream = .false. if(present(use_2stream)) using_2stream = use_2stream + + ! + ! Checking that optional arguements are consistent with one another and with optical properties + ! + select type (optical_props) + class is (ty_optical_props_1scl) + if (using_2stream) & + error_msg = "rte_lw: can't use two-stream methods with only absorption optical depth" + if (present(lw_Ds)) then + if(.not. extents_are(lw_Ds, ncol, ngpt)) & + error_msg = "rte_lw: lw_Ds inconsistently sized" + if(any_vals_less_than(lw_Ds, 1._wp)) & + error_msg = "rte_lw: one or more values of lw_Ds < 1." + if(n_quad_angs /= 1) & + error_msg = "rte_lw: providing lw_Ds incompatible with specifying n_gauss_angles" + end if + class is (ty_optical_props_2str) + if (present(lw_Ds)) & + error_msg = "rte_lw: lw_Ds not valid input for _2str class" + if (using_2stream .and. n_quad_angs /= 1) & + error_msg = "rte_lw: using_2stream=true incompatible with specifying n_gauss_angles" + if (using_2stream .and. (present(flux_up_Jac) .or. present(flux_up_Jac))) & + error_msg = "rte_lw: can't provide Jacobian of fluxes w.r.t surface temperature with 2-stream" + class default + call stop_on_err("rte_lw: lw_solver(...ty_optical_props_nstr...) not yet implemented") + end select + if(len_trim(error_msg) > 0) return + ! - ! Ensure values of tau, ssa, and g are reasonable + ! Ensure values of tau, ssa, and g are reasonable if using scattering ! error_msg = optical_props%validate() if(len_trim(error_msg) > 0) then @@ -183,11 +226,14 @@ function rte_lw(optical_props, top_at_1, & ! Lower boundary condition -- expand surface emissivity by band to gpoints ! allocate(gpt_flux_up (ncol, nlay+1, ngpt), gpt_flux_dn(ncol, nlay+1, ngpt)) + allocate(gpt_flux_upJac(ncol, nlay+1, ngpt)) allocate(sfc_emis_gpt(ncol, ngpt)) !!$acc enter data copyin(sources, sources%lay_source, sources%lev_source_inc, sources%lev_source_dec, sources%sfc_source) !$acc enter data copyin(optical_props) !$acc enter data create(gpt_flux_dn, gpt_flux_up) + !$acc enter data create(gpt_flux_upJac) !$acc enter data create(sfc_emis_gpt) + call expand_and_transpose(optical_props, sfc_emis, sfc_emis_gpt) ! ! Upper boundary condition @@ -203,6 +249,7 @@ function rte_lw(optical_props, top_at_1, & call apply_BC(ncol, nlay, ngpt, logical(top_at_1, wl), gpt_flux_dn) end if + ! ! Compute the radiative transfer... ! @@ -214,14 +261,35 @@ function rte_lw(optical_props, top_at_1, & !$acc enter data copyin(optical_props%tau) error_msg = optical_props%validate() if(len_trim(error_msg) > 0) return - call lw_solver_noscat_GaussQuad(ncol, nlay, ngpt, logical(top_at_1, wl), & - n_quad_angs, gauss_Ds(1:n_quad_angs,n_quad_angs), gauss_wts(1:n_quad_angs,n_quad_angs), & - optical_props%tau, & - sources%lay_source, sources%lev_source_inc, sources%lev_source_dec, & - sfc_emis_gpt, sources%sfc_source, & - gpt_flux_up, gpt_flux_dn) + + if (present(lw_Ds)) then + call lw_solver_noscat(ncol, nlay, ngpt, & + logical(top_at_1, wl), & + lw_Ds, gauss_wts(1,1), & + optical_props%tau, & + sources%lay_source, sources%lev_source_inc, sources%lev_source_dec, & + sfc_emis_gpt, sources%sfc_source, & + gpt_flux_up, gpt_flux_dn, sources%sfc_source_Jac, gpt_flux_upJac) + else + call lw_solver_noscat_GaussQuad(ncol, nlay, ngpt, & + logical(top_at_1, wl), & + n_quad_angs, & + gauss_Ds(1:n_quad_angs,n_quad_angs), & + gauss_wts(1:n_quad_angs,n_quad_angs), & + optical_props%tau, & + sources%lay_source, sources%lev_source_inc, & + sources%lev_source_dec, & + sfc_emis_gpt, sources%sfc_source, & + gpt_flux_up, gpt_flux_dn, sources%sfc_source_Jac, gpt_flux_upJac) + end if !$acc exit data delete(optical_props%tau) class is (ty_optical_props_2str) + if ((present(flux_dn_Jac))) then + if( .not. extents_are(flux_dn_Jac, ncol, nlay+1)) then + error_msg = "rte_lw: flux_dn_Jac inconsistently sized" + return + end if + endif if (using_2stream) then ! ! two-stream calculation with scattering @@ -239,6 +307,8 @@ function rte_lw(optical_props, top_at_1, & ! ! Re-scaled solution to account for scattering ! + allocate(gpt_flux_dnJac (ncol, nlay+1, ngpt)) + !$acc enter data create(gpt_flux_dnJac) !$acc enter data copyin(optical_props%tau, optical_props%ssa, optical_props%g) call lw_solver_1rescl_GaussQuad(ncol, nlay, ngpt, logical(top_at_1, wl), & n_quad_angs, gauss_Ds(1:n_quad_angs,n_quad_angs), & @@ -247,7 +317,8 @@ function rte_lw(optical_props, top_at_1, & sources%lay_source, sources%lev_source_inc, & sources%lev_source_dec, & sfc_emis_gpt, sources%sfc_source,& - gpt_flux_up, gpt_flux_dn) + gpt_flux_up, gpt_flux_dn, & + sources%sfc_source_Jac, gpt_flux_upJac, gpt_flux_dnJac) !$acc exit data delete(optical_props%tau, optical_props%ssa, optical_props%g) endif class is (ty_optical_props_nstr) @@ -262,6 +333,30 @@ function rte_lw(optical_props, top_at_1, & ! ...and reduce spectral fluxes to desired output quantities ! error_msg = fluxes%reduce(gpt_flux_up, gpt_flux_dn, optical_props, top_at_1) + if (error_msg /= '') return + + if (present(flux_up_Jac)) Jac_fluxes%flux_up => flux_up_Jac + select type (optical_props) + class is (ty_optical_props_1scl) + ! + ! gpoint Jacobian fluxes aren't defined for _1scl + ! + error_msg = Jac_fluxes%reduce(gpt_flux_upJac, gpt_flux_upJac, optical_props, top_at_1) + class is (ty_optical_props_2str) + ! + ! Compute Jacobians when using rescaling approach for scattering + ! + if(.not. using_2stream) then + if (present(flux_dn_Jac)) Jac_fluxes%flux_dn => flux_dn_Jac + error_msg = Jac_fluxes%reduce(gpt_flux_upJac, gpt_flux_dnJac, optical_props, top_at_1) + !$acc exit data delete(gpt_flux_dnJac) + deallocate(gpt_flux_dnJac) + end if + end select + + !$acc exit data delete(gpt_flux_upJac) + deallocate(gpt_flux_upJac) + !$acc exit data delete(sfc_emis_gpt) !$acc exit data delete(gpt_flux_up,gpt_flux_dn) !$acc exit data delete(optical_props) diff --git a/rte/mo_source_functions.F90 b/rte/mo_source_functions.F90 index af4dc4c79..859303386 100644 --- a/rte/mo_source_functions.F90 +++ b/rte/mo_source_functions.F90 @@ -32,15 +32,16 @@ module mo_source_functions ! Includes spectral weighting that accounts for state-dependent ! frequency to g-space mapping real(wp), allocatable, dimension(:,: ) :: sfc_source + real(wp), allocatable, dimension(:,: ) :: sfc_source_Jac ! surface source Jacobian contains generic, public :: alloc => alloc_lw, copy_and_alloc_lw procedure, private:: alloc_lw procedure, private:: copy_and_alloc_lw procedure, public :: is_allocated => is_allocated_lw - procedure, public :: finalize => finalize_lw - procedure, public :: get_subset => get_subset_range_lw - procedure, public :: get_ncol => get_ncol_lw - procedure, public :: get_nlay => get_nlay_lw + procedure, public :: finalize => finalize_lw + procedure, public :: get_subset => get_subset_range_lw + procedure, public :: get_ncol => get_ncol_lw + procedure, public :: get_nlay => get_nlay_lw ! validate? end type ty_source_func_lw ! ------------------------------------------------------------------------------------------------- @@ -92,14 +93,16 @@ function alloc_lw(this, ncol, nlay) result(err_message) err_message = "source_func_lw%alloc: must provide positive extents for ncol, nlay" if (err_message /= "") return - if(allocated(this%sfc_source)) deallocate(this%sfc_source) - if(allocated(this%lay_source)) deallocate(this%lay_source) + if(allocated(this%sfc_source)) deallocate(this%sfc_source) + if(allocated(this%sfc_source_Jac)) deallocate(this%sfc_source_Jac) + if(allocated(this%lay_source)) deallocate(this%lay_source) if(allocated(this%lev_source_inc)) deallocate(this%lev_source_inc) if(allocated(this%lev_source_dec)) deallocate(this%lev_source_dec) ngpt = this%get_ngpt() allocate(this%sfc_source (ncol, ngpt), this%lay_source (ncol,nlay,ngpt), & this%lev_source_inc(ncol,nlay,ngpt), this%lev_source_dec(ncol,nlay,ngpt)) + allocate(this%sfc_source_Jac(ncol, ngpt)) end function alloc_lw ! -------------------------------------------------------------- function copy_and_alloc_lw(this, ncol, nlay, spectral_desc) result(err_message) @@ -175,6 +178,7 @@ subroutine finalize_lw(this) if(allocated(this%lev_source_inc)) deallocate(this%lev_source_inc) if(allocated(this%lev_source_dec)) deallocate(this%lev_source_dec) if(allocated(this%sfc_source )) deallocate(this%sfc_source) + if(allocated(this%sfc_source_Jac)) deallocate(this%sfc_source_Jac) call this%ty_optical_props%finalize() end subroutine finalize_lw ! -------------------------------------------------------------- @@ -248,6 +252,7 @@ function get_subset_range_lw(full, start, n, subset) result(err_message) err_message = subset%alloc(n, full%get_nlay(), full) if(err_message /= "") return subset%sfc_source (1:n, :) = full%sfc_source (start:start+n-1, :) + subset%sfc_source_Jac(1:n, :) = full%sfc_source_Jac(start:start+n-1, :) subset%lay_source (1:n,:,:) = full%lay_source (start:start+n-1,:,:) subset%lev_source_inc(1:n,:,:) = full%lev_source_inc(start:start+n-1,:,:) subset%lev_source_dec(1:n,:,:) = full%lev_source_dec(start:start+n-1,:,:) diff --git a/tests/Makefile b/tests/Makefile new file mode 100644 index 000000000..dbc0906df --- /dev/null +++ b/tests/Makefile @@ -0,0 +1,66 @@ +# Location of RTE+RRTMGP libraries, module files. +RRTMGP_BUILD = $(RRTMGP_ROOT)/build +# Sets macros FC, FCFLAGS consistent with RTE+RRTMGP +-include $(RRTMGP_BUILD)/Makefile.conf + +# Location of netcdf C and Fortran libraries. Could specify with environment variables if file doesn't exist +-include $(RRTMGP_ROOT)/examples/rfmip-clear-sky/Makefile.libs +# +# RRTMGP library, module files +# +LDFLAGS += -L$(RRTMGP_BUILD) +LIBS += -lrrtmgp -lrte +FCINCLUDE += -I$(RRTMGP_BUILD) + +# +# netcdf library, module files +# C and Fortran interfaces respectively +# +FCINCLUDE += -I$(NFHOME)/include +LDFLAGS += -L$(NFHOME)/lib -L$(NCHOME)/lib +LIBS += -lnetcdff -lnetcdf + +VPATH = .:$(RRTMGP_ROOT)/examples:$(RRTMGP_ROOT)/examples/rfmip-clear-sky +VPATH += $(RRTMGP_ROOT)/extensions:$(RRTMGP_ROOT)/extensions/cloud_optics:$(RRTMGP_ROOT)/extensions/solar_variability + +# Compilation rules +%.o: %.F90 + $(FC) $(FCFLAGS) $(FCINCLUDE) -c $< +%: %.o + $(FC) $(FCFLAGS) -o $@ $^ $(LDFLAGS) $(LIBS) + + +# +# Extra sources -- extensions to RRTMGP classes, shared infrastructure, local sources +# +ADDITIONS = mo_heating_rates.o mo_compute_bc.o mo_rrtmgp_clr_all_sky.o +# File I/O +ADDITIONS += mo_load_coefficients.o mo_simple_netcdf.o mo_rfmip_io.o +ADDITIONS += mo_testing_io.o +# Cloud optics +CLOUDS += mo_cloud_sampling.o mo_cloud_optics.o mo_load_cloud_coefficients.o +# Solar variability +ADDITIONS += mo_solar_variability.o + +# Many codes will need to be updated if the library changes +LIB_DEPS = $(RRTMGP_BUILD)/librte.a $(RRTMGP_BUILD)/librrtmgp.a +# +# Targets +# +all: clear_sky_regression + +clear_sky_regression: $(ADDITIONS) $(LIB_DEPS) clear_sky_regression.o +clear_sky_regression.o: $(ADDITIONS) $(LIB_DEPS) clear_sky_regression.F90 + +mo_testing_io.o: $(LIB_DEPS) mo_simple_netcdf.o mo_testing_io.F90 + +mo_cloud_optics.o: $(LIB_DEPS) mo_cloud_optics.F90 +mo_load_cloud_coefficients.o: $(LIB_DEPS) mo_simple_netcdf.o mo_cloud_optics.o mo_load_cloud_coefficients.F90 +mo_cloud_sampling.o: $(LIB_DEPS) mo_cloud_sampling.F90 + +mo_load_coefficients.o: $(LIB_DEPS) mo_simple_netcdf.o mo_load_coefficients.F90 +mo_rfmip_io.o.o: $(LIB_DEPS) mo_simple_netcdf.o mo_rfmip_io.F90 +mo_simple_netcdf.o: $(LIB_DEPS) mo_simple_netcdf.F90 + +clean: + -rm clear_sky_regression *.o *.optrpt diff --git a/tests/clear_sky_regression.F90 b/tests/clear_sky_regression.F90 new file mode 100644 index 000000000..2d070cd5f --- /dev/null +++ b/tests/clear_sky_regression.F90 @@ -0,0 +1,748 @@ +subroutine stop_on_err(error_msg) + use iso_fortran_env, only : error_unit + character(len=*), intent(in) :: error_msg + + if(error_msg /= "") then + write (error_unit,*) trim(error_msg) + write (error_unit,*) "clear-sky regression tests stopping" + error stop 1 + end if +end subroutine stop_on_err +! ---------------------------------------------------------------------------------- +program rte_clear_sky_regression + use mo_rte_kind, only: wp + use mo_optical_props, only: ty_optical_props, & + ty_optical_props_arry, & + ty_optical_props_1scl, ty_optical_props_2str, ty_optical_props_nstr + use mo_gas_optics_rrtmgp, only: ty_gas_optics_rrtmgp + use mo_gas_concentrations, only: ty_gas_concs + use mo_source_functions, only: ty_source_func_lw + use mo_fluxes, only: ty_fluxes_broadband + use mo_rte_lw, only: rte_lw + use mo_rte_sw, only: rte_sw + use mo_load_coefficients, only: load_and_init + use mo_rfmip_io, only: read_size, read_and_block_pt, read_and_block_gases_ty, & + read_and_block_lw_bc, read_and_block_sw_bc, determine_gas_names + use mo_simple_netcdf, only: get_dim_size, read_field + use mo_heating_rates, only: compute_heating_rate + use mo_testing_io, only: write_broadband_field + implicit none + ! ---------------------------------------------------------------------------------- + ! Variables + ! ---------------------------------------------------------------------------------- + ! Arrays: dimensions (col, lay, exp), as read from RFMIP test cases + real(wp), dimension(:,:,:), allocatable :: p_lay_3d, t_lay_3d, p_lev_3d + real(wp), dimension(:,:,:), allocatable :: t_lev_3d + real(wp), dimension(:,:), allocatable :: sfc_t_3d + real(wp), dimension(:,:), allocatable :: bc_3d, tsi_3d ! boundary conditions (ncol, nexp) + real(wp), dimension(:,:), allocatable :: sza ! Solar zenith angle (ncol, nexp) + + ! + ! Local versions of variables - used + ! + real(wp), dimension(:,:), allocatable :: p_lay, t_lay, p_lev + ! + ! Longwave only + ! + real(wp), dimension(:,:), allocatable :: t_lev + real(wp), dimension(:), allocatable :: sfc_t + real(wp), dimension(:,:), allocatable :: sfc_emis ! First dimension is band + ! + ! Shortwave only + ! + real(wp), dimension(:), allocatable :: mu0 ! solar zenith angle, cosine of same + real(wp), dimension(:,:), allocatable :: sfc_alb_dir, sfc_alb_dif ! First dimension is band + logical, dimension(:), allocatable :: sun_up + ! + ! Source functions + ! + ! Longwave + type(ty_source_func_lw) :: lw_sources + ! Shortwave + real(wp), dimension(:,:), allocatable :: toa_flux + ! + ! Output variables + ! + real(wp), dimension(:,:), target, & + allocatable :: flux_up, flux_dn, flux_dir + ! + ! Derived types from the RTE and RRTMGP libraries + ! + type(ty_gas_optics_rrtmgp) :: k_dist + type(ty_gas_concs) :: gas_concs + type(ty_gas_concs), dimension(:), allocatable & + :: gas_conc_array + class(ty_optical_props_arry), & + allocatable :: atmos + type(ty_fluxes_broadband) :: fluxes + + ! + ! Inputs to RRTMGP + ! + logical :: top_at_1, is_sw, is_lw + + integer :: ncol, nlay, nbnd, ngpt, nexp + integer :: icol, ilay, ibnd, iloop, igas + + integer :: nUserArgs=0 + + character(len=32 ), & + dimension(:), allocatable :: kdist_gas_names, rfmip_gas_games + + character(len=256) :: input_file, k_dist_file + ! ---------------------------------------------------------------------------------- + ! Code + ! ---------------------------------------------------------------------------------- + ! + ! Parse command line for any file names, block size + ! + nUserArgs = command_argument_count() + if (nUserArgs < 2) call stop_on_err("Need to supply input_file k_distribution_file") + if (nUserArgs >= 1) call get_command_argument(1,input_file) + if (nUserArgs >= 2) call get_command_argument(2,k_dist_file) + if (nUserArgs > 3) print *, "Ignoring command line arguments beyond the first three..." + if(trim(input_file) == '-h' .or. trim(input_file) == "--help") then + call stop_on_err("clear_sky_regression input_file absorption_coefficients_file") + end if + ! + ! Read temperature, pressure, gas concentrations. + ! Arrays are allocated as they are read + ! + call read_size(input_file, ncol, nlay, nexp) + call determine_gas_names(input_file, k_dist_file, 1, kdist_gas_names, rfmip_gas_games) + call read_and_block_pt(input_file, ncol, p_lay_3d, p_lev_3d, t_lay_3d, t_lev_3d) + ! + ! Only do the first RFMIP experiment + ! + allocate(p_lay(ncol, nlay), p_lev(ncol, nlay+1), & + t_lay(ncol, nlay), t_lev(ncol, nlay+1)) + p_lay(:,:) = p_lay_3d(:,:,1) + p_lev(:,:) = p_lev_3d(:,:,1) + t_lay(:,:) = t_lay_3d(:,:,1) + t_lev(:,:) = t_lev_3d(:,:,1) + deallocate(p_lay_3d, p_lev_3d, t_lay_3d, t_lev_3d) + ! + ! Read the gas concentrations and surface properties + ! RFMIP I/O returns an array, we're going to use first ncol values = experiement 1 (present-day) + ! + call read_and_block_gases_ty(input_file, ncol*nexp, kdist_gas_names, rfmip_gas_games, gas_conc_array) + ! + ! All the profiles are in the first and only element of the array of gas concentration types + ! Extract the first ncol profiles (this is part of testing) + ! + call stop_on_err(gas_conc_array(1)%get_subset(1, ncol, gas_concs)) + call gas_conc_array(1)%reset() + deallocate(gas_conc_array) + ! ---------------------------------------------------------------------------- + ! load data into classes + call load_and_init(k_dist, k_dist_file, gas_concs) + is_sw = k_dist%source_is_external() + 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() + ! + ! Problem sizes + ! + nbnd = k_dist%get_nband() + ngpt = k_dist%get_ngpt() + ! + ! RRTMGP won't run with pressure less than its minimum. The top level in the RFMIP file + ! is set to 10^-3 Pa. Here we pretend the layer is just a bit less deep. + ! + top_at_1 = p_lay(1, 1) < p_lay(1, nlay) + if(top_at_1) then + p_lev(:,1) = k_dist%get_press_min() + epsilon(k_dist%get_press_min()) + else + p_lev(:,nlay+1) & + = k_dist%get_press_min() + epsilon(k_dist%get_press_min()) + end if + ! ---------------------------------------------------------------------------- + ! + ! Boundary conditions + ! + if(is_sw) then + allocate(toa_flux(ncol, ngpt)) + allocate(sfc_alb_dir(nbnd, ncol), sfc_alb_dif(nbnd, ncol), & + mu0(ncol), sun_up(ncol)) + call read_and_block_sw_bc(input_file, ncol, & + bc_3d, tsi_3d, sza) + ! + ! Surface albedo is spectrally uniform + ! + sfc_alb_dir(:,:) = spread(bc_3d(:,1), dim=1, ncopies=nbnd) + sfc_alb_dif(:,:) = sfc_alb_dir(:,:) + sun_up(:) = sza(:,1) < 87.5_wp ! Limits is from LBLRTM + mu0(:) = merge(cos(sza(:,1) * acos(-1._wp)/180._wp), & + 1._wp, & + sun_up) + else + allocate(sfc_t(ncol), sfc_emis(nbnd, ncol)) + call stop_on_err(lw_sources%alloc(ncol, nlay, k_dist)) + call read_and_block_lw_bc(input_file, ncol, bc_3d, sfc_t_3d) + ! + ! Surface emissivity is spectrally uniform + ! + sfc_emis(:,:) = spread(bc_3d(:,1), dim=1, ncopies=nbnd) + + sfc_t (:) = sfc_t_3d(:,1) + end if + ! ---------------------------------------------------------------------------- + ! + ! Fluxes + ! + allocate(flux_up(ncol,nlay+1), flux_dn(ncol,nlay+1), flux_dir(ncol,nlay+1)) + ! ---------------------------------------------------------------------------- + ! + ! Solvers + ! + fluxes%flux_up => flux_up(:,:) + fluxes%flux_dn => flux_dn(:,:) + if(is_lw) then + call make_optical_props_1scl + call atmos%set_name("gas only atmosphere") + call lw_clear_sky_default + call lw_clear_sky_notlev + call lw_clear_sky_3ang + call lw_clear_sky_optangle + call lw_clear_sky_jaco + call lw_clear_sky_subset + call lw_clear_sky_vr + call lw_clear_sky_incr + call make_optical_props_2str + call lw_clear_sky_2str + else + call make_optical_props_2str + call sw_clear_sky_default + call sw_clear_sky_tsi + call sw_clear_sky_vr + end if +contains + ! ---------------------------------------------------------------------------- + ! + ! Clear-sky longwave fluxes, + ! + subroutine lw_clear_sky_default + real(wp), dimension(ncol, nlay+1), target :: flux_net + real(wp), dimension(ncol, nlay) :: heating_rate + + fluxes%flux_net => flux_net + call stop_on_err(k_dist%gas_optics(p_lay, p_lev, & + t_lay, sfc_t, & + gas_concs, & + atmos, & + lw_sources, & + tlev = t_lev)) + call stop_on_err(rte_lw(atmos, top_at_1, & + lw_sources, & + sfc_emis, & + fluxes)) + call write_broadband_field(input_file, flux_up, "lw_flux_up", "LW flux up") + call write_broadband_field(input_file, flux_dn, "lw_flux_dn", "LW flux dn") + call write_broadband_field(input_file, flux_net, "lw_flux_net", "LW flux net") + + 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_default", "LW heating rate", vert_dim_name = "layer") + ! + ! Test for mo_fluxes_broadband for computing only net flux + ! + nullify(fluxes%flux_up) + nullify(fluxes%flux_dn) + call stop_on_err(rte_lw(atmos, top_at_1, & + lw_sources, & + sfc_emis, & + fluxes)) + call write_broadband_field(input_file, flux_net, "lw_flux_net_2", "LW flux net, direct") + fluxes%flux_up => flux_up + fluxes%flux_dn => flux_dn + nullify(fluxes%flux_net) + + end subroutine lw_clear_sky_default + ! ---------------------------------------------------------------------------- + ! + ! Clear-sky longwave fluxes, level temperatures provided + ! + subroutine lw_clear_sky_notlev + call stop_on_err(k_dist%gas_optics(p_lay, p_lev, & + t_lay, sfc_t, & + gas_concs, & + atmos, & + lw_sources)) + call stop_on_err(rte_lw(atmos, top_at_1, & + lw_sources, & + sfc_emis, & + fluxes)) + call write_broadband_field(input_file, flux_up, "lw_flux_up_notlev", "LW flux up, no level temperatures") + call write_broadband_field(input_file, flux_dn, "lw_flux_dn_notlev", "LW flux dn, no level temperatures") + end subroutine lw_clear_sky_notlev + ! ---------------------------------------------------------------------------- + ! + ! Clear-sky longwave fluxes, all info, three angles + ! + subroutine lw_clear_sky_3ang + call stop_on_err(k_dist%gas_optics(p_lay, p_lev, & + t_lay, sfc_t, & + gas_concs, & + atmos, & + lw_sources, & + tlev = t_lev)) + call stop_on_err(rte_lw(atmos, top_at_1, & + lw_sources, & + sfc_emis, & + fluxes, n_gauss_angles=3)) + call write_broadband_field(input_file, flux_up, "lw_flux_up_3ang", "LW flux up, three quadrature angles") + call write_broadband_field(input_file, flux_dn, "lw_flux_dn_3ang", "LW flux dn, three quadrature angles") + end subroutine lw_clear_sky_3ang + ! ---------------------------------------------------------------------------- + ! + ! Clear-sky longwave fluxes, all info, three angles + ! + subroutine lw_clear_sky_optangle + real(wp), dimension(ncol, ngpt) :: lw_Ds + call stop_on_err(k_dist%gas_optics(p_lay, p_lev, & + t_lay, sfc_t, & + gas_concs, & + atmos, & + lw_sources, & + tlev = t_lev)) + call stop_on_err(k_dist%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_optang", "LW flux up, single optimal angles") + call write_broadband_field(input_file, flux_dn, "lw_flux_dn_optang", "LW flux dn, single optimal angles") + end subroutine lw_clear_sky_optangle + ! ---------------------------------------------------------------------------- + ! + ! Clear-sky longwave fluxes, all info, computing surface Jacobian + ! + subroutine lw_clear_sky_jaco + real(wp), dimension(ncol,nlay+1) :: jFluxUp + + call stop_on_err(k_dist%gas_optics(p_lay, p_lev, & + t_lay, sfc_t, & + gas_concs, & + atmos, & + lw_sources, & + tlev=t_lev)) + call stop_on_err(rte_lw(atmos, top_at_1, & + lw_sources, & + sfc_emis, & + fluxes, & + flux_up_Jac = jFluxUp)) + call write_broadband_field(input_file, flux_up, "lw_flux_up_jaco", "LW flux up, computing Jaobians") + call write_broadband_field(input_file, flux_dn, "lw_flux_dn_jaco", "LW flux dn, computing Jaobians") + call write_broadband_field(input_file, jFluxUp, "lw_jaco_up" , "Jacobian of LW flux up to surface temperature") + + call stop_on_err(k_dist%gas_optics(p_lay, p_lev, & + t_lay, sfc_t + 1._wp, & + gas_concs, & + atmos, & + lw_sources, & + tlev=t_lev)) + call stop_on_err(rte_lw(atmos, top_at_1, & + lw_sources, & + sfc_emis, & + fluxes)) + call write_broadband_field(input_file, flux_up, "lw_flux_up_stp1", "LW flux up, surface T+1K") + call write_broadband_field(input_file, flux_dn, "lw_flux_dn_stp1", "LW flux dn, surface T+1K") + end subroutine lw_clear_sky_jaco + ! ---------------------------------------------------------------------------- + ! + ! Clear-sky longwave fluxes, half the columns at a time + ! We're counting on ncol being even + ! + subroutine lw_clear_sky_subset + type(ty_optical_props_1scl) :: atmos_subset + type(ty_source_func_lw) :: sources_subset + type(ty_fluxes_broadband) :: fluxes ! Use local variable + real(wp), dimension(ncol/2, nlay+1), target & + :: up, dn + integer :: i, colS, colE + + call stop_on_err(k_dist%gas_optics(p_lay, p_lev, & + t_lay, sfc_t, & + gas_concs, & + atmos, & + lw_sources, & + tlev=t_lev)) + ! Testing init, then alloc for 1scl + call stop_on_err(atmos_subset%init(atmos)) + fluxes%flux_up => up + fluxes%flux_dn => dn + do i = 1, 2 + colS = ((i-1) * ncol/2) + 1 + colE = i * ncol/2 + call stop_on_err(atmos%get_subset (colS, ncol/2, atmos_subset)) + call stop_on_err(lw_sources%get_subset(colS, ncol/2, sources_subset)) + call stop_on_err(rte_lw(atmos_subset, top_at_1, & + sources_subset, & + sfc_emis(:, colS:colE), & + fluxes)) + flux_up(colS:colE,:) = up + flux_dn(colS:colE,:) = dn + end do + + call write_broadband_field(input_file, flux_up, "lw_flux_up_subset", "LW flux up, done in two parts") + call write_broadband_field(input_file, flux_dn, "lw_flux_dn_subset", "LW flux dn, done in two parts") + end subroutine lw_clear_sky_subset + ! ---------------------------------------------------------------------------- + ! + ! Clear-sky longwave fluxes, all info, three angles + ! Reverse orientation in the vertical, compute, un-reverse + ! + subroutine lw_clear_sky_vr + type(ty_gas_concs) :: gas_concs_vr + integer :: i + real(wp), dimension(ncol,nlay) :: vmr + character(32), & + dimension(gas_concs%get_num_gases()) & + :: gc_gas_names + ! + ! Reverse the orientation of the problem + ! + p_lay (:,:) = p_lay (:, nlay :1:-1) + t_lay (:,:) = t_lay (:, nlay :1:-1) + p_lev (:,:) = p_lev (:,(nlay+1):1:-1) + t_lev (:,:) = t_lev (:,(nlay+1):1:-1) + top_at_1 = .not. top_at_1 + ! + ! No direct access to gas concentrations so use the classes + ! This also tests otherwise uncovered routines for ty_gas_concs + ! + gc_gas_names(:) = gas_concs%get_gas_names() + call stop_on_err(gas_concs_vr%init(gc_gas_names(:))) + do i = 1, gas_concs%get_num_gases() + call stop_on_err(gas_concs%get_vmr(gc_gas_names(i), vmr)) + vmr(:,:) = vmr(:,nlay:1:-1) + call stop_on_err(gas_concs_vr%set_vmr(gc_gas_names(i), vmr)) + end do + + call stop_on_err(k_dist%gas_optics(p_lay, p_lev, & + t_lay, sfc_t, & + gas_concs_vr, & + atmos, & + lw_sources, & + tlev=t_lev)) + call stop_on_err(rte_lw(atmos, top_at_1, & + lw_sources, & + sfc_emis, & + fluxes)) + flux_up(:,:) = flux_up(:,(nlay+1):1:-1) + flux_dn(:,:) = flux_dn(:,(nlay+1):1:-1) + call write_broadband_field(input_file, flux_up, "lw_flux_up_vr", "LW flux up, vertically-reversed") + call write_broadband_field(input_file, flux_dn, "lw_flux_dn_vr", "LW flux dn, vertically-reversed") + + p_lay (:,:) = p_lay (:, nlay :1:-1) + t_lay (:,:) = t_lay (:, nlay :1:-1) + p_lev (:,:) = p_lev (:,(nlay+1):1:-1) + t_lev (:,:) = t_lev (:,(nlay+1):1:-1) + top_at_1 = .not. top_at_1 + end subroutine lw_clear_sky_vr + ! ---------------------------------------------------------------------------- + ! + ! Clear-sky longwave fluxes, all info, three angles + ! Two variations on a scattering calculation but using purely emitting/absorbing layers + ! The first uses rescaling and should agree with _1scl answers + ! The second uses the two-stream solver + ! + subroutine lw_clear_sky_2str + call stop_on_err(k_dist%gas_optics(p_lay, p_lev, & + t_lay, sfc_t, & + gas_concs, & + atmos, & + lw_sources, & + tlev=t_lev)) + call stop_on_err(rte_lw(atmos, top_at_1, & + lw_sources, & + sfc_emis, & + fluxes)) + call write_broadband_field(input_file, flux_up, "lw_flux_up_1rescl", "LW flux up, clear-sky _1rescl") + call write_broadband_field(input_file, flux_dn, "lw_flux_dn_1rescl", "LW flux dn, clear-sky _1rescl") + + call stop_on_err(rte_lw(atmos, top_at_1, & + lw_sources, & + sfc_emis, & + fluxes, use_2stream=.true.)) + call write_broadband_field(input_file, flux_up, "lw_flux_up_2str", "LW flux up, clear-sky _2str") + call write_broadband_field(input_file, flux_dn, "lw_flux_dn_2str", "LW flux dn, clear-sky _2str") + + end subroutine lw_clear_sky_2str + ! ---------------------------------------------------------------------------- + ! + ! Tests for incrementing optical properties + ! + subroutine lw_clear_sky_incr + type(ty_optical_props_1scl) :: atmos2 + type(ty_optical_props_2str) :: atmos_2str + type(ty_optical_props_nstr) :: atmos_nstr + integer, parameter :: nmom = 4 + + call stop_on_err(k_dist%gas_optics(p_lay, p_lev, & + t_lay, sfc_t, & + gas_concs, & + atmos, & + lw_sources, & + tlev = t_lev)) + + ! Two sets of optical properties, each with half the original optical + ! depth, should produce the same fluxes + call stop_on_err(atmos2%alloc_1scl(ncol, nlay, atmos)) + atmos%tau (:,:,:) = atmos%tau(:,:,:) * 0.5_wp + atmos2%tau(:,:,:) = atmos%tau(:,:,:) + call stop_on_err(atmos2%increment(atmos)) + + call stop_on_err(rte_lw(atmos, top_at_1, & + lw_sources, & + sfc_emis, & + fluxes)) + call write_broadband_field(input_file, flux_up, "lw_flux_up_inc_1scl_with_1scl", "LW flux up, incr. 1scl with 1scl") + call write_broadband_field(input_file, flux_dn, "lw_flux_dn_inc_1scl_with_1scl", "LW flux dn, incr. 1scl with 1scl") + + ! + ! Create a transparent set of 2-stream optical properties + ! + call stop_on_err(atmos_2str%alloc_2str(ncol, nlay, k_dist)) + atmos_2str%tau(:,:,:) = 0._wp + atmos_2str%ssa(:,:,:) = 0._wp + atmos_2str%g (:,:,:) = 0._wp + ! + ! Add a transparent 2-stream atmosphere to no-scattering + ! + call stop_on_err(atmos_2str%increment(atmos)) + call stop_on_err(rte_lw(atmos, top_at_1, & + lw_sources, & + sfc_emis, & + fluxes)) + call write_broadband_field(input_file, flux_up, "lw_flux_up_inc_1scl_with_2str", "LW flux up, incr. 1scl with 2str") + call write_broadband_field(input_file, flux_dn, "lw_flux_dn_inc_1scl_with_2str", "LW flux dn, incr. 1scl with 2str") + + ! + ! Create a transparent set of n-stream optical properties + ! + call stop_on_err(atmos_nstr%alloc_nstr(nmom, ncol, nlay, k_dist)) + atmos_nstr%tau(:,:,:) = 0._wp + atmos_nstr%ssa(:,:,:) = 0._wp + atmos_nstr%p (:,:,:,:) = 0._wp + ! + ! Add a transparent n-stream atmosphere to no-scattering + ! + call stop_on_err(atmos_nstr%increment(atmos)) + call stop_on_err(rte_lw(atmos, top_at_1, & + lw_sources, & + sfc_emis, & + fluxes)) + call write_broadband_field(input_file, flux_up, "lw_flux_up_inc_1scl_with_nstr", "LW flux up, incr. 1scl with nstr") + call write_broadband_field(input_file, flux_dn, "lw_flux_dn_inc_1scl_with_nstr", "LW flux dn, incr. 1scl with nstr") + + end subroutine lw_clear_sky_incr + ! ---------------------------------------------------------------------------- + ! + ! Shortwave tests + ! + ! ---------------------------------------------------------------------------- + ! + ! Shortwave - default + ! + subroutine sw_clear_sky_default + real(wp), dimension(ncol, ngpt) :: rfmip_tsi_scale + real(wp) :: rrtmgp_tsi + type(ty_optical_props_2str) :: atmos2 + + call stop_on_err(k_dist%gas_optics(p_lay, p_lev, & + t_lay, & + gas_concs, & + atmos, & + toa_flux)) + ! + ! Scaling factor for column dependence of TSI in RFMIP + ! + rrtmgp_tsi = sum(toa_flux(1,:)) + rfmip_tsi_scale(:,:) = spread(tsi_3d(:,1)/rrtmgp_tsi, dim=2, ncopies=ngpt) + toa_flux(:,:) = toa_flux(:,:) * rfmip_tsi_scale(:,:) + + call stop_on_err(rte_sw(atmos, top_at_1, & + mu0, toa_flux, & + sfc_alb_dir, sfc_alb_dif, & + fluxes)) + ! + ! Fluxes were computed for all columns; here mask where sun is down + ! + where(spread(.not. sun_up, dim=2, ncopies=nlay+1)) + flux_up = 0._wp + flux_dn = 0._wp + end where + call write_broadband_field(input_file, flux_up, "sw_flux_up", "SW flux up") + call write_broadband_field(input_file, flux_dn, "sw_flux_dn", "SW flux dn") + + ! + ! Two sets of optical properties, each with half the original optical + ! depth, should produce the same fluxes + ! + call stop_on_err(atmos2%alloc_2str(ncol, nlay, atmos)) + atmos%tau(:,:,:) = atmos%tau(:,:,:) * 0.5_wp + atmos2%tau(:,:,:) = atmos%tau(:,:,:) + select type(atmos) + class is (ty_optical_props_2str) + atmos2%ssa(:,:,:) = atmos%ssa(:,:,:) + atmos2%g (:,:,:) = atmos%g (:,:,:) + class default + call stop_on_err("rte_rrtmgp_atmos: Don't recognize the kind of optical properties ") + end select + call stop_on_err(atmos2%increment(atmos)) + call stop_on_err(rte_sw(atmos, top_at_1, & + mu0, toa_flux, & + sfc_alb_dir, sfc_alb_dif, & + fluxes)) + where(spread(.not. sun_up, dim=2, ncopies=nlay+1)) + flux_up = 0._wp + flux_dn = 0._wp + end where + call write_broadband_field(input_file, flux_up, "sw_flux_up_incr", "SW flux up, incr. 2str with 2str") + call write_broadband_field(input_file, flux_dn, "sw_flux_dn_incr", "SW flux dn, incr. 2str with 2str") + end subroutine sw_clear_sky_default + ! ---------------------------------------------------------------------------- + ! + ! Shortwave - Set the total solar irradiance + ! + subroutine sw_clear_sky_tsi + real(wp), dimension(ncol, ngpt) :: rfmip_tsi_scale + real(wp) :: rrtmgp_tsi + real(wp), parameter :: tsi_scale = 0.5_wp + + call stop_on_err(k_dist%gas_optics(p_lay, p_lev, & + t_lay, & + gas_concs, & + atmos, & + toa_flux)) + ! + ! Scaling factor for column dependence of TSI in RFMIP + ! + rrtmgp_tsi = sum(toa_flux(1,:)) + rfmip_tsi_scale(:,:) = spread(tsi_3d(:,1)/rrtmgp_tsi, dim=2, ncopies=ngpt) + + ! Set TSI to half the default + call stop_on_err(k_dist%set_tsi(tsi_scale*rrtmgp_tsi)) + ! Redo gas optics + call stop_on_err(k_dist%gas_optics(p_lay, p_lev, & + t_lay, & + gas_concs, & + atmos, & + toa_flux)) + toa_flux(:,:) = toa_flux(:,:) * rfmip_tsi_scale(:,:) + call stop_on_err(rte_sw(atmos, top_at_1, & + mu0, toa_flux, & + sfc_alb_dir, sfc_alb_dif, & + fluxes)) + ! + ! Fluxes were computed for all columns; here mask where sun is down + ! + where(spread(.not. sun_up, dim=2, ncopies=nlay+1)) + flux_up = 0._wp + flux_dn = 0._wp + end where + call write_broadband_field(input_file, flux_up /tsi_scale, "sw_flux_up_tsi", "SW flux up, reset TSI") + call write_broadband_field(input_file, flux_dn /tsi_scale, "sw_flux_dn_tsi", "SW flux dn, reset TSI") + + call stop_on_err(k_dist%set_tsi(2.0_wp * sum(toa_flux(1,:)))) + end subroutine sw_clear_sky_tsi + ! ---------------------------------------------------------------------------- + ! + ! Shortwave - vertically reversed + ! + subroutine sw_clear_sky_vr + real(wp), dimension(ncol, ngpt) :: rfmip_tsi_scale + real(wp) :: rrtmgp_tsi + type(ty_gas_concs) :: gas_concs_vr + integer :: i + real(wp), dimension(ncol,nlay) :: vmr + character(32), & + dimension(gas_concs%get_num_gases()) & + :: gc_gas_names + + ! + ! Reverse the orientation of the problem + ! + p_lay (:,:) = p_lay (:, nlay :1:-1) + t_lay (:,:) = t_lay (:, nlay :1:-1) + p_lev (:,:) = p_lev (:,(nlay+1):1:-1) + top_at_1 = .not. top_at_1 + ! + ! No direct access to gas concentrations so use the classes + ! This also tests otherwise uncovered routines for ty_gas_concs + ! + gc_gas_names(:) = gas_concs%get_gas_names() + call stop_on_err(gas_concs_vr%init(gc_gas_names(:))) + do i = 1, gas_concs%get_num_gases() + call stop_on_err(gas_concs%get_vmr(gc_gas_names(i), vmr)) + vmr(:,:) = vmr(:,nlay:1:-1) + call stop_on_err(gas_concs_vr%set_vmr(gc_gas_names(i), vmr)) + end do + + call stop_on_err(k_dist%gas_optics(p_lay, p_lev, & + t_lay, & + gas_concs_vr, & + atmos, & + toa_flux)) + ! + ! Scaling factor for column dependence of TSI in RFMIP + ! + rrtmgp_tsi = sum(toa_flux(1,:)) + rfmip_tsi_scale(:,:) = spread(tsi_3d(:,1)/rrtmgp_tsi, dim=2, ncopies=ngpt) + toa_flux(:,:) = toa_flux(:,:) * rfmip_tsi_scale(:,:) + + call stop_on_err(rte_sw(atmos, top_at_1, & + mu0, toa_flux, & + sfc_alb_dir, sfc_alb_dif, & + fluxes)) + ! + ! Fluxes were computed for all columns; here mask where sun is down + ! + where(spread(.not. sun_up, dim=2, ncopies=nlay+1)) + flux_up = 0._wp + flux_dn = 0._wp + end where + flux_up (:,:) = flux_up (:,(nlay+1):1:-1) + flux_dn (:,:) = flux_dn (:,(nlay+1):1:-1) + call write_broadband_field(input_file, flux_up, "sw_flux_up_vr", "SW flux up, vertically-reversed") + call write_broadband_field(input_file, flux_dn, "sw_flux_dn_vr", "SW flux dn, vertically-reversed") + + p_lay (:,:) = p_lay (:, nlay :1:-1) + t_lay (:,:) = t_lay (:, nlay :1:-1) + p_lev (:,:) = p_lev (:,(nlay+1):1:-1) + top_at_1 = .not. top_at_1 +end subroutine sw_clear_sky_vr + ! ---------------------------------------------------------------------------- + subroutine make_optical_props_1scl + if(allocated(atmos )) deallocate(atmos) + allocate(ty_optical_props_1scl::atmos) + ! Clouds optical props are defined by band + ! + ! Allocate arrays for the optical properties themselves. + ! + select type(atmos) + class is (ty_optical_props_1scl) + call stop_on_err(atmos%alloc_1scl(ncol, nlay, k_dist)) + class default + call stop_on_err("rte_rrtmgp_atmos: Don't recognize the kind of optical properties ") + end select + end subroutine make_optical_props_1scl + ! ---------------------------------------------------------------------------- + subroutine make_optical_props_2str + if(allocated(atmos )) deallocate(atmos) + allocate(ty_optical_props_2str::atmos) + ! Clouds optical props are defined by band + ! + ! Allocate arrays for the optical properties themselves. + ! + select type(atmos) + class is (ty_optical_props_2str) + call stop_on_err(atmos%alloc_2str(ncol, nlay, k_dist)) + class default + call stop_on_err("rte_rrtmgp_atmos: Don't recognize the kind of optical properties ") + end select + end subroutine make_optical_props_2str + ! ---------------------------------------------------------------------------- +end program rte_clear_sky_regression diff --git a/tests/intel-codecov.sh b/tests/intel-codecov.sh new file mode 100644 index 000000000..7539c32db --- /dev/null +++ b/tests/intel-codecov.sh @@ -0,0 +1,75 @@ +#!/bin/bash +ulimit -s hard +export FC=ifort +export FCFLAGS="-m64 -prof-gen=srcpos -g -traceback -heap-arrays -assume realloc_lhs -extend-source 132" +# +# Intel specific - where will the profiling files be generated? +# +export RRTMGP_ROOT=`cd ..;pwd` +export RRTMGP_BUILD=${RRTMGP_ROOT}/build +export PROF_DIR=$PWD +# +# Environment variables for netCDF Fortran and C installations +# +export NFHOME=${HOME}/Applications/${FC} +export NCHOME=/opt/local + +# +# An Anaconda environent with modules needed for other python scripts +# (xarray, matplotlib, ...) +# +source activate pangeo +cd ${PROF_DIR} +rm -rf *.dyn pgopti.* CODE_COVERAGE.html CodeCoverage/ +# +# Build RTE+RRTMGP librarues +# +make -C ${RRTMGP_BUILD} clean +make -C ${RRTMGP_BUILD} -j 4 || exit 1 + +# +# Build and run RFMIP examples +# +cd ${RRTMGP_ROOT}/examples/rfmip-clear-sky || exit 1 +export FCFLAGS+=" -I${RRTMGP_BUILD} -I${NFHOME}/include" +export LDFLAGS+=" -L${RRTMGP_BUILD} -L${NFHOME}/lib -L${NCHOME}/lib -lrte -lrrtmgp -lnetcdff -lnetcdf" +make clean || exit 1 +make -j 4 || exit 1 +python ./stage_files.py +python ./run-rfmip-examples.py +python ./compare-to-reference.py --fail=7.e-4 +make clean +# +# Build and run all-sky examples +# +cd ${RRTMGP_ROOT}/examples/all-sky || exit 1 +make clean || exit 1 +make -j 4 || exit 1 +python ./run-allsky-example.py +python ./compare-to-reference.py +make clean +# +# Build and run regression tests +# +cd ${PROF_DIR} || exit 1 +make clean || exit 1 +make -j 4 || exit 1 +cp ${RRTMGP_ROOT}/examples/rfmip-clear-sky/multiple_input4MIPs_radiation_RFMIP_UColorado-RFMIP-1-2_none.nc test_atmospheres.nc +./clear_sky_regression test_atmospheres.nc ${RRTMGP_ROOT}/rrtmgp/data/rrtmgp-data-lw-g256-2018-12-04.nc +./clear_sky_regression test_atmospheres.nc ${RRTMGP_ROOT}/rrtmgp/data/rrtmgp-data-sw-g224-2018-12-04.nc +# Need to repeat for openacc-kernels + +# +# Merge +# +cd ${PROF_DIR} +profmerge -a +echo " +mo_ +~mo_simple_netcdf +~mo_rfmip_io +~mo_testing_io +~mo_garand_atmos_io +~mo_load" > intel_codecov_filter.txt +codecov -prj rte-rrtmgp -comp intel_codecov_filter.txt +rm intel_codecov_filter.txt diff --git a/tests/mo_testing_io.F90 b/tests/mo_testing_io.F90 new file mode 100644 index 000000000..8444cac56 --- /dev/null +++ b/tests/mo_testing_io.F90 @@ -0,0 +1,60 @@ +!-------------------------------------------------------------------------------------------------------------------- +module mo_testing_io + ! + ! RTE+RRTMGP modules + ! + use mo_rte_kind, only: wp + ! + ! NetCDF I/O routines, shared with other RTE+RRTMGP examples + ! + use mo_simple_netcdf, only: write_field, create_dim, create_var + use netcdf + implicit none + private + public :: write_broadband_field +contains + !-------------------------------------------------------------------------------------------------------------------- + subroutine write_broadband_field(fileName, field, field_name, field_description, col_dim_name, vert_dim_name) + ! + ! Write a field defined by column and some vertical dimension (lev or lay)) + ! + character(len=*), intent(in) :: fileName + real(wp), dimension(:,:), intent(in) :: field + character(len=*), intent(in) :: field_name, field_description + character(len=*), optional, & + intent(in) ::col_dim_name, vert_dim_name + ! ------------------- + integer :: ncid, varid, ncol, nlev + ! + ! Names of column (first) and vertical (second) dimension. + ! Because they are used in an array constuctor the need to have the same number of characters + ! + character(len=32) :: cdim, vdim + ! ------------------- + cdim = "site " + vdim = "level" + if(present(col_dim_name)) cdim = col_dim_name + if(present(vert_dim_name)) vdim = vert_dim_name + + if(nf90_open(trim(fileName), NF90_WRITE, ncid) /= NF90_NOERR) & + call stop_on_err("write_fluxes: can't open file " // trim(fileName)) + + ncol = size(field, dim=1) + nlev = size(field, dim=2) + + call create_dim(ncid, trim(cdim), ncol) + call create_var(ncid, trim(field_name), [cdim, vdim], [ncol, nlev]) + call stop_on_err(write_field(ncid, trim(field_name), field)) + ! + ! Adding descriptive text as an attribute means knowing the varid + ! + if(nf90_inq_varid(ncid, trim(field_name), varid) /= NF90_NOERR) & + call stop_on_err("Can't find variable " // trim(field_name)) + if(nf90_put_att(ncid, varid, "description", trim(field_description)) /= NF90_NOERR) & + call stop_on_err("Can't write 'description' attribute to variable " // trim(field_name)) + + ncid = nf90_close(ncid) + + end subroutine write_broadband_field + !-------------------------------------------------------------------------------------------------------------------- +end module mo_testing_io diff --git a/tests/validation-plots.py b/tests/validation-plots.py new file mode 100644 index 000000000..0b223cb0b --- /dev/null +++ b/tests/validation-plots.py @@ -0,0 +1,138 @@ +import numpy as np +import xarray as xr +import matplotlib.pyplot as plt +from matplotlib.backends.backend_pdf import PdfPages +import seaborn as sb +import colorcet as cc + +def mae(diff, col_dim): + # + # Mean absolute error + # + return(xr.ufuncs.fabs(diff).mean(dim=col_dim)) +def rms(diff, col_dim): + # + # Root mean square error + # + return(xr.ufuncs.sqrt(xr.ufuncs.square(diff).mean(dim=col_dim))) + +def make_comparison_plot(variants, labels, reference, vscale, col_dim="site", lay_dim="layer"): + # + # Make a plot comparing differences with respect to reference + # + if type(variants) is not list: + make_comparison_plot([variants], [labels], reference, vscale) + else: + for i in np.arange(0, len(variants)): + delta = variants[i] - reference + plt.plot(mae(delta, col_dim), + vscale, '-',\ + color=cols[i], label=labels[i]) + plt.plot(rms(delta, col_dim), + vscale, '--',\ + color=cols[i]) + # Reverse vertical ordering + plt.legend() + # Reverse vertical ordering + plt.ylim(vscale.max(), vscale.min()) + +def construct_lbl_esgf_name(var): + # + # For a given variable name, provide the OpenDAP URL for the LBLRM line-by-line results + # + prefix = "http://esgf3.dkrz.de/thredds/dodsC/cmip6/RFMIP/AER/LBLRTM-12-8/rad-irf/r1i1p1f1/Efx/" + return(prefix + var + "/gn/v20190514/" + var + "_Efx_LBLRTM-12-8_rad-irf_r1i1p1f1_gn.nc") + +######################################################################## +if __name__ == '__main__': + gp = xr.open_dataset("test_atmospheres.nc") + # + # Does the flux plus the Jacobian equal a calculation with perturbed surface temperature? + # + gp['lw_flux_up_from_deriv'] = gp.lw_flux_up + gp.lw_jaco_up + gp.lw_flux_up_from_deriv.attrs = {"description":"LW flux up, surface T+1K, computed from Jacobian"} + lbl = xr.open_mfdataset([construct_lbl_esgf_name(v) for v in ["rsd", "rsu", "rld", "rlu"]], + combine = "by_coords").sel(expt=0) + ######################################################################## + # + # The RFMIP cases are on an irregular pressure grid so we can't compute errors + # as a function of pressure + # Interpolate onto a uniform pressure level - every 10 hPa + # plev varies from site to site - could the interpolation be made more concise? + # + plevs = np.arange(0, lbl.plev.max(), 1000) + plevs[0] = lbl.plev.min().values + plev = xr.DataArray(plevs, dims=['plev'], coords={'plev': plevs}) + lbli = xr.concat([lbl.sel(site=i).reset_coords().swap_dims({"level":"plev"}).interp(plev=plev) for i in np.arange(0, lbl.site.size)], dim = 'site') + gpi = xr.concat([ gp.sel(site=i).rename({"pres_level":"plev"}).\ + reset_coords().swap_dims({"level":"plev"}).interp(plev=plev) for i in np.arange(0, gp.site.size)], dim = 'site') + + cols = cc.glasbey_dark + with PdfPages('validation-figures.pdf') as pdf: + ######################################################################## + # Longwave + ######################################################################## + # + # Accuracy - 3-angle and single-angle + # + variants = [[gpi.lw_flux_dn, gpi.lw_flux_dn_optang, gpi.lw_flux_dn_3ang, gpi.lw_flux_dn_2str], + [gpi.lw_flux_up, gpi.lw_flux_up_optang, gpi.lw_flux_up_3ang, gpi.lw_flux_up_2str], + [gpi.lw_flux_net, + gpi.lw_flux_dn_optang - gpi.lw_flux_up_optang, + gpi.lw_flux_dn_3ang - gpi.lw_flux_up_3ang, + gpi.lw_flux_dn_2str - gpi.lw_flux_up_2str]] + 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", "optimal-angle", "3-angle", "2-stream"], \ + reference = r, \ + vscale = plev/100.) + plt.ylabel("Pressure (Pa)") + plt.xlabel("Error (W/m2), solid=mean, dash=RMS") + plt.title(t) + pdf.savefig() + plt.close() + + # + # Variants on the default - should be near but not precisely 0 + # + # Level temperatures not provided + # Using scattering representations of optical properties + # to do a no-scattering calculation + # + variants = [[gpi.lw_flux_dn_notlev, gpi.lw_flux_dn_2str], + [gpi.lw_flux_up_notlev, gpi.lw_flux_up_2str]] + refs = [gpi.lw_flux_dn, gpi.lw_flux_up] + titles = ["Variants: LW down", "Variants: LW up"] + for v, r, t in zip(variants, refs, titles): + make_comparison_plot(v, \ + labels = ["no-tlev", "2str"], \ + reference = r, \ + vscale = plev/100.) + plt.ylabel("Pressure (Pa)") + plt.xlabel("Difference (W/m2), solid=mean, dash=RMS") + plt.title(t) + pdf.savefig() + plt.close() + ######################################################################## + # Shortwave + # These are in W/m2 so profiles with high run count for more. + # Would they be better scaled to incoming solar flux? + ######################################################################## + # + # Accuracy comparison + # + variants = [gpi.sw_flux_dn, gpi.sw_flux_up] + refs = [lbli.rsd, lbli.rsu] + titles = ["Accuracy: SW down", "Accuracy: SW up"] + for v, r, t in zip(variants, refs, titles): + make_comparison_plot(v, \ + labels = "default", \ + reference = r, \ + vscale = plev/100.) + plt.ylabel("Pressure (Pa)") + plt.xlabel("Error (W/m2), solid=mean, dash=RMS") + plt.title(t) + pdf.savefig() + plt.close() diff --git a/tests/verification.py b/tests/verification.py new file mode 100644 index 000000000..d1e792c41 --- /dev/null +++ b/tests/verification.py @@ -0,0 +1,77 @@ +import xarray as xr +import argparse +import sys + +def assert_equal(variants, reference): + # + # Computes difference between reference and each variant, reports if differences + # exceed a threshold + # + passed = True + if type(variants) is not list: + assert_equal([variants], reference) + else: + print('Using %s as reference:'%(reference.description)) + for v in variants: + diff = xr.ufuncs.fabs(v - reference) + print(' %s'%(v.description)) + if diff.max() > report_threshold: + print(' differs from reference by as much as %e'%(diff.max())) + passed = passed and diff.max() <= failure_threshold + print('') + return(passed) + + +######################################################################## +if __name__ == '__main__': + parser = argparse.ArgumentParser(description="Compares all-sky example output to file in reference directory") + parser.add_argument("--report_threshold", type=float, default=1.e-10, + help="Threshold for reporting differences") + parser.add_argument("--failure_threshold", type=float, default=1.e-5, + help="Threshold at which differences cause failure") + args = parser.parse_args() + report_threshold = args.report_threshold + failure_threshold = args.failure_threshold + + gp = xr.open_dataset("test_atmospheres.nc") + ######################################################################## + # + # Some variants we expect to be essentially equal to the default, so we'll check + # for those + ############################### + # + # Longwave + # + gp['lw_flux_net_from_updn'] = gp.lw_flux_dn - gp.lw_flux_up + gp.lw_flux_net_from_updn.attrs = {"description":"LW flux net, computed externally from dn-up"} + + passed = assert_equal([gp.lw_flux_dn_vr, gp.lw_flux_dn_jaco, gp.lw_flux_dn_subset], gp.lw_flux_dn) + passed = assert_equal([gp.lw_flux_up_vr, gp.lw_flux_up_jaco, gp.lw_flux_up_subset], gp.lw_flux_up) and passed + passed = assert_equal([gp.lw_flux_net, gp.lw_flux_net_2], gp.lw_flux_net_from_updn) and passed + # + # Does the flux plus the Jacobian equal a calculation with perturbed surface temperature? + # + gp['lw_flux_up_from_deriv'] = gp.lw_flux_up_jaco + gp.lw_jaco_up + gp.lw_flux_up_from_deriv.attrs = {"description":"LW flux up, surface T+1K, computed from Jacobian"} + passed = assert_equal(gp.lw_flux_up_from_deriv, gp.lw_flux_up_stp1) and passed + ############################### + # + # Shortwave + # + passed = assert_equal([gp.sw_flux_dn_vr, gp.sw_flux_dn_tsi], gp.sw_flux_dn) and passed + passed = assert_equal([gp.sw_flux_up_vr, gp.sw_flux_up_tsi], gp.sw_flux_up) and passed + + print('Incrementing') + passed = assert_equal([gp.lw_flux_dn_inc_1scl_with_1scl, gp.lw_flux_dn_inc_1scl_with_2str, gp.lw_flux_dn_inc_1scl_with_nstr], + gp.lw_flux_dn) and passed + passed = assert_equal([gp.lw_flux_up_inc_1scl_with_1scl, gp.lw_flux_up_inc_1scl_with_2str, gp.lw_flux_up_inc_1scl_with_nstr], + gp.lw_flux_up) and passed + # passed = assert_equal(gp.lw_flux_dn_inc_2str_with_1scl, gp.lw_flux_dn_2str) and passed + # passed = assert_equal(gp.lw_flux_up_inc_2str_with_1scl, gp.lw_flux_up_2str) and passed + passed = assert_equal([gp.sw_flux_dn_incr], + gp.sw_flux_dn) and passed + passed = assert_equal([gp.sw_flux_up_incr], + gp.sw_flux_up) and passed + + if not passed: print("Something is terribly, terribly wrong") + sys.exit(0) if passed else sys.exit(1)