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) &