From 22f3b92b8411d9f4e935f2295bd5c45cfe316f39 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Tue, 26 May 2020 21:29:59 -0400 Subject: [PATCH] Minor bug fixes, changes to continuous integration, input sanitizing made optional. --- .github/workflows/main.yml | 110 +++++++++++ .travis.yml | 177 ------------------ README.md | 1 + azure-pipelines.yml | 54 +++--- examples/all-sky/rrtmgp_allsky.F90 | 8 +- examples/rfmip-clear-sky/Makefile | 1 + examples/rfmip-clear-sky/stage_files.py | 3 +- extensions/cloud_optics/mo_cloud_optics.F90 | 17 +- extensions/mo_compute_bc.F90 | 22 +-- extensions/mo_fluxes_byband.F90 | 34 ++-- extensions/mo_heating_rates.F90 | 17 +- rrtmgp/Make.depends | 4 +- .../kernels-openacc/mo_gas_optics_kernels.F90 | 1 + rrtmgp/kernels/mo_gas_optics_kernels.F90 | 1 + .../mo_rrtmgp_util_reorder_kernels.F90 | 1 + rrtmgp/mo_gas_concentrations.F90 | 15 +- rrtmgp/mo_gas_optics_rrtmgp.F90 | 101 +++++----- rrtmgp/mo_rrtmgp_constants.F90 | 1 + rte/Make.depends | 8 +- rte/kernels-openacc/mo_rte_solver_kernels.F90 | 40 ++-- rte/kernels/mo_rte_solver_kernels.F90 | 10 +- rte/mo_fluxes.F90 | 65 +++---- rte/mo_optical_props.F90 | 10 +- rte/mo_rte_config.F90 | 45 +++++ rte/mo_rte_kind.F90 | 1 + rte/mo_rte_lw.F90 | 75 ++++---- rte/mo_rte_sw.F90 | 54 +++--- 27 files changed, 451 insertions(+), 425 deletions(-) create mode 100644 .github/workflows/main.yml delete mode 100644 .travis.yml create mode 100644 rte/mo_rte_config.F90 diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml new file mode 100644 index 000000000..7dcdfce66 --- /dev/null +++ b/.github/workflows/main.yml @@ -0,0 +1,110 @@ +name: Continuous Integration +on: [push, pull_request] + +jobs: + CI: + runs-on: ubuntu-18.04 + strategy: + matrix: + fortran-compiler: [gfortran-8, gfortran-9, ifort] + rte-kernels: [default, openacc] + env: + FC: ${{ matrix.fortran-compiler }} + NCHOME: /usr + NFHOME: /usr + # I would like to use this variable within the run scripts + # i.e. to source files and change FCFLAGS, but can't figure out how + USING_IFORT: ${{ contains(matrix.fortran-compiler, 'ifort') }} + steps: + # + # Checks-out repository under $GITHUB_WORKSPACE + # + - uses: actions/checkout@v2 + # + # Set up Python and dependencies + # + - name: Set up Python + uses: actions/setup-python@v2 + with: + python-version: 3.7 + - name: Setup conda + uses: s-weigand/setup-conda@v1 + with: + python-version: 3.7 + - name: Install python dependencies + run: conda install --yes urllib3 netcdf4 xarray dask scipy + # + # Install NetCDF library + # + - name: Install netcdf C library + run: sudo apt-get install libnetcdf-dev + # + # Intel compilers and libraries if needed + # https://software.intel.com/content/www/us/en/develop/articles/oneapi-repo-instructions.html#aptpkg + # + - name: Install Intel compilers libraries + if: contains(matrix.fortran-compiler, 'ifort') + run: | + wget https://apt.repos.intel.com/intel-gpg-keys/GPG-PUB-KEY-INTEL-SW-PRODUCTS-2023.PUB + sudo apt-key add GPG-PUB-KEY-INTEL-SW-PRODUCTS-2023.PUB + sudo add-apt-repository "deb https://apt.repos.intel.com/oneapi all main" + sudo apt-get update + sudo apt-get install intel-oneapi-common-licensing + sudo apt-get install intel-oneapi-common-vars + sudo apt-get install intel-oneapi-dev-utilities + sudo apt-get install intel-oneapi-ifort + # + # NetCDF FORTRAN library + # + - name: Build NetCDF Fortran library + run: | + if [[ -e /opt/intel/inteloneapi/setvars.sh ]]; then source /opt/intel/inteloneapi/setvars.sh; fi + git clone https://github.com/Unidata/netcdf-fortran.git --branch v4.4.4 + cd netcdf-fortran + ./configure --prefix=${NFHOME} F77=${FC} + make + sudo make install + - name: Stage RFMIP files + run: | + export RRTMGP_ROOT=${GITHUB_WORKSPACE} + cd ${RRTMGP_ROOT}/examples/rfmip-clear-sky + python ./stage_files.py + - name: Make library, examples, tests + # Compiler flags for gfortran 8 and 9. Over-ridden in run script if using ifort + env: + FCFLAGS: "-ffree-line-length-none -m64 -std=f2008 -march=native -fbounds-check -finit-real=nan -DUSE_CBOOL" + run: | + export RRTMGP_ROOT=${GITHUB_WORKSPACE} + export RTE_KERNELS=${{ matrix.rte-kernels }} + if [[ -e /opt/intel/inteloneapi/setvars.sh ]]; then + source /opt/intel/inteloneapi/setvars.sh; + export FCFLAGS="-m64 -g -traceback -heap-arrays -assume realloc_lhs -extend-source 132 -check bounds,uninit,pointers,stack -stand f08"; + fi + cd ${RRTMGP_ROOT} + ${FC} --version + make -C build -j 2 + make -C tests -j 1 + make -C examples/all-sky -j 2 + export RRTMGP_BUILD=${RRTMGP_ROOT}/build + make -C examples/rfmip-clear-sky -j 2 + - name: Run examples, tests + run: | + export RRTMGP_ROOT=${GITHUB_WORKSPACE} + if [[ -e /opt/intel/inteloneapi/setvars.sh ]]; then source /opt/intel/inteloneapi/setvars.sh; fi + cd ${RRTMGP_ROOT}/examples/rfmip-clear-sky + python ./run-rfmip-examples.py --block_size 8 + cd ${RRTMGP_ROOT}/examples/all-sky + 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 + ./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 + - name: Comparison + run: | + export RRTMGP_ROOT=${GITHUB_WORKSPACE} + 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 diff --git a/.travis.yml b/.travis.yml deleted file mode 100644 index f2cc6ad31..000000000 --- a/.travis.yml +++ /dev/null @@ -1,177 +0,0 @@ -sudo: required -dist: xenial -language: cpp -# Could run both sets of kernels in parallel but setup is the majority of the time -#env: -# - RTE_KERNELS="" -# - RTE_KERNELS="openacc" -matrix: - include: - - name: "PGI community edition " - compiler: gcc - env: - - FC=pgfortran - - FCFLAGS="-Mallocatable=03 -Mstandard -Mbounds -Mchkptr -Kieee -Mchkstk" - - TEST_PGI=yes - addons: - apt: - sources: - - ubuntu-toolchain-r-test - packages: - - gcc-8 - - gfortran-8 - - libnetcdf-dev - - name: "gcc-9" - compiler: gcc - env: - - FC=gfortran-9 - - FCFLAGS="-ffree-line-length-none -m64 -std=f2008 -march=native -fbounds-check -finit-real=nan -DUSE_CBOOL" - - TEST_PGI=no - - NCHOME=/usr/ - addons: - apt: - sources: - - ubuntu-toolchain-r-test - packages: - - gcc-9 - - gfortran-9 - - libnetcdf-dev - - name: "gcc-8" - compiler: gcc - env: - - FC=gfortran-8 - - FCFLAGS="-ffree-line-length-none -m64 -std=f2008 -march=native -fbounds-check -finit-real=nan -DUSE_CBOOL" - - TEST_PGI=no - - NCHOME=/usr/ - addons: - apt: - sources: - - ubuntu-toolchain-r-test - packages: - - gcc-8 - - gfortran-8 - - libnetcdf-dev -# cache: -# timeout: 1000 -# directories: -# - "$HOME/pgi/" -before_install: -# Install PGI Community edition if not present in cache -- if [[ "${TEST_PGI}" == "yes" ]]; then - wget -q 'https://raw.githubusercontent.com/nemequ/pgi-travis/master/install-pgi.sh'; - chmod +x install-pgi.sh; - ./install-pgi.sh --nvidia; - pgfortran --version; - echo "switch -pthread is" > .userrc; - echo "append(LDLIB1=-lpthread);" >> .userrc; - 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 -- make -- make install -- cd .. || exit 1 - -# Install CONDA -- wget http://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -O miniconda.sh -- chmod +x miniconda.sh -- ./miniconda.sh -b -- export PATH=/home/travis/miniconda3/bin:$PATH -- conda update --yes conda - -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 scipy - -before_script: -- export RRTMGP_ROOT=$PWD -- 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: -# -# Default kernels -# -- cd ${RRTMGP_ROOT} -# -# Build library -# -- make -C build clean || exit 1 -- make -C build || exit 1 -# -# Build and run RFMIP examples -# -- 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=7.e-4 -- make clean -# -# Build and run all-sky example -# -- cd ${RRTMGP_ROOT}/examples/all-sky || exit 1 -- make clean -- make || exit 1 -- python3 ./run-allsky-example.py -- python3 ./compare-to-reference.py -- 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 -# -# 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 ${RRTMGP_ROOT}/examples/rfmip-clear-sky || exit 1 -- make clean -- make || exit 1 -- python3 ./run-rfmip-examples.py -- python3 ./compare-to-reference.py --fail=7.e-4 -- make clean -# -# Build and run all-sky example -# -- cd ${RRTMGP_ROOT}/examples/all-sky || exit 1 -- make clean -- make || exit 1 -- python3 ./run-allsky-example.py -- python3 ./compare-to-reference.py -- 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 - on_failure: never diff --git a/README.md b/README.md index 486c9a182..13c42f8ea 100644 --- a/README.md +++ b/README.md @@ -13,6 +13,7 @@ Example programs and documentation are evolving - please see examples/ in the re 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. +4. A new module, `mo_rte_config`, contains two logical variables that indicate whether arguments to routines are to be checked for correct extents and/or valid values. These variables can be changed via calls to `rte_config_checks()`. Setting the values to `.false.` removes the checks. Invalid values may cause incorrect results, crashes, or other mayhem 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. diff --git a/azure-pipelines.yml b/azure-pipelines.yml index 48bb0b79b..166ff69e8 100644 --- a/azure-pipelines.yml +++ b/azure-pipelines.yml @@ -4,38 +4,40 @@ jobs: strategy: matrix: pgi_19_9_gpu: + compiler_base: pgi compiler_module: PGI/19.9 - PrgEnv: PrgEnv-pgi accel_module: cudatoolkit FCFLAGS: "-O3 -ta=tesla:cc60,cuda10.1 -Mallocatable=03 -gopt -Minline,reshape,maxsize:40" RTE_KERNELS: openacc + RUN_CMD: "srun -C gpu -A c15 -p cscsci" pgi_default_gpu: + compiler_base: pgi compiler_module: pgi - PrgEnv: PrgEnv-pgi accel_module: craype-accel-nvidia60 # Generic accelerator flag FCFLAGS: "-O3 -acc -Mallocatable=03 -gopt" RTE_KERNELS: openacc + RUN_CMD: "srun -C gpu -A c15 -p cscsci" pgi_19_10_cpu: + compiler_base: pgi 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_base: pgi 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 + cce-cpu-icon-production: + compiler_base: cray + compiler_module: ccce-icon/9.0.2-classic accel_module: - # Error checking flags - FCFLAGS: "-Mallocatable=03 -Mstandard -Mbounds -Mchkptr -Kieee -Mchkstk" + # Production flags for Icon model + FCFLAGS: "-hadd_paren -r am -Ktrap=divz,ovf,inv -hflex_mp=intolerant -hfp1 -hnoacc -O1,cache0" maxParallel: 2 workspace: @@ -46,31 +48,30 @@ jobs: set -e echo " - export PATH=$CRAY_BINUTILS_BIN:$PATH module load daint-gpu - module swap PrgEnv-cray $(PrgEnv) - module load cray-netcdf cray-hdf5 - module swap pgi $(compiler_module) + export PATH=$CRAY_BINUTILS_BIN:$PATH + module swap PrgEnv-cray PrgEnv-$(compiler_base) + module swap $(compiler_base) $(compiler_module) module load $(accel_module) - module load cray-python/3.6.5.7 + module load cray-netcdf cray-hdf5 export LD_LIBRARY_PATH=$CRAY_LD_LIBRARY_PATH:$LD_LIBRARY_PATH - - echo Environment: + echo Compiler Environment: module list - " > modules + echo LD_LIBRARY_PATH is: + echo $LD_LIBRARY_PATH + " > compiler_modules displayName: 'Create module environment' - script: | set -e - source modules - # This module will unload some of the build modules, so load the files separately + module load daint-gpu module load netcdf-python cd examples/rfmip-clear-sky python ./stage_files.py displayName: 'Stage files' - script: | set -e - source modules + source compiler_modules export RRTMGP_ROOT=$PWD export RRTMGP_BUILD=$PWD/build export FC=ftn @@ -86,20 +87,21 @@ jobs: displayName: 'Make' - script: | set -e - source modules + source compiler_modules + module load cray-python/3.6.5.7 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 + ${RUN_CMD} python ./run-rfmip-examples.py --block_size 1800 cd ${RRTMGP_ROOT}/examples/all-sky - srun -C gpu -A c15 -p cscsci python ./run-allsky-example.py + ${RUN_CMD} 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 + ${RUN_CMD} ./clear_sky_regression test_atmospheres.nc ${RRTMGP_ROOT}/rrtmgp/data/rrtmgp-data-lw-g256-2018-12-04.nc + ${RUN_CMD} ./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 + module load daint-gpu export RRTMGP_ROOT=$PWD # This module will unload some of the build modules, so do the checks separately module load netcdf-python diff --git a/examples/all-sky/rrtmgp_allsky.F90 b/examples/all-sky/rrtmgp_allsky.F90 index a87ef084b..a1c3ea666 100755 --- a/examples/all-sky/rrtmgp_allsky.F90 +++ b/examples/all-sky/rrtmgp_allsky.F90 @@ -31,7 +31,7 @@ subroutine vmr_2d_to_1d(gas_concs, gas_concs_garand, name, sz1, sz2) end subroutine vmr_2d_to_1d ! ---------------------------------------------------------------------------------- program rte_rrtmgp_clouds - use mo_rte_kind, only: wp + use mo_rte_kind, only: wp, i8 use mo_optical_props, only: ty_optical_props, & ty_optical_props_arry, ty_optical_props_1scl, ty_optical_props_2str use mo_gas_optics_rrtmgp, only: ty_gas_optics_rrtmgp @@ -112,9 +112,9 @@ program rte_rrtmgp_clouds ! ! Timing variables ! - integer(kind=8) :: start, finish, start_all, finish_all, clock_rate - real(wp) :: avg - integer(kind=8), allocatable :: elapsed(:) + integer(kind=i8) :: start, finish, start_all, finish_all, clock_rate + real(wp) :: avg + integer(kind=i8), allocatable :: elapsed(:) !$omp threadprivate( lw_sources, toa_flux, flux_up, flux_dn, flux_dir ) ! ---------------------------------------------------------------------------------- ! Code diff --git a/examples/rfmip-clear-sky/Makefile b/examples/rfmip-clear-sky/Makefile index 1a47d561d..74d59c44c 100644 --- a/examples/rfmip-clear-sky/Makefile +++ b/examples/rfmip-clear-sky/Makefile @@ -2,6 +2,7 @@ # Here set variables RRTMGP_BUILD, NCHOME, NFHOME, TIME_DIR (for GPTL) # or have those variables set in the environment # +RRTMGP_BUILD=$(RRTMGP_ROOT)/build -include Makefile.libs -include $(RRTMGP_BUILD)/Makefile.conf # diff --git a/examples/rfmip-clear-sky/stage_files.py b/examples/rfmip-clear-sky/stage_files.py index 0009bc48e..72ed54475 100755 --- a/examples/rfmip-clear-sky/stage_files.py +++ b/examples/rfmip-clear-sky/stage_files.py @@ -2,6 +2,7 @@ # # This script downloads and creates files needed for the RFMIP off-line test cases # +import sys import os import subprocess import glob @@ -30,4 +31,4 @@ urllib.request.urlretrieve(conds_url, conds_file) print("Dowloading scripts for generating output templates") urllib.request.urlretrieve(templ_scr_url, templ_scr) -subprocess.run(["python3", templ_scr, "--source_id", "RTE-RRTMGP-181204"]) +subprocess.run([sys.executable, templ_scr, "--source_id", "RTE-RRTMGP-181204"]) diff --git a/extensions/cloud_optics/mo_cloud_optics.F90 b/extensions/cloud_optics/mo_cloud_optics.F90 index f8eafafac..bd812ed88 100644 --- a/extensions/cloud_optics/mo_cloud_optics.F90 +++ b/extensions/cloud_optics/mo_cloud_optics.F90 @@ -20,6 +20,7 @@ module mo_cloud_optics use mo_rte_kind, only: wp, wl + use mo_rte_config, only: check_values, check_extents use mo_rte_util_array,only: any_vals_less_than, any_vals_outside, extents_are use mo_optical_props, only: ty_optical_props, & ty_optical_props_arry, & @@ -433,13 +434,15 @@ function cloud_optics(this, & ! ! Particle size, liquid/ice water paths ! - if(any_vals_outside(reliq, liqmsk, this%radliq_lwr, this%radliq_upr)) & - error_msg = 'cloud optics: liquid effective radius is out of bounds' - if(any_vals_outside(reice, icemsk, this%radice_lwr, this%radice_upr)) & - error_msg = 'cloud optics: ice effective radius is out of bounds' - if(any_vals_less_than(clwp, liqmsk, 0._wp) .or. any_vals_less_than(ciwp, icemsk, 0._wp)) & - error_msg = 'cloud optics: negative clwp or ciwp where clouds are supposed to be' - if(error_msg == "") then + if(check_values) then + if(any_vals_outside(reliq, liqmsk, this%radliq_lwr, this%radliq_upr)) & + error_msg = 'cloud optics: liquid effective radius is out of bounds' + if(any_vals_outside(reice, icemsk, this%radice_lwr, this%radice_upr)) & + error_msg = 'cloud optics: ice effective radius is out of bounds' + if(any_vals_less_than(clwp, liqmsk, 0._wp) .or. any_vals_less_than(ciwp, icemsk, 0._wp)) & + error_msg = 'cloud optics: negative clwp or ciwp where clouds are supposed to be' + if(error_msg == "") then + end if ! ! ! ---------------------------------------- diff --git a/extensions/mo_compute_bc.F90 b/extensions/mo_compute_bc.F90 index c44a338b1..ad692f1a0 100644 --- a/extensions/mo_compute_bc.F90 +++ b/extensions/mo_compute_bc.F90 @@ -20,6 +20,7 @@ module mo_compute_bc ! The boundary condition is on diffuse flux in the LW and direct flux in the SW ! ------------------------------------------------------------------------------------------------- use mo_rte_kind, only: wp, wl + use mo_rte_config, only: check_extents use mo_rte_util_array, only: extents_are use mo_source_functions, only: ty_source_func_lw use mo_gas_concentrations, only: ty_gas_concs @@ -89,19 +90,16 @@ function compute_bc(k_dist, & ncol = size(play, dim=1) nlay = size(play, dim=2) ngpt = k_dist%get_ngpt() - if(.not. extents_are(plev, ncol, nlay+1)) then - error_msg = "compute_bc: array plev has wrong dimensions" - return - end if - if(.not. extents_are(tlay, ncol, nlay )) then - error_msg = "compute_bc: array tlay has wrong dimensions" - return - end if - if(present(mu0)) then - if(size(mu0) /= ncol) then - error_msg = "compute bc: array mu0 has wrong dimensions" - return + if(check_extents) then + if(.not. extents_are(plev, ncol, nlay+1)) & + error_msg = "compute_bc: array plev has wrong dimensions" + if(.not. extents_are(tlay, ncol, nlay )) & + error_msg = "compute_bc: array tlay has wrong dimensions" + if(present(mu0)) then + if(size(mu0) /= ncol) & + error_msg = "compute bc: array mu0 has wrong dimensions" end if + if(error_msg /= "") return end if ! diff --git a/extensions/mo_fluxes_byband.F90 b/extensions/mo_fluxes_byband.F90 index be08a2d72..c330f63a6 100644 --- a/extensions/mo_fluxes_byband.F90 +++ b/extensions/mo_fluxes_byband.F90 @@ -17,6 +17,7 @@ ! module mo_fluxes_byband use mo_rte_kind, only: wp + use mo_rte_config, only: check_extents use mo_rte_util_array,only: extents_are use mo_fluxes, only: ty_fluxes, ty_fluxes_broadband use mo_optical_props, only: ty_optical_props @@ -69,29 +70,24 @@ function reduce_byband(this, gpt_flux_up, gpt_flux_dn, spectral_disc, top_at_1, end if ! Check sizes of output arrays - if(associated(this%bnd_flux_up)) then - if(.not. extents_are(this%bnd_flux_up, ncol, nlev, nbnd)) then - error_msg = "reduce: bnd_flux_up array incorrectly sized (can't compute net flux either)" - return + if(check_extents) then + if(associated(this%bnd_flux_up)) then + if(.not. extents_are(this%bnd_flux_up, ncol, nlev, nbnd)) & + error_msg = "reduce: bnd_flux_up array incorrectly sized (can't compute net flux either)" end if - end if - if(associated(this%bnd_flux_dn)) then - if(.not. extents_are(this%bnd_flux_dn, ncol, nlev, nbnd)) then - error_msg = "reduce: bnd_flux_dn array incorrectly sized (can't compute net flux either)" - return + if(associated(this%bnd_flux_dn)) then + if(.not. extents_are(this%bnd_flux_dn, ncol, nlev, nbnd)) & + error_msg = "reduce: bnd_flux_dn array incorrectly sized (can't compute net flux either)" end if - end if - if(associated(this%bnd_flux_dn_dir)) then - if(.not. extents_are(this%bnd_flux_dn_dir, ncol, nlev, nbnd)) then - error_msg = "reduce: bnd_flux_dn_dir array incorrectly sized" - return + if(associated(this%bnd_flux_dn_dir)) then + if(.not. extents_are(this%bnd_flux_dn_dir, ncol, nlev, nbnd)) & + error_msg = "reduce: bnd_flux_dn_dir array incorrectly sized" end if - end if - if(associated(this%bnd_flux_net)) then - if(.not. extents_are(this%bnd_flux_net, ncol, nlev, nbnd)) then - error_msg = "reduce: bnd_flux_net array incorrectly sized (can't compute net flux either)" - return + if(associated(this%bnd_flux_net)) then + if(.not. extents_are(this%bnd_flux_net, ncol, nlev, nbnd)) & + error_msg = "reduce: bnd_flux_net array incorrectly sized (can't compute net flux either)" end if + if(error_msg /= "") return end if ! ! Self-consistency -- shouldn't be asking for direct beam flux if it isn't supplied diff --git a/extensions/mo_heating_rates.F90 b/extensions/mo_heating_rates.F90 index 19736552e..5ee14daba 100644 --- a/extensions/mo_heating_rates.F90 +++ b/extensions/mo_heating_rates.F90 @@ -13,6 +13,7 @@ module mo_heating_rates use mo_rte_kind, only: wp, wl + use mo_rte_config, only: check_extents use mo_rte_util_array, only: extents_are use mo_rrtmgp_constants, only: cp_dry, grav ! Only needed for heating rate calculation implicit none @@ -34,13 +35,15 @@ function compute_heating_rate(flux_up, flux_dn, plev, heating_rate) result(error ncol = size(flux_up, 1) nlay = size(flux_up, 2) - 1 - if(.not. extents_are(flux_dn, ncol, nlay+1)) & - error_msg = "heating_rate: flux_dn array inconsistently sized." - if(.not. extents_are(plev, ncol, nlay+1)) & - error_msg = "heating_rate: plev array inconsistently sized." - if(.not. extents_are(heating_rate, ncol, nlay)) & - error_msg = "heating_rate: heating_rate array inconsistently sized." - if(error_msg /= "") return + if(check_extents) then + if(.not. extents_are(flux_dn, ncol, nlay+1)) & + error_msg = "heating_rate: flux_dn array inconsistently sized." + if(.not. extents_are(plev, ncol, nlay+1)) & + error_msg = "heating_rate: plev array inconsistently sized." + if(.not. extents_are(heating_rate, ncol, nlay)) & + error_msg = "heating_rate: heating_rate array inconsistently sized." + if(error_msg /= "") return + end if do ilay = 1, nlay heating_rate(1:ncol,ilay) = (flux_up(1:ncol,ilay+1) - flux_up(1:ncol,ilay) - & diff --git a/rrtmgp/Make.depends b/rrtmgp/Make.depends index d051ee0e6..1395b71c5 100644 --- a/rrtmgp/Make.depends +++ b/rrtmgp/Make.depends @@ -24,13 +24,13 @@ mo_rrtmgp_util_reorder.o: mo_rte_kind.o mo_rrtmgp_util_reorder_kernels.o mo_rrtm # # Gas concentrations # -mo_gas_concentrations.o: mo_rte_kind.o mo_rte_util_array.o mo_rrtmgp_util_string.o mo_gas_concentrations.F90 +mo_gas_concentrations.o: mo_rte_kind.o mo_rte_config.o mo_rte_util_array.o mo_rrtmgp_util_string.o mo_gas_concentrations.F90 # # Gas optics # mo_gas_optics_kernels.o: mo_rte_kind.o mo_gas_optics_kernels.F90 -mo_gas_optics.o: mo_rte_kind.o mo_gas_concentrations.o \ +mo_gas_optics.o: mo_rte_kind.o mo_rte_config.o mo_gas_concentrations.o \ mo_optical_props.o mo_source_functions.o \ mo_gas_optics.F90 diff --git a/rrtmgp/kernels-openacc/mo_gas_optics_kernels.F90 b/rrtmgp/kernels-openacc/mo_gas_optics_kernels.F90 index e6e7996ed..bbfc2ba36 100644 --- a/rrtmgp/kernels-openacc/mo_gas_optics_kernels.F90 +++ b/rrtmgp/kernels-openacc/mo_gas_optics_kernels.F90 @@ -17,6 +17,7 @@ module mo_gas_optics_kernels use mo_rte_kind, only: wp, wl implicit none + public contains ! -------------------------------------------------------------------------------------- ! Compute interpolation coefficients diff --git a/rrtmgp/kernels/mo_gas_optics_kernels.F90 b/rrtmgp/kernels/mo_gas_optics_kernels.F90 index 1c11ab36c..4f064aa6c 100644 --- a/rrtmgp/kernels/mo_gas_optics_kernels.F90 +++ b/rrtmgp/kernels/mo_gas_optics_kernels.F90 @@ -17,6 +17,7 @@ module mo_gas_optics_kernels use mo_rte_kind, only : wp, wl implicit none + public contains ! -------------------------------------------------------------------------------------- ! Compute interpolation coefficients diff --git a/rrtmgp/kernels/mo_rrtmgp_util_reorder_kernels.F90 b/rrtmgp/kernels/mo_rrtmgp_util_reorder_kernels.F90 index a952ff221..9d5d850c2 100644 --- a/rrtmgp/kernels/mo_rrtmgp_util_reorder_kernels.F90 +++ b/rrtmgp/kernels/mo_rrtmgp_util_reorder_kernels.F90 @@ -16,6 +16,7 @@ module mo_rrtmgp_util_reorder_kernels use mo_rte_kind, only: wp implicit none + public contains ! ---------------------------------------------------------------------------- subroutine reorder_123x312_kernel(d1, d2, d3, array_in, array_out) & diff --git a/rrtmgp/mo_gas_concentrations.F90 b/rrtmgp/mo_gas_concentrations.F90 index a7df8616f..2c6e14879 100755 --- a/rrtmgp/mo_gas_concentrations.F90 +++ b/rrtmgp/mo_gas_concentrations.F90 @@ -31,6 +31,7 @@ module mo_gas_concentrations use mo_rte_kind, only: wp + use mo_rte_config, only: check_values use mo_rrtmgp_util_string, only: lower_case use mo_rte_util_array, only: any_vals_outside implicit none @@ -167,9 +168,10 @@ function set_vmr_1d(this, gas, w) result(error_msg) ! --------- error_msg = '' - if (any_vals_outside(w, 0._wp, 1._wp)) then - error_msg = 'ty_gas_concs%set_vmr: concentrations should be >= 0, <= 1' - endif + if (check_values) then + if (any_vals_outside(w, 0._wp, 1._wp)) & + error_msg = 'ty_gas_concs%set_vmr: concentrations should be >= 0, <= 1' + end if if(this%nlay > 0) then if(size(w) /= this%nlay) error_msg = 'ty_gas_concs%set_vmr: different dimension (nlay)' else @@ -217,9 +219,10 @@ function set_vmr_2d(this, gas, w) result(error_msg) ! --------- error_msg = '' - if (any_vals_outside(w, 0._wp, 1._wp)) then - error_msg = 'ty_gas_concs%set_vmr: concentrations should be >= 0, <= 1' - endif + if (check_values) then + if (any_vals_outside(w, 0._wp, 1._wp)) & + error_msg = 'ty_gas_concs%set_vmr: concentrations should be >= 0, <= 1' + end if if(this%ncol > 0 .and. size(w, 1) /= this%ncol) then error_msg = 'ty_gas_concs%set_vmr: different dimension (ncol)' diff --git a/rrtmgp/mo_gas_optics_rrtmgp.F90 b/rrtmgp/mo_gas_optics_rrtmgp.F90 index f2e8f2d1c..3982f0978 100644 --- a/rrtmgp/mo_gas_optics_rrtmgp.F90 +++ b/rrtmgp/mo_gas_optics_rrtmgp.F90 @@ -22,14 +22,14 @@ ! ------------------------------------------------------------------------------------------------- module mo_gas_optics_rrtmgp use mo_rte_kind, only: wp, wl - use mo_rrtmgp_constants, only: avogad, m_dry, m_h2o, grav + use mo_rte_config, only: check_extents, check_values use mo_rte_util_array, only: zero_array, any_vals_less_than, any_vals_outside, extents_are use mo_optical_props, only: ty_optical_props use mo_source_functions, only: ty_source_func_lw use mo_gas_optics_kernels, only: interpolation, & compute_tau_absorption, compute_tau_rayleigh, compute_Planck_source, & combine_and_reorder_2str, combine_and_reorder_nstr - + use mo_rrtmgp_constants, only: avogad, m_dry, m_h2o, grav use mo_rrtmgp_util_string, only: lower_case, string_in_array, string_loc_in_array use mo_gas_concentrations, only: ty_gas_concs use mo_optical_props, only: ty_optical_props_arry, ty_optical_props_1scl, ty_optical_props_2str, ty_optical_props_nstr @@ -270,27 +270,31 @@ function gas_optics_int(this, & ! External source -- check arrays sizes and values ! input data sizes and values ! - !$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)) & - error_msg = "gas_optics(): array tsfc has values outside range" - if(error_msg /= '') return + !$acc enter data copyin(tsfc, tlev) ! Should be fine even if tlev is not supplied - 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)) & - error_msg = "gas_optics(): array tlev has values outside range" - if(error_msg /= '') return + if(check_extents) then + if(.not. extents_are(tsfc, ncol)) & + error_msg = "gas_optics(): array tsfc has wrong size" + if(present(tlev)) then + if(.not. extents_are(tlev, ncol, nlay+1)) & + error_msg = "gas_optics(): array tlev has wrong size" + end if + ! + ! output extents + ! + if(any([sources%get_ncol(), sources%get_nlay(), sources%get_ngpt()] /= [ncol, nlay, ngpt])) & + error_msg = "gas_optics%gas_optics: source function arrays inconsistently sized" end if + if(error_msg /= '') return - ! - ! output extents - ! - if(any([sources%get_ncol(), sources%get_nlay(), sources%get_ngpt()] /= [ncol, nlay, ngpt])) & - error_msg = "gas_optics%gas_optics: source function arrays inconsistently sized" + if(check_values) then + if(any_vals_outside(tsfc, this%temp_ref_min, this%temp_ref_max)) & + error_msg = "gas_optics(): array tsfc has values outside range" + if(present(tlev)) then + if(any_vals_outside(tlev, this%temp_ref_min, this%temp_ref_max)) & + error_msg = "gas_optics(): array tlev has values outside range" + end if + end if if(error_msg /= '') return ! @@ -376,8 +380,10 @@ function gas_optics_ext(this, & ! External source function is constant ! !$acc enter data create(toa_src) - if(.not. extents_are(toa_src, ncol, ngpt)) & - error_msg = "gas_optics(): array toa_src has wrong size" + if(check_extents) then + if(.not. extents_are(toa_src, ncol, ngpt)) & + error_msg = "gas_optics(): array toa_src has wrong size" + end if if(error_msg /= '') return !$acc parallel loop collapse(2) @@ -463,33 +469,38 @@ function compute_gas_taus(this, & ! Check input data sizes and values ! !$acc enter data copyin(play,plev,tlay) - if(.not. extents_are(play, ncol, nlay )) & - error_msg = "gas_optics(): array play has wrong size" - if(.not. extents_are(tlay, ncol, nlay )) & - 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(check_extents) then + if(.not. extents_are(play, ncol, nlay )) & + error_msg = "gas_optics(): array play has wrong size" + if(.not. extents_are(tlay, ncol, nlay )) & + 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(present(col_dry)) then + if(.not. extents_are(col_dry, ncol, nlay)) & + error_msg = "gas_optics(): array col_dry has wrong size" + end if + end if if(error_msg /= '') return - if(any_vals_outside(play, this%press_ref_min,this%press_ref_max)) & - error_msg = "gas_optics(): array play has values outside range" - if(any_vals_outside(plev, this%press_ref_min,this%press_ref_max)) & - error_msg = "gas_optics(): array plev has values outside range" - if(any_vals_outside(tlay, this%temp_ref_min, this%temp_ref_max)) & - error_msg = "gas_optics(): array tlay has values outside range" + if(check_values) then + if(any_vals_outside(play, this%press_ref_min,this%press_ref_max)) & + error_msg = "gas_optics(): array play has values outside range" + if(any_vals_outside(plev, this%press_ref_min,this%press_ref_max)) & + error_msg = "gas_optics(): array plev has values outside range" + if(any_vals_outside(tlay, this%temp_ref_min, this%temp_ref_max)) & + error_msg = "gas_optics(): array tlay has values outside range" + if(present(col_dry)) then + if(any_vals_less_than(col_dry, 0._wp)) & + error_msg = "gas_optics(): array col_dry has values outside range" + end if + end if if(error_msg /= '') return - if(present(col_dry)) then - if(.not. extents_are(col_dry, ncol, nlay)) & - error_msg = "gas_optics(): array col_dry has wrong size" - if(any_vals_less_than(col_dry, 0._wp)) & - error_msg = "gas_optics(): array col_dry has values outside range" - if(error_msg /= '') return - end if ! ---------------------------------------------------------- ngas = this%get_ngas() diff --git a/rrtmgp/mo_rrtmgp_constants.F90 b/rrtmgp/mo_rrtmgp_constants.F90 index 5c98aeb05..10293c60b 100644 --- a/rrtmgp/mo_rrtmgp_constants.F90 +++ b/rrtmgp/mo_rrtmgp_constants.F90 @@ -21,6 +21,7 @@ ! ------------------------------------------------------------------------------------------------- module mo_rrtmgp_constants use mo_rte_kind, only: wp + public ! ----------------------------------------- ! Physical constants, 2018 SI defintion of metric system diff --git a/rte/Make.depends b/rte/Make.depends index 8edcea3ac..47cb2e2df 100644 --- a/rte/Make.depends +++ b/rte/Make.depends @@ -1,5 +1,6 @@ RTE_SRC = \ mo_rte_kind.o \ + mo_rte_config.o \ mo_rte_util_array.o \ mo_optical_props_kernels.o \ mo_optical_props.o \ @@ -15,6 +16,9 @@ RTE_SRC = \ ################################## # # +mo_rte_config.o: mo_rte_config.F90 mo_rte_kind.o +# +# mo_rte_util_array.o: mo_rte_util_array.F90 mo_rte_kind.o # # Optical properties @@ -29,13 +33,14 @@ mo_source_functions.o: mo_rte_kind.o mo_optical_props.o mo_source_functions # Flux reduction # mo_fluxes_broadband_kernels.o : mo_rte_kind.o mo_fluxes_broadband_kernels.F90 -mo_fluxes.o: mo_rte_kind.o mo_fluxes_broadband_kernels.o mo_optical_props.o mo_rte_util_array.o mo_fluxes.F90 +mo_fluxes.o: mo_rte_kind.o mo_fluxes_broadband_kernels.o mo_rte_config.o mo_optical_props.o mo_rte_util_array.o mo_fluxes.F90 # # Radiative transfer # mo_rte_solver_kernels.o: mo_rte_kind.o mo_rte_solver_kernels.F90 mo_rte_lw.o: mo_rte_kind.o \ + mo_rte_config.o \ mo_rte_util_array.o \ mo_optical_props.o \ mo_source_functions.o \ @@ -44,6 +49,7 @@ mo_rte_lw.o: mo_rte_kind.o \ mo_rte_lw.F90 mo_rte_sw.o: mo_rte_kind.o \ + mo_rte_config.o \ mo_rte_util_array.o \ mo_optical_props.o \ mo_source_functions.o \ diff --git a/rte/kernels-openacc/mo_rte_solver_kernels.F90 b/rte/kernels-openacc/mo_rte_solver_kernels.F90 index f234e1500..775c57520 100644 --- a/rte/kernels-openacc/mo_rte_solver_kernels.F90 +++ b/rte/kernels-openacc/mo_rte_solver_kernels.F90 @@ -525,9 +525,11 @@ subroutine lw_source_noscat_stencil(ncol, nlay, ngpt, icol, ilay, igpt, ! is of order epsilon (smallest difference from 1. in working precision) ! Thanks to Peter Blossey ! - fact = merge((1._wp - trans(icol,ilay,igpt))/tau(icol,ilay,igpt) - trans(icol,ilay,igpt), & - tau(icol,ilay,igpt) * ( 0.5_wp - 1._wp/3._wp*tau(icol,ilay,igpt) ), & - tau(icol,ilay,igpt) > tau_thresh) + if(tau(icol,ilay,igpt) > tau_thresh) then + fact = (1._wp - trans(icol,ilay,igpt))/tau(icol,ilay,igpt) - trans(icol,ilay,igpt) + else + fact = tau(icol, ilay,igpt) * (0.5_wp - 1._wp/3._wp*tau(icol,ilay,igpt)) + end if ! ! Equation below is developed in Clough et al., 1992, doi:10.1029/92JD01419, Eq 13 ! @@ -1234,13 +1236,13 @@ end subroutine apply_BC_0 ! b) adds adustment term based on cloud properties (lw_transport_1rescl) ! adustment terms is computed based on solution of the Tang equations ! for "linear-in-tau" internal source (not in the paper) -! +! ! Attention: ! use must prceompute scaling before colling the function ! ! Implemented based on the paper ! Tang G, et al, 2018: https://doi.org/10.1175/JAS-D-18-0014.1 -! +! ! ------------------------------------------------------------------------------------------------- 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, & @@ -1354,12 +1356,12 @@ subroutine lw_solver_1rescl(ncol, nlay, ngpt, top_at_1, D, ! ! Transport ! - ! compute no-scattering fluxes + ! 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,& source_sfcJac, rad_up_Jac) - ! make adjustment + ! make adjustment call lw_transport_1rescl(ncol, nlay, ngpt, top_at_1, & tau_loc, trans, & sfc_albedo, source_dn, source_up, & @@ -1380,9 +1382,9 @@ end subroutine lw_solver_1rescl ! ! Similar to lw_solver_noscat_GaussQuad. ! It is main solver to use the Tang approximation for fluxes -! In addition to the no scattering input parameters the user must provide -! scattering related properties (ssa and g) that the solver uses to compute scaling -! +! In addition to the no scattering input parameters the user must provide +! scattering related properties (ssa and g) that the solver uses to compute scaling +! ! --------------------------------------------------------------- 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, & @@ -1498,7 +1500,7 @@ end subroutine lw_solver_1rescl_GaussQuad ! ! Computes Tang scaling of layer optical thickness and scaling parameter ! unsafe if ssa*g =1. -! +! ! --------------------------------------------------------------- pure subroutine scaling_1rescl(ncol, nlay, ngpt, tauLoc, scaling, tau, ssa, g) integer , intent(in) :: ncol @@ -1522,9 +1524,9 @@ pure subroutine scaling_1rescl(ncol, nlay, ngpt, tauLoc, scaling, tau, ssa, g) ssal = ssa(icol, ilay, igpt) wb = ssal*(1._wp - g(icol, ilay, igpt)) / 2._wp scaleTau = (1._wp - ssal + wb ) - + tauLoc(icol, ilay, igpt) = scaleTau * tau(icol, ilay, igpt) ! Eq.15 of the paper - ! + ! ! here scaling is used to store parameter wb/[1-w(1-b)] of Eq.21 of the Tang's paper ! actually it is in line of parameter rescaling defined in Eq.7 ! potentialy if g=ssa=1 then wb/scaleTau = NaN @@ -1540,7 +1542,7 @@ end subroutine scaling_1rescl ! ! Computes Tang scaling of layer optical thickness and scaling parameter ! Safe implementation -! +! ! --------------------------------------------------------------- pure subroutine scaling_1rescl_safe(ncol, nlay, ngpt, tauLoc, scaling, tau, ssa, g) integer , intent(in) :: ncol @@ -1564,9 +1566,9 @@ pure subroutine scaling_1rescl_safe(ncol, nlay, ngpt, tauLoc, scaling, tau, ssa, ssal = ssa(icol, ilay, igpt) wb = ssal*(1._wp - g(icol, ilay, igpt)) / 2._wp scaleTau = (1._wp - ssal + wb ) - + tauLoc(icol, ilay, igpt) = scaleTau * tau(icol, ilay, igpt) ! Eq.15 of the paper - ! + ! ! here scaling is used to store parameter wb/[1-w(1-b)] of Eq.21 of the Tang's paper ! actually it is in line of parameter rescaling defined in Eq.7 if (scaleTau < 1e-6_wp) then @@ -1602,7 +1604,7 @@ subroutine lw_transport_1rescl(ncol, nlay, ngpt, top_at_1, & 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 @@ -1627,7 +1629,7 @@ subroutine lw_transport_1rescl(ncol, nlay, ngpt, top_at_1, & source_dn(icol,ilev,igpt) *trans(icol,ilev,igpt ) - & source_up(icol,ilev,igpt)) radn_up(icol,ilev,igpt) = radn_up(icol,ilev,igpt) + adjustmentFactor - enddo + 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) @@ -1640,7 +1642,7 @@ subroutine lw_transport_1rescl(ncol, nlay, ngpt, top_at_1, & 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 enddo else diff --git a/rte/kernels/mo_rte_solver_kernels.F90 b/rte/kernels/mo_rte_solver_kernels.F90 index 7c8604472..a88391536 100644 --- a/rte/kernels/mo_rte_solver_kernels.F90 +++ b/rte/kernels/mo_rte_solver_kernels.F90 @@ -423,9 +423,11 @@ subroutine lw_source_noscat(ncol, nlay, lay_source, lev_source_up, lev_source_dn ! is of order epsilon (smallest difference from 1. in working precision) ! Thanks to Peter Blossey ! - fact = merge((1._wp - trans(icol,ilay))/tau(icol,ilay) - trans(icol,ilay), & - tau(icol, ilay) * ( 0.5_wp - 1._wp/3._wp*tau(icol, ilay) ), & - tau(icol, ilay) > tau_thresh) + if(tau(icol, ilay) > tau_thresh) then + fact = (1._wp - trans(icol,ilay))/tau(icol,ilay) - trans(icol,ilay) + else + fact = tau(icol, ilay) * (0.5_wp - 1._wp/3._wp*tau(icol, ilay)) + end if ! ! Equation below is developed in Clough et al., 1992, doi:10.1029/92JD01419, Eq 13 ! @@ -1123,7 +1125,7 @@ subroutine lw_solver_1rescl_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_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] diff --git a/rte/mo_fluxes.F90 b/rte/mo_fluxes.F90 index f607b8a6b..48e2e8c45 100644 --- a/rte/mo_fluxes.F90 +++ b/rte/mo_fluxes.F90 @@ -18,6 +18,7 @@ ! ------------------------------------------------------------------------------------------------- module mo_fluxes use mo_rte_kind, only: wp + use mo_rte_config, only: check_extents use mo_rte_util_array, only: extents_are use mo_optical_props, only: ty_optical_props use mo_fluxes_broadband_kernels, & @@ -110,49 +111,43 @@ function reduce_broadband(this, gpt_flux_up, gpt_flux_dn, spectral_disc, top_at_ ngpt = size(gpt_flux_up, DIM=3) error_msg = "" - ! - ! Check array sizes - ! Input arrays - ! - if(.not. extents_are(gpt_flux_dn, ncol, nlev, ngpt)) then - error_msg = "reduce: gpt_flux_dn array incorrectly sized" - return - end if - if(present(gpt_flux_dn_dir)) then - if(.not. extents_are(gpt_flux_dn_dir, ncol, nlev, ngpt)) then - error_msg = "reduce: gpt_flux_dn_dir array incorrectly sized" - return + if(check_extents) then + ! + ! Check array sizes + ! Input arrays + ! + if(.not. extents_are(gpt_flux_dn, ncol, nlev, ngpt)) & + error_msg = "reduce: gpt_flux_dn array incorrectly sized" + + if(present(gpt_flux_dn_dir)) then + if(.not. extents_are(gpt_flux_dn_dir, ncol, nlev, ngpt)) & + error_msg = "reduce: gpt_flux_dn_dir array incorrectly sized" end if - end if - ! - ! Output arrays - ! - if(associated(this%flux_up)) then - if(.not. extents_are(this%flux_up, ncol, nlev)) then - error_msg = 'reduce: flux_up array incorrectly sized' - return + ! + ! Output arrays + ! + if(associated(this%flux_up)) then + if(.not. extents_are(this%flux_up, ncol, nlev)) & + error_msg = 'reduce: flux_up array incorrectly sized' end if - end if - if(associated(this%flux_dn)) then - if(.not. extents_are(this%flux_dn, ncol, nlev)) then - error_msg = 'reduce: flux_dn array incorrectly sized' - return + if(associated(this%flux_dn)) then + if(.not. extents_are(this%flux_dn, ncol, nlev)) & + error_msg = 'reduce: flux_dn array incorrectly sized' end if - end if - if(associated(this%flux_net)) then - if(.not. extents_are(this%flux_net, ncol, nlev)) then - error_msg = 'reduce: flux_net array incorrectly sized' - return + if(associated(this%flux_net)) then + if(.not. extents_are(this%flux_net, ncol, nlev)) & + error_msg = 'reduce: flux_net array incorrectly sized' end if - end if - if(associated(this%flux_dn_dir)) then - if(.not. extents_are(this%flux_dn_dir, ncol, nlev)) then - error_msg = 'reduce: flux_dn_dir array incorrectly sized' - return + if(associated(this%flux_dn_dir)) then + if(.not. extents_are(this%flux_dn_dir, ncol, nlev)) & + error_msg = 'reduce: flux_dn_dir array incorrectly sized' end if + + if(error_msg /= "") return end if ! ! Self-consistency -- shouldn't be asking for direct beam flux if it isn't supplied + ! if(associated(this%flux_dn_dir) .and. .not. present(gpt_flux_dn_dir)) then error_msg = "reduce: requesting direct downward flux but this hasn't been supplied" return diff --git a/rte/mo_optical_props.F90 b/rte/mo_optical_props.F90 index ffc8a3d7a..3aa2bf548 100644 --- a/rte/mo_optical_props.F90 +++ b/rte/mo_optical_props.F90 @@ -158,7 +158,7 @@ end function subset_range_abstract ! phase function moments (index 1 = g) for use with discrete ordinate methods ! ! ------------------------------------------------------------------------------------------------- - type, extends(ty_optical_props_arry) :: ty_optical_props_1scl + type, public, extends(ty_optical_props_arry) :: ty_optical_props_1scl contains procedure, public :: validate => validate_1scalar procedure, public :: get_subset => subset_1scl_range @@ -171,7 +171,7 @@ end function subset_range_abstract end type ! --- 2 stream ------------------------------------------------------------------------ - type, extends(ty_optical_props_arry) :: ty_optical_props_2str + type, public, extends(ty_optical_props_arry) :: ty_optical_props_2str real(wp), dimension(:,:,:), allocatable :: ssa ! single-scattering albedo (ncol, nlay, ngpt) real(wp), dimension(:,:,:), allocatable :: g ! asymmetry parameter (ncol, nlay, ngpt) contains @@ -186,7 +186,7 @@ end function subset_range_abstract end type ! --- n stream ------------------------------------------------------------------------ - type, extends(ty_optical_props_arry) :: ty_optical_props_nstr + type, public, extends(ty_optical_props_arry) :: ty_optical_props_nstr real(wp), dimension(:,:,:), allocatable :: ssa ! single-scattering albedo (ncol, nlay, ngpt) real(wp), dimension(:,:,:,:), allocatable :: p ! phase-function moments (nmom, ncol, nlay, ngpt) contains @@ -519,7 +519,7 @@ function delta_scale_2str(this, for) result(err_message) err_message = "delta_scale: dimension of 'for' don't match optical properties arrays" return end if - if(any(for < 0._wp .or. for > 1._wp)) then + if(any_vals_outside(for, 0._wp, 1._wp)) then err_message = "delta_scale: values of 'for' out of bounds [0,1]" return end if @@ -817,7 +817,7 @@ function increment(op_in, op_io) result(err_message) 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_config.F90 b/rte/mo_rte_config.F90 new file mode 100644 index 000000000..30dd2305e --- /dev/null +++ b/rte/mo_rte_config.F90 @@ -0,0 +1,45 @@ +! This code is part of RRTM for GCM Applications - Parallel (RRTMGP) +! +! Contacts: Robert Pincus and Eli Mlawer +! email: rrtmgp@aer.com +! +! Copyright 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 +! ------------------------------------------------------------------------------------------------- +! +! Control over input sanitization in Fortan front-end +! Module variables can be changed only by calling one of the included subroutine +! +! ------------------------------------------------------------------------------------------------- +module mo_rte_config + use mo_rte_kind, only: wl + implicit none + private + + logical(wl), protected, public :: check_extents = .true. + logical(wl), protected, public :: check_values = .true. + + interface rte_config_checks + module procedure rte_config_checks_each, rte_config_checks_all + end interface + public :: rte_config_checks +contains + ! -------------------------------------------------------------- + subroutine rte_config_checks_each(extents, values) + logical(wl), intent(in) :: extents, values + + check_extents = extents + check_values = values + end subroutine rte_config_checks_each + ! -------------------------------------------------------------- + subroutine rte_config_checks_all(do_checks) + logical(wl), intent(in) :: do_checks + + check_extents = do_checks + check_values = do_checks + end subroutine rte_config_checks_all + ! -------------------------------------------------------------- +end module mo_rte_config diff --git a/rte/mo_rte_kind.F90 b/rte/mo_rte_kind.F90 index daa09262f..de9e13e90 100644 --- a/rte/mo_rte_kind.F90 +++ b/rte/mo_rte_kind.F90 @@ -21,6 +21,7 @@ module mo_rte_kind use, intrinsic :: iso_c_binding, only: c_float, c_double, c_long, c_int, c_bool implicit none + public integer, parameter :: dp = c_double, sp = c_float, i8 = c_long, i4 = c_int ! ! Floating point working precision diff --git a/rte/mo_rte_lw.F90 b/rte/mo_rte_lw.F90 index 7ff091b2a..488f5ec34 100644 --- a/rte/mo_rte_lw.F90 +++ b/rte/mo_rte_lw.F90 @@ -35,6 +35,7 @@ ! ------------------------------------------------------------------------------------------------- module mo_rte_lw use mo_rte_kind, only: wp, wl + use mo_rte_config, only: check_extents, check_values use mo_rte_util_array,only: any_vals_less_than, any_vals_outside, extents_are use mo_optical_props, only: ty_optical_props, & ty_optical_props_arry, ty_optical_props_1scl, ty_optical_props_2str, ty_optical_props_nstr @@ -77,9 +78,9 @@ function rte_lw(optical_props, top_at_1, & 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) + target, optional, intent(inout) :: flux_up_Jac ! surface temperature flux Jacobian [W/m2/K] (ncol, nlay+1) real(wp), dimension(:,:), & - target, optional, intent(inout) :: flux_dn_Jac ! surface temperature flux Jacobian [W/m2/K] (ncol, ngpts) + target, optional, intent(inout) :: flux_dn_Jac ! surface temperature flux Jacobian [W/m2/K] (ncol, nlay+1) character(len=128) :: error_msg ! If empty, calculation was successful ! -------------------------------- ! @@ -133,7 +134,8 @@ function rte_lw(optical_props, top_at_1, & if(.not. fluxes%are_desired()) & error_msg = "rte_lw: no space allocated for fluxes" - if (present(flux_up_Jac)) then + + if (present(flux_up_Jac) .and. check_extents) then if( .not. extents_are(flux_up_Jac, ncol, nlay+1)) & error_msg = "rte_lw: flux Jacobian inconsistently sized" endif @@ -141,27 +143,40 @@ function rte_lw(optical_props, top_at_1, & ! ! Source functions ! - if(any([sources%get_ncol(), sources%get_nlay(), sources%get_ngpt()] /= [ncol, nlay, ngpt])) & - error_msg = "rte_lw: sources and optical properties inconsistently sized" + if (check_extents) then + if(any([sources%get_ncol(), sources%get_nlay(), sources%get_ngpt()] /= [ncol, nlay, ngpt])) & + error_msg = "rte_lw: sources and optical properties inconsistently sized" + end if ! Also need to validate - ! - ! Surface emissivity - ! - if(.not. extents_are(sfc_emis, nband, ncol)) & - error_msg = "rte_lw: sfc_emis inconsistently sized" - if(any_vals_outside(sfc_emis, 0._wp, 1._wp)) & - error_msg = "rte_lw: sfc_emis has values < 0 or > 1" - if(len_trim(error_msg) > 0) return + if (check_extents) then + ! + ! Surface emissivity + ! + if(.not. extents_are(sfc_emis, nband, ncol)) & + error_msg = "rte_lw: sfc_emis inconsistently sized" + ! + ! Incident flux, if present + ! + if(present(inc_flux)) then + if(.not. extents_are(inc_flux, ncol, ngpt)) & + error_msg = "rte_lw: inc_flux inconsistently sized" + end if + end if - ! - ! Incident flux, if present - ! - if(present(inc_flux)) then - if(.not. extents_are(inc_flux, ncol, ngpt)) & - error_msg = "rte_lw: inc_flux inconsistently sized" - if(any_vals_less_than(inc_flux, 0._wp)) & - error_msg = "rte_lw: inc_flux has values < 0" + if(check_values) then + if(any_vals_outside(sfc_emis, 0._wp, 1._wp)) & + error_msg = "rte_lw: sfc_emis has values < 0 or > 1" + if(present(inc_flux)) then + if(any_vals_less_than(inc_flux, 0._wp)) & + error_msg = "rte_lw: inc_flux has values < 0" + end if + 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" + if(n_gauss_angles < 1) & + error_msg = "rte_lw: have to ask for at least one quadrature point for no-scattering calculation" + end if end if if(len_trim(error_msg) > 0) return @@ -169,15 +184,7 @@ function rte_lw(optical_props, top_at_1, & ! Number of quadrature points for no-scattering calculation ! 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 RT calculation" - if(n_gauss_angles < 1) & - 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 - + if(present(n_gauss_angles)) n_quad_angs = n_gauss_angles ! ! Optionally - use 2-stream methods when low-order scattering properties are provided? ! @@ -207,14 +214,15 @@ function rte_lw(optical_props, top_at_1, & 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") + error_msg = "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 if using scattering ! - error_msg = optical_props%validate() + if(check_values) error_msg = optical_props%validate() + if(len_trim(error_msg) > 0) then if(len_trim(optical_props%get_name()) > 0) & error_msg = trim(optical_props%get_name()) // ': ' // trim(error_msg) @@ -284,7 +292,8 @@ function rte_lw(optical_props, top_at_1, & end if !$acc exit data delete(optical_props%tau) class is (ty_optical_props_2str) - if ((present(flux_dn_Jac))) then + + if (present(flux_dn_Jac) .and. check_extents) then if( .not. extents_are(flux_dn_Jac, ncol, nlay+1)) then error_msg = "rte_lw: flux_dn_Jac inconsistently sized" return diff --git a/rte/mo_rte_sw.F90 b/rte/mo_rte_sw.F90 index aece7971d..18169b3e3 100644 --- a/rte/mo_rte_sw.F90 +++ b/rte/mo_rte_sw.F90 @@ -29,6 +29,7 @@ ! ------------------------------------------------------------------------------------------------- module mo_rte_sw use mo_rte_kind, only: wp, wl + use mo_rte_config, only: check_extents, check_values use mo_rte_util_array,only: any_vals_less_than, any_vals_outside, extents_are use mo_optical_props, only: ty_optical_props, & ty_optical_props_arry, ty_optical_props_1scl, ty_optical_props_2str, ty_optical_props_nstr @@ -84,32 +85,41 @@ function rte_sw(atmos, top_at_1, & end if ! - ! Sizes and values of input arrays + ! Sizes of input arrays ! - if(.not. extents_are(mu0, ncol)) & - error_msg = "rte_sw: mu0 inconsistently sized" - if(any_vals_outside(mu0, 0._wp, 1._wp)) & - error_msg = "rte_sw: one or more mu0 <= 0 or > 1" + if(check_extents) then + if(.not. extents_are(mu0, ncol)) & + error_msg = "rte_sw: mu0 inconsistently sized" + if(.not. extents_are(inc_flux, ncol, ngpt)) & + error_msg = "rte_sw: inc_flux inconsistently sized" + if(.not. extents_are(sfc_alb_dir, nband, ncol)) & + error_msg = "rte_sw: sfc_alb_dir inconsistently sized" + if(.not. extents_are(sfc_alb_dif, nband, ncol)) & + error_msg = "rte_sw: sfc_alb_dif inconsistently sized" + if(present(inc_flux_dif)) then + if(.not. extents_are(inc_flux_dif, ncol, ngpt)) & + error_msg = "rte_sw: inc_flux_dif inconsistently sized" + end if + end if - if(.not. extents_are(inc_flux, ncol, ngpt)) & - error_msg = "rte_sw: inc_flux inconsistently sized" - if(any_vals_less_than(inc_flux, 0._wp)) & - error_msg = "rte_sw: one or more inc_flux < 0" - if(present(inc_flux_dif)) then - if(.not. extents_are(inc_flux_dif, ncol, ngpt)) & - error_msg = "rte_sw: inc_flux_dif inconsistently sized" - if(any_vals_less_than(inc_flux_dif, 0._wp)) & - error_msg = "rte_sw: one or more inc_flux_dif < 0" + ! + ! Values of input arrays + ! + if(check_values) then + if(any_vals_outside(mu0, 0._wp, 1._wp)) & + error_msg = "rte_sw: one or more mu0 <= 0 or > 1" + if(any_vals_less_than(inc_flux, 0._wp)) & + error_msg = "rte_sw: one or more inc_flux < 0" + if(any_vals_outside(sfc_alb_dir, 0._wp, 1._wp)) & + error_msg = "rte_sw: sfc_alb_dir out of bounds [0,1]" + if(any_vals_outside(sfc_alb_dif, 0._wp, 1._wp)) & + error_msg = "rte_sw: sfc_alb_dif out of bounds [0,1]" + if(present(inc_flux_dif)) then + if(any_vals_less_than(inc_flux_dif, 0._wp)) & + error_msg = "rte_sw: one or more inc_flux_dif < 0" + end if end if - if(.not. extents_are(sfc_alb_dir, nband, ncol)) & - error_msg = "rte_sw: sfc_alb_dir inconsistently sized" - if(any_vals_outside(sfc_alb_dir, 0._wp, 1._wp)) & - error_msg = "rte_sw: sfc_alb_dir out of bounds [0,1]" - if(.not. extents_are(sfc_alb_dif, nband, ncol)) & - error_msg = "rte_sw: sfc_alb_dif inconsistently sized" - if(any_vals_outside(sfc_alb_dif, 0._wp, 1._wp)) & - error_msg = "rte_sw: sfc_alb_dif out of bounds [0,1]" if(len_trim(error_msg) > 0) then if(len_trim(atmos%get_name()) > 0) &