diff --git a/.cmake-format.py b/.cmake-format.py new file mode 100644 index 000000000..f6ee981d9 --- /dev/null +++ b/.cmake-format.py @@ -0,0 +1,25 @@ +with section("parse"): + additional_commands = {"check_python3_package": {"pargs": 1, "kwargs": {"CODE": 1}}} + +with section("format"): + dangle_parens = True + max_lines_hwrap = 0 + keyword_case = "upper" + autosort = True + +with section("lint"): + # The formatter sometimes fails to fit the code into the line limit (C0301) and can + # disagree with the linter regarding the indentation (C0307): + disabled_codes = ["C0301", "C0307"] + # Names of local variables must be in lowercase but sometimes we need to + # override standard CMake variables: + local_var_pattern = "CMAKE_[0-9A-Z_]+|[a-z][0-9a-z_]+" + # The standard names of the languages in CMake are C and Fortran. Names of + # private variables must be in lowercase but can have substings "C" and + # "Fortran": + private_var_pattern = ( + "([a-z_][0-9a-z_]*_)?(C|Fortran)(_[a-z_][0-9a-z_]*)?|[a-z_][0-9a-z_]+" + ) + # The standard name of the language in CMake is Fortran. Names of public + # variables must be in uppercase but can have substring "Fortran": + public_var_pattern = "([A-Z][0-9A-Z_]*_)?Fortran(_[A-Z][0-9A-Z_]*)?|[A-Z][0-9A-Z_]+" diff --git a/.git-blame-ignore-revs b/.git-blame-ignore-revs new file mode 100644 index 000000000..593947014 --- /dev/null +++ b/.git-blame-ignore-revs @@ -0,0 +1,2 @@ +# Formatted entire CMake code base with cmake-format +45b43632309cff022326472a8be0fdd5efc8f5c8 diff --git a/.github/workflows/check-api.yml b/.github/workflows/check-api.yml index b3681e666..4b50a0355 100644 --- a/.github/workflows/check-api.yml +++ b/.github/workflows/check-api.yml @@ -14,11 +14,11 @@ jobs: API: runs-on: ubuntu-22.04 env: - # Core variables: FC: gfortran-12 - FCFLAGS: "-ffree-line-length-none -m64 -std=f2008 -march=native -fbounds-check -fmodule-private -fimplicit-none -finit-real=nan -g -DRTE_USE_CBOOL" - RRTMGP_ROOT: ${{ github.workspace }} + FFLAGS: "-m64 -std=f2008 -march=native -fbounds-check -fmodule-private -fimplicit-none -finit-real=nan" RTE_KERNELS: extern + CMAKE_BUILD_PARALLEL_LEVEL: 8 + VERBOSE: steps: # # Check out repository under $GITHUB_WORKSPACE @@ -31,4 +31,7 @@ jobs: - name: Build libraries run: | $FC --version - make -j4 libs \ No newline at end of file + cmake -S . -B build \ + -DCMAKE_BUILD_TYPE=RelWithDebInfo \ + -DKERNEL_MODE=$RTE_KERNELS + cmake --build build diff --git a/.github/workflows/containerized-ci.yml b/.github/workflows/containerized-ci.yml index 3a931c60b..cf6ffad2d 100644 --- a/.github/workflows/containerized-ci.yml +++ b/.github/workflows/containerized-ci.yml @@ -18,24 +18,31 @@ jobs: fortran-compiler: [ifort, ifx, nvfortran] rte-kernels: [default, accel] fpmodel: [DP, SP] + exclude: + # Fails with error #5633: **Internal compiler error: segmentation violation signal raised** + - fortran-compiler: ifx + rte-kernels: accel + # fcflags: -debug -traceback -O0 -heap-arrays -assume realloc_lhs -extend-source 132 -stand f08 -fiopenmp -fopenmp-targets=spir64 + # build-type: None include: # Set flags for Intel Fortran Compiler Classic - fortran-compiler: ifort - fcflags: -m64 -g -traceback -heap-arrays -assume realloc_lhs -extend-source 132 -check bounds,uninit,pointers,stack -stand f08 -diag-disable=10448 + fcflags: -m64 -traceback -heap-arrays -assume realloc_lhs -extend-source 132 -check bounds,uninit,pointers,stack -stand f08 -diag-disable=10448 + build-type: RelWithDebInfo # Set flags for Intel Fortran Compiler - fortran-compiler: ifx rte-kernels: default fcflags: -debug -traceback -O0 -heap-arrays -assume realloc_lhs -extend-source 132 -stand f08 - - fortran-compiler: ifx - rte-kernels: accel - fcflags: -debug -traceback -O0 -heap-arrays -assume realloc_lhs -extend-source 132 -stand f08 -fiopenmp -fopenmp-targets=spir64 + build-type: None # Set flags for NVIDIA Fortran compiler - fortran-compiler: nvfortran rte-kernels: default fcflags: -Mallocatable=03 -Mstandard -Mbounds -Mchkptr -Kieee -Mchkstk + build-type: None - fortran-compiler: nvfortran rte-kernels: accel fcflags: -Mallocatable=03 -Mstandard -Mbounds -Mchkptr -Kieee -Mchkstk -acc + build-type: None # Set container images - fortran-compiler: ifort image: ghcr.io/earth-system-radiation/rte-rrtmgp-ci:oneapi @@ -46,86 +53,53 @@ jobs: container: image: ${{ matrix.image }} env: - # Core variables: FC: ${{ matrix.fortran-compiler }} - FCFLAGS: ${{ matrix.fcflags }} -DRTE_USE_${{ matrix.fpmodel}} - # Make variables: - NFHOME: /opt/netcdf-fortran - RRTMGP_ROOT: ${{ github.workspace }} - RRTMGP_DATA: ${{ github.workspace }}/rrtmgp-data - RTE_KERNELS: ${{ matrix.rte-kernels }} - RUN_CMD: + FFLAGS: ${{ matrix.fcflags }} + NetCDF_Fortran_ROOT: /opt/netcdf-fortran + CMAKE_BUILD_PARALLEL_LEVEL: 8 + VERBOSE: + CTEST_PARALLEL_LEVEL: 8 + CTEST_OUTPUT_ON_FAILURE: 1 # https://github.com/earth-system-radiation/rte-rrtmgp/issues/194 OMP_TARGET_OFFLOAD: DISABLED - FAILURE_THRESHOLD: 7.e-4 steps: # - # Checks-out repository under $GITHUB_WORKSPACE - # - - uses: actions/checkout@v4 - # - # Check out data + # Check out repository under $GITHUB_WORKSPACE # - - name: Check out data + - name: Check out code uses: actions/checkout@v4 - with: - repository: earth-system-radiation/rrtmgp-data - path: rrtmgp-data - ref: v1.8.1 # - # Build libraries, examples and tests (expect success) + # Install required tools # - - name: Build libraries, examples and tests (expect success) - id: build-success - if: matrix.fortran-compiler != 'ifx' || matrix.rte-kernels != 'accel' + - name: Install Required Tools run: | - $FC --version - make -j4 libs + apt-get update + apt-get install -y git cmake ninja-build # - # Build libraries, examples and tests (expect failure) + # Build libraries, examples and tests # - - name: Build libraries, examples and tests (expect failure) - if: steps.build-success.outcome == 'skipped' - shell: bash + - name: Build libraries, examples and tests run: | - $FC --version - make -j4 libs 2> >(tee make.err >&2) && { - echo "Unexpected success" - exit 1 - } || { - grep make.err -e 'Internal compiler error' && { - echo "Expected failure" - } || { - echo "Unexpected failure" - exit 1 - } - } + cmake -S . -B build -G "Ninja" \ + -DCMAKE_BUILD_TYPE=${{ matrix.build-type }} \ + -DRTE_ENABLE_SP=`test 'x${{ matrix.fpmodel }}' = xSP && echo ON || echo OFF` \ + -DKERNEL_MODE=${{ matrix.rte-kernels }} \ + -DBUILD_TESTING=ON + cmake --build build # # Run examples and tests # - name: Run examples and tests - if: steps.build-success.outcome != 'skipped' - run: make -j4 tests - # - # Relax failure thresholds for single precision - # - - name: Relax failure threshold for single precision - if: matrix.fpmodel == 'SP' && steps.build-success.outcome != 'skipped' - run: echo "FAILURE_THRESHOLD=3.5e-1" >> $GITHUB_ENV - # - # Compare the results - # - - name: Compare the results - if: steps.build-success.outcome != 'skipped' - run: make -j4 check + working-directory: build + run: ctest # # Generate validation plots # - name: Generate validation plots if: matrix.fortran-compiler == 'ifort' && matrix.rte-kernels == 'default' && matrix.fpmodel == 'DP' - working-directory: tests - run: python validation-plots.py + run: | + cmake --build build --target validation-plots # # Upload validation plots # @@ -134,4 +108,4 @@ jobs: uses: actions/upload-artifact@v4 with: name: valdiation-plot - path: tests/validation-figures.pdf + path: build/tests/validation-figures.pdf diff --git a/.github/workflows/continuous-integration.yml b/.github/workflows/continuous-integration.yml index d8a52f32d..d9e346e9e 100644 --- a/.github/workflows/continuous-integration.yml +++ b/.github/workflows/continuous-integration.yml @@ -1,12 +1,13 @@ name: Continuous Integration + on: push: branches: - - main - - develop + - main + - develop pull_request: branches-ignore: - - documentation + - documentation workflow_dispatch: defaults: @@ -16,88 +17,115 @@ defaults: jobs: CI: - runs-on: ubuntu-22.04 strategy: fail-fast: false matrix: - fortran-compiler: [gfortran-10, gfortran-11, gfortran-12] + os: [ubuntu-24.04, macos-13, macos-latest, windows-2022] + gfortran-version: [12, 13, 14] + gfortran-from: [system, conda] fpmodel: [DP, SP] + exclude: + - os: ubuntu-24.04 + gfortran-from: conda + - os: macos-13 + gfortran-from: system + - os: macos-13 + gfortran-version: 14 + gfortran-from: conda + - os: macos-latest + gfortran-from: system + - os: macos-latest + gfortran-version: 14 + gfortran-from: conda + - os: windows-2022 + include: + - os: ubuntu-24.04 + gfortran-version: 13 + gfortran-from: conda + fpmodel: DP + - os: windows-2022 + gfortran-version: 13 + gfortran-from: conda + fpmodel: DP + - os: windows-2022 + gfortran-version: 13 + gfortran-from: conda + fpmodel: SP + - os: windows-2022 + gfortran-version: 14 + gfortran-from: conda + fpmodel: DP + - os: windows-2022 + gfortran-version: 14 + gfortran-from: conda + fpmodel: SP env: - # Core variables: - FC: ${{ matrix.fortran-compiler }} - FCFLAGS: "-ffree-line-length-none -m64 -std=f2008 -march=native -fbounds-check -fmodule-private -fimplicit-none -finit-real=nan -g -DRTE_USE_CBOOL -DRTE_USE_${{ matrix.fpmodel }}" - # Make variables: - FCINCLUDE: -I/usr/include - RRTMGP_ROOT: ${{ github.workspace }} - RRTMGP_DATA: ${{ github.workspace }}/rrtmgp-data - RUN_CMD: - FAILURE_THRESHOLD: 7.e-4 + FC: gfortran-${{ matrix.gfortran-version }} + FFLAGS: "-m64 -std=f2008 -march=native -fbounds-check -fmodule-private -fimplicit-none -finit-real=nan" + CMAKE_BUILD_PARALLEL_LEVEL: 8 + VERBOSE: + CTEST_PARALLEL_LEVEL: 8 + CTEST_OUTPUT_ON_FAILURE: 1 + runs-on: ${{ matrix.os }} steps: - # - # Relax failure thresholds for single precision - # - - name: Relax failure threshold for single precision - if: matrix.fpmodel == 'SP' - run: echo "FAILURE_THRESHOLD=3.5e-1" >> $GITHUB_ENV - # - # Check out repository under $GITHUB_WORKSPACE - # - - name: Check out code - uses: actions/checkout@v4 - # - # Check out data - # - - name: Check out data - uses: actions/checkout@v4 - with: - repository: earth-system-radiation/rrtmgp-data - path: rrtmgp-data - ref: v1.8.1 - # - # Synchronize the package index - # - - name: Synchronize the package index - run: sudo apt-get update - # - # Install NetCDF-Fortran (compatible with all compilers) - # - - name: Install NetCDF-Fortran - run: sudo apt-get install libnetcdff-dev - # - # Cache Conda packages - # - - name: Cache Conda packages - uses: actions/cache@v4 - with: - path: ~/conda_pkgs_dir - key: conda-pkgs - # - # Set up Conda - # - - name: Set up Conda - uses: conda-incubator/setup-miniconda@v3 - with: - miniforge-version: latest - activate-environment: rte_rrtmgp_test - environment-file: environment-noplots.yml - python-version: 3.11 - auto-activate-base: false - # Use the cache properly: - use-only-tar-bz2: true - # - # Build libraries, examples and tests - # - - name: Build libraries - run: | - $FC --version - make -j4 libs - # - # Run examples and tests - # - - name: Build and run examples and tests - run: make -j4 tests - # - # Compare the results - # - - name: Compare the results - run: make -j4 check + # + # Check out repository under $GITHUB_WORKSPACE + # + - name: Check out code + uses: actions/checkout@v4 + # + # Cache Conda packages + # + - name: Cache Conda packages + uses: actions/cache@v4 + with: + path: ~/conda_pkgs_dir + key: conda-pkgs-${{ matrix.os }} + # + # Set up Conda + # + - name: Set up Conda + uses: conda-incubator/setup-miniconda@v3 + with: + miniforge-version: latest + activate-environment: rte_rrtmgp_test + environment-file: environment-noplots.yml + python-version: 3.11 + auto-activate-base: false + # Use the cache properly: + use-only-tar-bz2: false + # + # Install compiler + # + - name: Install compiler + if: matrix.gfortran-from == 'conda' + run: | + conda install -c conda-forge gfortran=${{ matrix.gfortran-version }} -y + echo "FC=gfortran" >> $GITHUB_ENV + # + # Install dependencies + # + - name: Install dependencies + run: conda install -c conda-forge netcdf-fortran ninja -y + # + # Adjust toolchain + # + - name: Adjust toolchain + if: matrix.os == 'windows-2022' + run: echo "FC=${FC}.exe" >> $GITHUB_ENV + # + # Build libraries, examples, and tests + # + - name: Build libraries and tests + run: | + cmake -S . -B build -G "Ninja" \ + -DCMAKE_BUILD_TYPE=RelWithDebInfo \ + -DRTE_ENABLE_SP=`test 'x${{ matrix.fpmodel }}' = xSP && echo ON || echo OFF` \ + -DBUILD_TESTING=ON + cmake --build build + # + # Run examples, tests and checks + # + - name: Run examples, tests and checks + working-directory: build + run: ctest diff --git a/.github/workflows/doc-deployment.yml b/.github/workflows/doc-deployment.yml index 4bc492e9b..2c3a5c216 100644 --- a/.github/workflows/doc-deployment.yml +++ b/.github/workflows/doc-deployment.yml @@ -70,7 +70,7 @@ jobs: # Deploy documentation # - name: Deploy API Documentation - uses: JamesIves/github-pages-deploy-action@v4.6.1 + uses: JamesIves/github-pages-deploy-action@v4.7.2 if: ${{ github.event_name == 'push' && github.ref == 'refs/heads/documentation' }} with: branch: gh-pages diff --git a/.github/workflows/gitlab-ci.yml b/.github/workflows/gitlab-ci.yml index 2a46e8557..feacb5f7f 100644 --- a/.github/workflows/gitlab-ci.yml +++ b/.github/workflows/gitlab-ci.yml @@ -25,6 +25,7 @@ jobs: github.event.pull_request.user.login != 'dependabot[bot]' )) runs-on: ubuntu-latest outputs: + ref-type: ${{ steps.g-push-rev.outputs.ref-type }} ref-name: ${{ steps.g-push-rev.outputs.ref-name }} pipeline-id: ${{ steps.gl-trigger-pipeline.outputs.pipeline-id }} steps: @@ -94,7 +95,7 @@ jobs: with: remote-url: ${{ vars.DKRZ_GITLAB_SERVER }}/${{ vars.DKRZ_GITLAB_PROJECT }}.git password: ${{ secrets.DKRZ_GITLAB_TOKEN }} - ref-type: tag + ref-type: ${{ needs.levante-init.outputs.ref-type }} ref-name: ${{ needs.levante-init.outputs.ref-name }} force: true # @@ -102,12 +103,14 @@ jobs: # lumi-init: if: | + false && github.repository_owner == 'earth-system-radiation' && ( github.event_name != 'pull_request' || ( github.event.pull_request.head.repo.owner.login == github.repository_owner && github.event.pull_request.user.login != 'dependabot[bot]' )) runs-on: ubuntu-latest outputs: + ref-type: ${{ steps.g-push-rev.outputs.ref-type }} ref-name: ${{ steps.g-push-rev.outputs.ref-name }} pipeline-id: ${{ steps.gl-create-pipeline.outputs.pipeline-id }} steps: @@ -130,8 +133,7 @@ jobs: rev-id: ${{ github.sha }} rev-signing-format: ssh rev-signing-key: ${{ secrets.GITLAB_SIGNING_KEY }} - ref-type: tag - ref-message: ${{ github.server_url }}/${{ github.repository }}/actions/runs/${{ github.run_id }} + ref-type: branch force-push: true # # Create GitLab CI/CD Pipeline @@ -189,11 +191,12 @@ jobs: project-name: ${{ vars.GITLAB_PROJECT }} token: ${{ secrets.GITLAB_TOKEN }} pipeline-id: ${{ needs.lumi-init.outputs.pipeline-id }} + force: true - uses: "skosukhin/git-ci-hub-lab/gl-delete-ref@v1" with: server-url: ${{ vars.GITLAB_SERVER }} project-name: ${{ vars.GITLAB_PROJECT }} token: ${{ secrets.GITLAB_TOKEN }} - ref-type: tag + ref-type: ${{ needs.lumi-init.outputs.ref-type }} ref-name: ${{ needs.lumi-init.outputs.ref-name }} force: true diff --git a/.github/workflows/self-hosted-ci.yml b/.github/workflows/self-hosted-ci.yml deleted file mode 100644 index 7077eda6f..000000000 --- a/.github/workflows/self-hosted-ci.yml +++ /dev/null @@ -1,129 +0,0 @@ -name: Self-hosted CI -on: - push: - branches: - - main - - develop - pull_request: - branches-ignore: - - documentation - workflow_dispatch: - -defaults: - run: - shell: bash - -jobs: - CI: - if: github.repository == 'earth-system-radiation/rte-rrtmgp' - runs-on: - labels: cscs-ci - strategy: - fail-fast: false - matrix: - config-name: [nvidia-gpu-openacc, cce-cpu-icon-production, cce-gpu-openmp] - fpmodel: [DP, SP] - include: - - config-name: nvidia-gpu-openacc - rte-kernels: accel - compiler-modules: "PrgEnv-nvidia nvidia craype-accel-nvidia60 cdt-cuda/21.09 !cray-libsci_acc" - # Generic accelerator flag - fcflags: "-O3 -acc -Mallocatable=03 -gopt" - - config-name: cce-cpu-icon-production - rte-kernels: default - compiler-modules: "PrgEnv-cray" - # Production flags for Icon model - fcflags: "-hadd_paren -r am -Ktrap=divz,ovf,inv -hflex_mp=intolerant -hfp1 -hnoacc -O1,cache0" - - config-name: cce-gpu-openmp - rte-kernels: accel - compiler-modules: "PrgEnv-cray craype-accel-nvidia60 cdt-cuda/22.05 cudatoolkit/11.2.0_3.39-2.1__gf93aa1c" - # OpenMP flags from Nichols Romero (Argonne) - fcflags: "-hnoacc -homp -O0" - env: - # Core variables: - FC: ftn - FCFLAGS: ${{ matrix.fcflags }} -DRTE_USE_${{ matrix.fpmodel }} - # Make variables: - RRTMGP_ROOT: ${{ github.workspace }} - RRTMGP_DATA: ${{ github.workspace }}/rrtmgp-data - RTE_KERNELS: ${{ matrix.rte-kernels }} - RUN_CMD: "srun -C gpu -A d56 -p cscsci -t 15:00" - FAILURE_THRESHOLD: 7.e-4 - steps: - # - # Checks-out repository under $GITHUB_WORKSPACE - # - - uses: actions/checkout@v3 - # - # Check out data - # - - name: Check out data - uses: actions/checkout@v3 - with: - repository: earth-system-radiation/rrtmgp-data - path: rrtmgp-data - ref: v1.8.1 - # - # Finalize build environment - # - - name: Finalize build environment - run: | - # There are significant limitations on what can go into ${GITHUB_ENV}, - # therefore, we use ${BASH_ENV} but only when necessary: - BASH_ENV="${GITHUB_WORKSPACE}/.bash" - echo "source '${GITHUB_WORKSPACE}/.github/workflows/module_switcher'" >> "${BASH_ENV}" - echo "switch_for_module daint-gpu ${{ matrix.compiler-modules }} cray-netcdf cray-hdf5" >> "${BASH_ENV}" - # Use custom Python environment: - # The environment can be re-generated as follows: - # module load cray-python - # python3 -m venv /scratch/snx3000/rpincus/rte-rrtmgp-python - # /scratch/snx3000/rpincus/rte-rrtmgp-python/bin/pip3 install --upgrade pip - # /scratch/snx3000/rpincus/rte-rrtmgp-python/bin/pip3 install dask[array] netCDF4 numpy xarray - echo 'PATH="/scratch/snx3000/rpincus/rte-rrtmgp-python/bin:${PATH}"' >> "${BASH_ENV}" - # Make bash run the above on startup: - echo "BASH_ENV=${BASH_ENV}" >> "${GITHUB_ENV}" - # Compiler needs more temporary space than normally available: - tmpdir='${{ github.workspace }}/tmp' - mkdir "${tmpdir}" && echo "TMPDIR=${tmpdir}" >> "${GITHUB_ENV}" - # We use the "non-default products" for the tests - # (see https://support.hpe.com/hpesc/public/docDisplay?docId=a00113984en_us&page=Modify_Linking_Behavior_to_Use_Non-default_Libraries.html): - echo 'LD_LIBRARY_PATH="${CRAY_LD_LIBRARY_PATH}:${LD_LIBRARY_PATH}"' >> "${BASH_ENV}" - # SLURM jobs, user home directories and HDF5 file locking are - # incompatible on Daint: - echo 'HDF5_USE_FILE_LOCKING=FALSE' >> "${GITHUB_ENV}" - # - # Build libraries, examples and tests - # - - name: Build libraries - run: | - $FC --version - make -j8 libs - # - # Run examples and tests (expect success) - # - - name: Build and run examples and tests (expect success) - id: run-success - if: matrix.config-name != 'cce-gpu-openmp' - run: make -j8 tests - # - # Run examples and tests (expect failure) - # - - name: Build and run examples and tests (expect failure) - if: steps.run-success.outcome == 'skipped' - run: | - make -j8 tests && { - echo "Unexpected success" - exit 1 - } || echo "Expected failure" - # - # Relax failure thresholds for single precision - # - - name: Relax failure threshold for single precision - if: matrix.fpmodel == 'SP' && steps.run-success.outcome != 'skipped' - run: echo "FAILURE_THRESHOLD=3.5e-1" >> $GITHUB_ENV - # - # Compare the results - # - - name: Compare the results - if: steps.run-success.outcome != 'skipped' - run: make -j8 check diff --git a/.github/workflows/style.yml b/.github/workflows/style.yml new file mode 100644 index 000000000..a156e744e --- /dev/null +++ b/.github/workflows/style.yml @@ -0,0 +1,72 @@ +name: Check Style +on: + push: + branches: + - main + - develop + pull_request: + branches-ignore: + - documentation + workflow_dispatch: + +env: + FIND_CMAKE_FILES_CMD: "find '${{ github.workspace }}' -name 'CMakeLists.txt' -o -name '*.cmake' -o -name '*.cmake.in'" + +jobs: + Format: + runs-on: ubuntu-22.04 + env: + DEFAULT: '\033[0m' + RED: '\033[0;31m' + GREEN: '\033[0;32m' + FORMAT_PATCH: '${{ github.workspace }}/format.patch' + steps: + - name: Check out code + uses: actions/checkout@v4 + - name: Install Python + uses: actions/setup-python@v5 + with: + python-version: '>=3.8' + - name: Install required tools + run: python -m pip install cmake-format + - name: Format CMake + run: cmake-format -i $(eval "${FIND_CMAKE_FILES_CMD}") + - name: Check if patching is required + id: patch-required + run: | + git -C '${{ github.workspace }}' diff --patch-with-raw > "${FORMAT_PATCH}" + test -s "${FORMAT_PATCH}" && { + printf "${RED}The source code does not meet the format requirements. \ + Please, apply the patch (see artifacts).${DEFAULT}\n" + + printf "${RED}Note that the result of the formatting might depend \ + on the versions of the formatting tools. In this project, whatever \ + formatting this CI job produces if the correct one. If it expects \ + you to reformat parts of the source code that you did not modify, do \ + so in a separate commit, which must not be squashed, and list the \ + commit in the '.git-blame-ignore-revs' file.${DEFAULT}\n" + + exit 1 + } || { + printf "${GREEN}The source code meets the format requirements.${DEFAULT}\n" + rm -rf "${FORMAT_PATCH}" + } + - name: Upload the patch file + if: always() && steps.patch-required.outcome == 'failure' + uses: actions/upload-artifact@v4 + with: + name: format-patch + path: ${{ env.FORMAT_PATCH }} + Lint: + runs-on: ubuntu-22.04 + steps: + - name: Check out code + uses: actions/checkout@v4 + - name: Install Python + uses: actions/setup-python@v5 + with: + python-version: '>=3.8' + - name: Install required tools + run: python -m pip install cmake-format + - name: Lint CMake + run: cmake-lint $(eval "${FIND_CMAKE_FILES_CMD}") diff --git a/.gitignore b/.gitignore index 2ae1d6a51..e32360bf6 100644 --- a/.gitignore +++ b/.gitignore @@ -33,8 +33,14 @@ examples/*/*/*.nc *.html *.gif +# Python virtual environments +.venv*/ + # gh-pages directory public +# build +build/ + # Ruby stuff **/Gemfile.lock diff --git a/.gitlab/common.yml b/.gitlab/common.yml index 83a267840..81aac69e4 100644 --- a/.gitlab/common.yml +++ b/.gitlab/common.yml @@ -1,37 +1,29 @@ .dp: variables: - FPMODEL: DP - FAILURE_THRESHOLD: "7.e-4" + RTE_ENABLE_SP: OFF .sp: variables: - FPMODEL: SP - FAILURE_THRESHOLD: "3.5e-1" + RTE_ENABLE_SP: ON .common: variables: - # Make variables: - MAKEFLAGS: -j8 - RRTMGP_ROOT: ${CI_PROJECT_DIR} - RRTMGP_DATA: ${CI_PROJECT_DIR}/rrtmgp-data - # Convenience variables: - RRTMGP_DATA_REPO: https://github.com/earth-system-radiation/rrtmgp-data.git - RRTMGP_DATA_TAG: v1.8.1 + CMAKE_BUILD_PARALLEL_LEVEL: 8 + VERBOSE: + CTEST_PARALLEL_LEVEL: 8 + CTEST_OUTPUT_ON_FAILURE: 1 script: # # Build libraries, examples and tests # - - ${FC} ${VERSION_FCFLAGS} - - make libs + - | + cmake -S . -B build \ + -DCMAKE_BUILD_TYPE=None \ + -DRTE_ENABLE_SP=$RTE_ENABLE_SP \ + -DKERNEL_MODE=$KERNEL_MODE \ + -DBUILD_TESTING=ON + - cmake --build build # - # Check out data + # Run examples, tests and checks # - - git clone --depth 1 ${RRTMGP_DATA_TAG:+--branch "${RRTMGP_DATA_TAG}"} "${RRTMGP_DATA_REPO}" "${RRTMGP_DATA}" - # - # Run examples and tests - # - - make tests - # - # Compare the results - # - - make check + - ctest --test-dir build diff --git a/.gitlab/levante.yml b/.gitlab/levante.yml index be6f253b2..677cf4039 100644 --- a/.gitlab/levante.yml +++ b/.gitlab/levante.yml @@ -29,21 +29,13 @@ variables: .nvhpc: variables: - # Core variables: - FC: /sw/spack-levante/nvhpc-22.5-v4oky3/Linux_x86_64/22.5/compilers/bin/nvfortran - # Convenience variables: - VERSION_FCFLAGS: --version - NFHOME: /sw/spack-levante/netcdf-fortran-4.5.4-syv4qr - NCHOME: /sw/spack-levante/netcdf-c-4.9.0-gc7kgj + FC: /sw/spack-levante/nvhpc-24.9-p7iohv/Linux_x86_64/24.9/compilers/bin/nvfortran + NetCDF_Fortran_ROOT: /sw/spack-levante/netcdf-fortran-4.6.1-4wu5wt .nag: variables: - # Core variables: FC: /sw/spack-levante/nag-7.1-lqjbej/bin/nagfor - # Convenience variables: - VERSION_FCFLAGS: -V - NFHOME: /sw/spack-levante/netcdf-fortran-4.5.3-5di6qe - NCHOME: /sw/spack-levante/netcdf-c-4.8.1-vbnli5 + NetCDF_Fortran_ROOT: /sw/spack-levante/netcdf-fortran-4.5.3-5di6qe .common-levante: extends: .common @@ -51,15 +43,11 @@ variables: PYHOME: /sw/spack-levante/mambaforge-22.9.0-2-Linux-x86_64-kptncg # Suppress an irrelevant but annoying error message: PROJ_LIB: ${PYHOME}/share/proj - # Make variables: - FCINCLUDE: -I${NFHOME}/include - LDFLAGS: -L${NFHOME}/lib -L${NCHOME}/lib before_script: - module purge - module load git # Extend the existing environment variables: - export PATH="${PYHOME}/bin:${PATH}" - - export LD_LIBRARY_PATH="${NFHOME}/lib:${NCHOME}/lib:${LD_LIBRARY_PATH-}" # Some tests require a large stack: - ulimit -s unlimited @@ -70,8 +58,8 @@ variables: - .common-levante variables: # Compiler flags used for ICON model: - FCFLAGS: -g -O2 -Mrecursive -Mallocatable=03 -Mstack_arrays -Minfo=accel,inline -acc=gpu,verystrict -gpu=cc80,cuda11.7 -DRTE_USE_${FPMODEL} - RTE_KERNELS: accel + FFLAGS: -g -O2 -Mrecursive -Mallocatable=03 -Mstack_arrays -Minfo=accel,inline -acc=gpu,verystrict -gpu=cc80,cuda12.6 + KERNEL_MODE: accel .nag-cpu: extends: @@ -80,17 +68,17 @@ variables: - .common-levante variables: # Compiler flags used for ICON model: - FCFLAGS: -Wc=/sw/spack-levante/gcc-11.2.0-bcn7mb/bin/gcc -f2008 -colour -w=uep -g -gline -O0 -float-store -nan -Wc,-g -Wc,-pipe -Wc,--param,max-vartrack-size=200000000 -Wc,-mno-fma -C=all -DRTE_USE_CBOOL -DRTE_USE_${FPMODEL} + FFLAGS: -Wc=/sw/spack-levante/gcc-11.2.0-bcn7mb/bin/gcc -f2008 -colour -w=uep -g -gline -O0 -float-store -nan -Wc,-g -Wc,-pipe -Wc,--param,max-vartrack-size=200000000 -Wc,-mno-fma -C=all .nag-cpu-default: extends: .nag-cpu variables: - RTE_KERNELS: default + KERNEL_MODE: default .nag-cpu-accel: extends: .nag-cpu variables: - RTE_KERNELS: accel + KERNEL_MODE: accel nvhpc-gpu-openacc-DP: extends: diff --git a/.gitlab/lumi.yml b/.gitlab/lumi.yml index 9e2c1ad4c..4f351f2e2 100644 --- a/.gitlab/lumi.yml +++ b/.gitlab/lumi.yml @@ -38,10 +38,7 @@ variables: .cce: variables: - # Core variables: FC: ftn - # Convenience variables: - VERSION_FCFLAGS: -V COMPILER_MODULES: PrgEnv-cray cce/16.0.1 craype-x86-milan # @@ -91,8 +88,8 @@ setup-python: - .common-lumi variables: # Compiler flags used for ICON model: - FCFLAGS: -hacc -hadd_paren -Ktrap=divz,ovf,inv -hflex_mp=intolerant -hfp1 -g -DRTE_USE_${FPMODEL} - RTE_KERNELS: accel + FFLAGS: -hacc -hadd_paren -Ktrap=divz,ovf,inv -hflex_mp=intolerant -hfp1 -g + KERNEL_MODE: accel # Convenience variables: EXTRA_COMPILER_MODULES: craype-accel-amd-gfx90a rocm diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 000000000..d6bb9a102 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,157 @@ +cmake_minimum_required(VERSION 3.18) + +project( + rte-rrtmgp + VERSION 1.9.0 + LANGUAGES Fortran +) + +option(BUILD_TESTING "Build tests" OFF) + +option(RTE_ENABLE_SP "Enable single-precision floating-point model" OFF) + +set(PREFERRED_KERNEL_MODES "default" "accel" "extern") +set(KERNEL_MODE + "default" + CACHE STRING "Select the kernel mode: ${PREFERRED_KERNEL_MODES}" +) +set_property(CACHE KERNEL_MODE PROPERTY STRINGS ${PREFERRED_KERNEL_MODES}) + +list(APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake) + +# GNUInstallDirs issues a warning if CMAKE_SIZEOF_VOID_P is not defined, which +# is the case with NAG. One way to circumvent that is to enable C language for +# the project: +if(CMAKE_Fortran_COMPILER_ID STREQUAL NAG) + enable_language(C) +endif() +include(GNUInstallDirs) + +add_compile_options( + $<$:-ffree-line-length-none> +) + +set(CMAKE_Fortran_MODULE_DIRECTORY ${PROJECT_BINARY_DIR}/modules) + +add_subdirectory(rte-kernels) +add_subdirectory(rte-frontend) +add_subdirectory(rrtmgp-kernels) +add_subdirectory(rrtmgp-frontend) + +include(CTest) +if(BUILD_TESTING) + find_package(Python3 REQUIRED COMPONENTS Interpreter) + + if(NOT _RTE_RRTMGP_HAVE_PY_PACKAGES) + include(CheckPython3Package) + check_python3_package(numpy) + check_python3_package( + "netCDF4 or h5netcdf+scipy" + CODE "try: + import netCDF4 +except: + import h5netcdf + import scipy +" + ) + check_python3_package( + "xarray>=0.12.2" + CODE "import xarray +exit(tuple(map(int, xarray.__version__.split('.'))) < (0, 12, 2))" + ) + check_python3_package(dask.array) + set(_RTE_RRTMGP_HAVE_PY_PACKAGES + TRUE + CACHE INTERNAL + "RTE-RRTMGP found all Python packages required for testing" + ) + endif() + + find_package(NetCDF_Fortran REQUIRED) + + if(RRTMGP_DATA) + add_test(NAME fetch_rrtmgp_data COMMAND ${CMAKE_COMMAND} -E true) + message( + NOTICE + "Using an external dataset from ${RRTMGP_DATA}: the data files will not be installed" + ) + else() + set(RRTMGP_DATA "${PROJECT_BINARY_DIR}/rrtmgp-data") + + include(ExternalProject) + ExternalProject_Add( + rrtmgp-data + GIT_REPOSITORY https://github.com/earth-system-radiation/rrtmgp-data.git + GIT_TAG "v1.8.2" + GIT_SHALLOW True + EXCLUDE_FROM_ALL True + PREFIX rrtmgp-data-cmake + SOURCE_DIR ${RRTMGP_DATA} + CONFIGURE_COMMAND "" + BUILD_COMMAND "" + INSTALL_COMMAND "" + ) + + set(fetch_rrtmgp_data_command + ${CMAKE_COMMAND} --build ${CMAKE_CURRENT_BINARY_DIR} --config + "$" --target rrtmgp-data + ) + + add_test(NAME fetch_rrtmgp_data COMMAND ${fetch_rrtmgp_data_command}) + + install(CODE "execute_process(COMMAND ${fetch_rrtmgp_data_command})") + install( + FILES # cmake-format: sort + ${RRTMGP_DATA}/rrtmgp-aerosols-merra-lw.nc + ${RRTMGP_DATA}/rrtmgp-aerosols-merra-sw.nc + ${RRTMGP_DATA}/rrtmgp-clouds-lw.nc + ${RRTMGP_DATA}/rrtmgp-clouds-sw.nc + ${RRTMGP_DATA}/rrtmgp-gas-lw-g128.nc + ${RRTMGP_DATA}/rrtmgp-gas-lw-g256.nc + ${RRTMGP_DATA}/rrtmgp-gas-sw-g112.nc + ${RRTMGP_DATA}/rrtmgp-gas-sw-g224.nc + DESTINATION ${CMAKE_INSTALL_DATAROOTDIR}/rte-rrtmgp/ + ) + endif() + + set_tests_properties( + fetch_rrtmgp_data PROPERTIES FIXTURES_SETUP fetch_rrtmgp_data + ) + + add_subdirectory(examples) + add_subdirectory(tests) +else() + # Allow for 'make test' even if the tests are disabled: + enable_testing() +endif() + +export( + EXPORT rte-rrtmgp-targets FILE ${PROJECT_BINARY_DIR}/rte-rrtmgp-targets.cmake +) + +include(CMakePackageConfigHelpers) +configure_package_config_file( + ${PROJECT_SOURCE_DIR}/cmake/config.cmake.in + ${PROJECT_BINARY_DIR}/rte-rrtmgp-config.cmake + INSTALL_DESTINATION ${CMAKE_INSTALL_LIBDIR}/rte-rrtmgp/cmake + NO_SET_AND_CHECK_MACRO NO_CHECK_REQUIRED_COMPONENTS_MACRO +) + +write_basic_package_version_file( + ${PROJECT_BINARY_DIR}/rte-rrtmgp-config-version.cmake + VERSION ${PROJECT_VERSION} + COMPATIBILITY SameMinorVersion +) + +install(DIRECTORY ${CMAKE_Fortran_MODULE_DIRECTORY}/ TYPE INCLUDE) + +install( + EXPORT rte-rrtmgp-targets + DESTINATION ${CMAKE_INSTALL_LIBDIR}/rte-rrtmgp/cmake +) + +install( + FILES ${PROJECT_BINARY_DIR}/rte-rrtmgp-config.cmake + ${PROJECT_BINARY_DIR}/rte-rrtmgp-config-version.cmake + DESTINATION ${CMAKE_INSTALL_LIBDIR}/rte-rrtmgp/cmake +) diff --git a/Makefile b/Makefile deleted file mode 100644 index 96c1e2edd..000000000 --- a/Makefile +++ /dev/null @@ -1,28 +0,0 @@ -# -# Top-level Makefile -# -.PHONY: libs tests check -all: libs tests check - -libs: - $(MAKE) -C build $@ - -tests: - $(MAKE) -C tests $@ - $(MAKE) -C examples/rfmip-clear-sky $@ - $(MAKE) -C examples/all-sky $@ - -check: - $(MAKE) -C tests $@ - $(MAKE) -C examples/rfmip-clear-sky $@ - $(MAKE) -C examples/all-sky $@ - -docs: - @cd doc; ./build_documentation.sh - -clean: - $(MAKE) -C build $@ - $(MAKE) -C tests $@ - $(MAKE) -C examples/rfmip-clear-sky $@ - $(MAKE) -C examples/all-sky $@ - rm -rf public diff --git a/build/Makefile b/build/Makefile deleted file mode 100644 index efd0341d2..000000000 --- a/build/Makefile +++ /dev/null @@ -1,68 +0,0 @@ -#!/usr/bin/env make - -RTE_DIR = ../rte-frontend -GAS_OPTICS_DIR = ../gas-optics -RRTMGP_DIR = ../rrtmgp-frontend -RTE_KERNEL_DIR = ../rte-kernels -RRTMGP_KERNEL_DIR = ../rrtmgp-kernels -# -# Compiler variables FC, FCFLAGS must be set in the environment -# -# Make all the libraries though we'll only use the interface + kernels -all: librte.a librrtmgp.a \ - librtekernels.a librtef.a librrtmgpkernels.a librrtmgpf.a -separate-libs: librtekernels.a librtef.a librrtmgpkernels.a librrtmgpf.a -libs: all - -COMPILE = $(FC) $(FCFLAGS) $(FCINCLUDE) -c -%.o: %.F90 - $(COMPILE) $< - -include $(RTE_DIR)/Make.depends -include $(RRTMGP_DIR)/Make.depends -include $(RTE_KERNEL_DIR)/Make.depends -include $(RRTMGP_KERNEL_DIR)/Make.depends - -# -# If using OpenACC/OpenMP files in *-kernels/accel take precendence -# -ifeq ($(RTE_KERNELS), accel) - VPATH = $(RTE_KERNEL_DIR)/accel:$(RRTMGP_KERNEL_DIR)/accel -endif -# -# If using external libraries just compile the interfaces -# -ifeq ($(RTE_KERNELS), extern) - VPATH = $(RTE_KERNEL_DIR)/api:$(RRTMGP_KERNEL_DIR)/api -endif -VPATH += $(RTE_DIR):$(RTE_KERNEL_DIR):$(RRTMGP_DIR):$(RRTMGP_KERNEL_DIR):$(GAS_OPTICS_DIR) - -# -# Complete library - kernels plus Fortran front end -# -librte.a: $(RTE_FORTRAN_KERNELS) $(RTE_FORTRAN_INTERFACE) - ar -rvs librte.a $(RTE_FORTRAN_KERNELS) $(RTE_FORTRAN_INTERFACE) -# -# Library with just the kernels... -# -librtekernels.a: $(RTE_FORTRAN_KERNELS) - ar -rvs librtekernels.a $(RTE_FORTRAN_KERNELS) -# -# ... and just the Fortran front-end -# -librtef.a: $(RTE_FORTRAN_INTERFACE) - ar -rvs librtef.a $(RTE_FORTRAN_INTERFACE) -# -# As with RTE, libraries with Fortran front-end and kernels, separate and combined -# -librrtmgp.a: $(RRTMGP_FORTRAN_KERNELS) $(RRTMGP_FORTRAN_INTERFACE) - ar -rvs librrtmgp.a $(RRTMGP_FORTRAN_KERNELS) $(RRTMGP_FORTRAN_INTERFACE) - -librrtmgpkernels.a: $(RRTMGP_FORTRAN_KERNELS) - ar -rvs librrtmgpkernels.a $(RRTMGP_FORTRAN_KERNELS) - -librrtmgpf.a: $(RRTMGP_FORTRAN_INTERFACE) - ar -rvs librrtmgpf.a $(RRTMGP_FORTRAN_INTERFACE) - -clean: - rm -f *.optrpt *.mod *.o lib*.a diff --git a/cmake/CheckFortranNeedsCBool.cmake b/cmake/CheckFortranNeedsCBool.cmake new file mode 100644 index 000000000..1780d21e6 --- /dev/null +++ b/cmake/CheckFortranNeedsCBool.cmake @@ -0,0 +1,72 @@ +# ~~~ +# check_fortran_needs_c_bool() +# ~~~ +# Checks whether the Fortran compiler requires the c_bool kind of the logical +# type. Sets the cache to the result of the check. +# +function(check_fortran_needs_cbool var) + if(DEFINED ${var}) + return() + endif() + + if(NOT CMAKE_REQUIRED_QUIET) + message(CHECK_START "Checking if the Fortran compiler requires C_BOOL") + endif() + + set(CMAKE_TRY_COMPILE_TARGET_TYPE STATIC_LIBRARY) + + set(check_source_code + " + subroutine conftest_foo(a) bind(C) + use iso_c_binding + implicit none +#ifdef RTE_USE_CBOOL + integer, parameter :: wl = c_bool +#else + integer, parameter :: wl = kind(.true.) +#endif + logical(wl) :: a + end subroutine +" + ) + + set(check_source_file + "${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/CMakeTmp/src.F90" + ) + + unset(check_result) + file(WRITE "${check_source_file}" "${check_source_code}") + try_compile(try_result "${CMAKE_BINARY_DIR}" "${check_source_file}") + if(try_result) + set(check_result FALSE) + else() + file(WRITE "${check_source_file}" "${check_source_code}") + try_compile( + try_result "${CMAKE_BINARY_DIR}" + "${check_source_file}" + COMPILE_DEFINITIONS "-DRTE_USE_CBOOL" + ) + if(try_result) + set(check_result TRUE) + endif() + endif() + + if(NOT CMAKE_REQUIRED_QUIET) + if(NOT DEFINED check_result) + message(CHECK_FAIL "unknown (assuming no)") + elseif(check_result) + message(CHECK_PASS "yes") + else() + message(CHECK_PASS "no") + endif() + endif() + + if(NOT DEFINED check_result) + set(check_result FALSE) + endif() + + set(${var} + "${check_result}" + CACHE BOOL "Whether the Fortran compiler requires CBOOL type" + ) +endfunction() diff --git a/cmake/CheckPython3Package.cmake b/cmake/CheckPython3Package.cmake new file mode 100644 index 000000000..c70241edd --- /dev/null +++ b/cmake/CheckPython3Package.cmake @@ -0,0 +1,43 @@ +# ~~~ +# check_python3_package( +# [CODE ]) +# ~~~ +# Checks whether the Python3 is available by running the with +# ${Python3_EXECUTABLE} (defaults to 'import '). Fails the +# configuration if the result is negative. +# +function(check_python3_package package) + cmake_parse_arguments(PARSE_ARGV 1 arg "" "CODE" "") + + if(DEFINED ${var}) + return() + endif() + + if(NOT CMAKE_REQUIRED_QUIET) + message( + CHECK_START "Checking if the Python3 package ${package} is available" + ) + endif() + + if(NOT arg_CODE) + set(arg_CODE "import ${package}") + endif() + + execute_process( + COMMAND ${Python3_EXECUTABLE} -c "${arg_CODE}" + RESULT_VARIABLE exit_status + OUTPUT_QUIET ERROR_QUIET + ) + + if(NOT CMAKE_REQUIRED_QUIET) + if(NOT exit_status) + message(CHECK_PASS "yes") + else() + message(CHECK_FAIL "no") + endif() + endif() + + if(exit_status) + message(FATAL_ERROR "Required Python3 package ${package} is not available") + endif() +endfunction() diff --git a/cmake/FindNetCDF_Fortran.cmake b/cmake/FindNetCDF_Fortran.cmake new file mode 100644 index 000000000..7dd8da726 --- /dev/null +++ b/cmake/FindNetCDF_Fortran.cmake @@ -0,0 +1,19 @@ +find_library( + NetCDF_Fortran_LIBRARY + NAMES netcdff + DOC "NetCDF-Fortran library" +) +mark_as_advanced(NetCDF_Fortran_LIBRARY) + +find_path( + NetCDF_Fortran_INCLUDE_DIR + NAMES netcdf.mod NETCDF.mod + DOC "NetCDF_Fortran include directory" +) +mark_as_advanced(NetCDF_Fortran_INCLUDE_DIR) + +include(FindPackageHandleStandardArgs) +find_package_handle_standard_args( + NetCDF_Fortran + REQUIRED_VARS NetCDF_Fortran_LIBRARY NetCDF_Fortran_INCLUDE_DIR +) diff --git a/cmake/config.cmake.in b/cmake/config.cmake.in new file mode 100644 index 000000000..b2de3db4d --- /dev/null +++ b/cmake/config.cmake.in @@ -0,0 +1,6 @@ +@PACKAGE_INIT@ + +include(${CMAKE_CURRENT_LIST_DIR}/rte-rrtmgp-targets.cmake) + +include(FindPackageHandleStandardArgs) +find_package_handle_standard_args(rte-rrtmgp REQUIRED_VARS rte-rrtmgp_DIR) diff --git a/doc/jekyll_site/how-tos/build-and-test.md b/doc/jekyll_site/how-tos/build-and-test.md index 07b4dc186..bec95a35b 100644 --- a/doc/jekyll_site/how-tos/build-and-test.md +++ b/doc/jekyll_site/how-tos/build-and-test.md @@ -5,22 +5,26 @@ title: "How to build and run tests" How to build the libraries, tests, and examples, run the tests, and verify the results ## In a nutshell -In the root directory: -- `make libs` makes the RTE and RRTMGP libraries, the unit tests, and the examples -- `make tests` runs the tests -- `make check` uses Python to verify results against reference calculations -- `make` invoked without a target in the top level attempts all three steps. +RTE+RRTMGP uses `CMake`. In the root directory: +`cmake -S . -B build` will guide you through configuration options. + +Environment variables can also be set and passed to `CMake` as in this example, which +builds and runs the tests: + +` +cmake -S . -B build \ + -DCMAKE_Fortran_COMPILER=$FC \ + -DCMAKE_Fortran_FLAGS=$FCFLAGS \ + -DRRTMGP_DATA_VERSION=$RRTMGP_DATA_VERSION \ + -DPRECISION=$FP_MODEL \ + -DUSE_C_BOOL=$RTE_CBOOL \ + -DKERNEL_MODE=$RTE_KERNELS \ + -DENABLE_TESTS=ON \ + -DFAILURE_THRESHOLD=$FAILURE_THRESHOLD +` Evaluating the results of the tests requires `Python` and the packages described in `environment*.yml`. -## Building and testing using the handbuilt Makefiles - -Before using the Makefiles supplied with the `RTE+RRTMGP` repository, the environment variables `FC` and -`FCFLAGS`, identifying the Fortran compiler and flags passed to it, need to be set. - -To build any of the examples in `examples/` or `tests` the locations of the C and Fortran netCDF libraries and the -location of the netCDF Fortran module file (`netcdf.mod`) must be in the search path. -Non-standard paths can also be added via macros `FCINCLUDE` and/or `LDFLAGS`. ## Building and testing using (Gnu) autotools diff --git a/environment-dev.yml b/environment-dev.yml index b9fa161d9..affe2b0b2 100644 --- a/environment-dev.yml +++ b/environment-dev.yml @@ -19,6 +19,9 @@ dependencies: - colorcet - gfortran - netcdf-fortran + - cmake + - cmake-format + - ninja variables: FC: gfortran # Debugging flags below diff --git a/environment-noplots.yml b/environment-noplots.yml index dfa6536fc..50a628def 100644 --- a/environment-noplots.yml +++ b/environment-noplots.yml @@ -2,10 +2,13 @@ # Python modules below are needed to run tests and check results # name: rte_rrtmgp_test_noplots - +channels: + - conda-forge + - nodefaults dependencies: - python=3.11 - netcdf4 - xarray - dask - numpy + - cmake diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt new file mode 100644 index 000000000..1b4a7c96a --- /dev/null +++ b/examples/CMakeLists.txt @@ -0,0 +1,32 @@ +if(RTE_ENABLE_SP) + set(default_FAILURE_THRESHOLD "3.5e-1") +else() + set(default_FAILURE_THRESHOLD "7.e-4") +endif() + +set(FAILURE_THRESHOLD + "${default_FAILURE_THRESHOLD}" + CACHE STRING "Default failure threshold for tests" +) + +set(CMAKE_Fortran_MODULE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/modules) + +add_library( + examples_utils STATIC # cmake-format: sort + mo_garand_atmos_io.F90 mo_load_coefficients.F90 mo_simple_netcdf.F90 +) + +target_include_directories( + examples_utils + PUBLIC + $:${CMAKE_Fortran_MODULE_DIRECTORY}>> + ${NetCDF_Fortran_INCLUDE_DIR} +) + +target_link_libraries( + examples_utils + PUBLIC rte-rrtmgp::rrtmgp rte-rrtmgp::rte ${NetCDF_Fortran_LIBRARY} +) + +add_subdirectory(all-sky) +add_subdirectory(rfmip-clear-sky) diff --git a/examples/all-sky/CMakeLists.txt b/examples/all-sky/CMakeLists.txt new file mode 100644 index 000000000..185900070 --- /dev/null +++ b/examples/all-sky/CMakeLists.txt @@ -0,0 +1,89 @@ +add_library( + all_sky_utils STATIC # cmake-format: sort + mo_load_aerosol_coefficients.F90 mo_load_cloud_coefficients.F90 +) + +target_link_libraries(all_sky_utils PUBLIC examples_utils) + +add_executable(rrtmgp_allsky rrtmgp_allsky.F90) +target_link_libraries(rrtmgp_allsky PRIVATE all_sky_utils) + +add_test( + NAME run_allsky_lw + COMMAND + rrtmgp_allsky 24 72 1 rrtmgp-allsky-lw.nc + ${RRTMGP_DATA}/rrtmgp-gas-lw-g256.nc ${RRTMGP_DATA}/rrtmgp-clouds-lw.nc + ${RRTMGP_DATA}/rrtmgp-aerosols-merra-lw.nc +) +set_tests_properties( + run_allsky_lw + PROPERTIES FIXTURES_REQUIRED fetch_rrtmgp_data FIXTURES_SETUP run_allsky +) + +add_test( + NAME run_allsky_sw + COMMAND + rrtmgp_allsky 24 72 1 rrtmgp-allsky-sw.nc + ${RRTMGP_DATA}/rrtmgp-gas-sw-g224.nc ${RRTMGP_DATA}/rrtmgp-clouds-sw.nc + ${RRTMGP_DATA}/rrtmgp-aerosols-merra-sw.nc +) +set_tests_properties( + run_allsky_sw + PROPERTIES FIXTURES_REQUIRED fetch_rrtmgp_data FIXTURES_SETUP run_allsky +) + +add_test( + NAME check_allsky_lw_sw + COMMAND + ${Python3_EXECUTABLE} ${CMAKE_SOURCE_DIR}/examples/compare-to-reference.py + --ref_dir ${RRTMGP_DATA}/examples/all-sky/reference --tst_dir + ${CMAKE_CURRENT_BINARY_DIR} --variables lw_flux_up lw_flux_dn sw_flux_up + sw_flux_dn sw_flux_dir --file_names rrtmgp-allsky-lw.nc rrtmgp-allsky-sw.nc + --failure_threshold ${FAILURE_THRESHOLD} +) +set_tests_properties( + check_allsky_lw_sw + PROPERTIES FIXTURES_REQUIRED "fetch_rrtmgp_data;run_allsky" +) + +add_test( + NAME run_allsky_no_aerosols_lw + COMMAND + rrtmgp_allsky 24 72 1 rrtmgp-allsky-lw-no-aerosols.nc + ${RRTMGP_DATA}/rrtmgp-gas-lw-g256.nc ${RRTMGP_DATA}/rrtmgp-clouds-lw.nc +) +set_tests_properties( + run_allsky_no_aerosols_lw + PROPERTIES FIXTURES_REQUIRED + fetch_rrtmgp_data + FIXTURES_SETUP + run_allsky_no_aerosols +) + +add_test( + NAME run_allsky_no_aerosols_sw + COMMAND + rrtmgp_allsky 24 72 1 rrtmgp-allsky-sw-no-aerosols.nc + ${RRTMGP_DATA}/rrtmgp-gas-sw-g224.nc ${RRTMGP_DATA}/rrtmgp-clouds-sw.nc +) +set_tests_properties( + run_allsky_no_aerosols_sw + PROPERTIES FIXTURES_REQUIRED + fetch_rrtmgp_data + FIXTURES_SETUP + run_allsky_no_aerosols +) + +add_test( + NAME check_allsky_no_aerosols_lw_sw + COMMAND + ${Python3_EXECUTABLE} ${CMAKE_SOURCE_DIR}/examples/compare-to-reference.py + --ref_dir ${RRTMGP_DATA}/examples/all-sky/reference --tst_dir + ${CMAKE_CURRENT_BINARY_DIR} --variables lw_flux_up lw_flux_dn sw_flux_up + sw_flux_dn sw_flux_dir --file_names rrtmgp-allsky-lw-no-aerosols.nc + rrtmgp-allsky-sw-no-aerosols.nc --failure_threshold ${FAILURE_THRESHOLD} +) +set_tests_properties( + check_allsky_no_aerosols_lw_sw + PROPERTIES FIXTURES_REQUIRED "fetch_rrtmgp_data;run_allsky_no_aerosols" +) diff --git a/examples/all-sky/Makefile b/examples/all-sky/Makefile deleted file mode 100644 index d404181c1..000000000 --- a/examples/all-sky/Makefile +++ /dev/null @@ -1,62 +0,0 @@ -#!/usr/bin/env make -# -# Location of RTE+RRTMGP libraries, module files. -# -RRTMGP_BUILD = $(RRTMGP_ROOT)/build -# -# RRTMGP library, module files -# -LDFLAGS += -L$(RRTMGP_BUILD) -LIBS += -lrrtmgp -lrte -FCINCLUDE += -I$(RRTMGP_BUILD) - -# netcdf Fortran module files has to be in the search path or added via environment variable FCINCLUDE e.g. -#FCINCLUDE += -I$(NFHOME)/include - -# netcdf C and Fortran libraries have to be in the search path or added via environment variable LDFLAGS e.g. -#LDFLAGS += -L$(NFHOME)/lib -L$(NCHOME)/lib -LIBS += -lnetcdff -lnetcdf - -VPATH = ../:$(RRTMGP_ROOT)/rrtmgp-frontend # Needed for cloud_optics and aerosol_optics - -# Compilation rules -%.o: %.F90 - $(FC) $(FCFLAGS) $(FCINCLUDE) -c $< -%: %.o - $(FC) $(FCFLAGS) -o $@ $^ $(LDFLAGS) $(LIBS) - -# -# Extra sources -- extensions to RRTMGP classes, shared infrastructure, local sources -# -ADDITIONS = mo_load_coefficients.o mo_simple_netcdf.o -ADDITIONS += mo_cloud_optics_rrtmgp.o mo_load_cloud_coefficients.o -ADDITIONS += mo_aerosol_optics_rrtmgp_merra.o mo_load_aerosol_coefficients.o - -# -# Targets -# -all: rrtmgp_allsky - -rrtmgp_allsky: $(ADDITIONS) rrtmgp_allsky.o - -rrtmgp_allsky.o: $(ADDITIONS) rrtmgp_allsky.F90 - -mo_cloud_optics_rrtmgp.o: mo_cloud_optics_rrtmgp.F90 -mo_aerosol_optics_rrtmgp_merra.o: mo_aerosol_optics_rrtmgp_merra.F90 -mo_load_coefficients.o: mo_simple_netcdf.o mo_load_coefficients.F90 -mo_load_cloud_coefficients.o: mo_simple_netcdf.o mo_cloud_optics_rrtmgp.o mo_load_cloud_coefficients.F90 -mo_load_aerosol_coefficients.o: mo_simple_netcdf.o mo_aerosol_optics_rrtmgp_merra.o mo_load_aerosol_coefficients.F90 - -tests: rrtmgp_allsky - $(RUN_CMD) bash all_tests.sh - -check: - $${PYTHON-python} ${RRTMGP_ROOT}/examples/compare-to-reference.py --ref_dir ${RRTMGP_DATA}/examples/all-sky/reference --tst_dir ${RRTMGP_ROOT}/examples/all-sky \ - --var lw_flux_up lw_flux_dn sw_flux_up sw_flux_dn sw_flux_dir \ - --file rrtmgp-allsky-lw.nc rrtmgp-allsky-sw.nc - $${PYTHON-python} ${RRTMGP_ROOT}/examples/compare-to-reference.py --ref_dir ${RRTMGP_DATA}/examples/all-sky/reference --tst_dir ${RRTMGP_ROOT}/examples/all-sky \ - --var lw_flux_up lw_flux_dn sw_flux_up sw_flux_dn sw_flux_dir \ - --file rrtmgp-allsky-lw-no-aerosols.nc rrtmgp-allsky-sw-no-aerosols.nc - -clean: - -rm rrtmgp_allsky *.o *.optrpt ../*.optrpt *.mod *.nc diff --git a/examples/all-sky/all_tests.sh b/examples/all-sky/all_tests.sh deleted file mode 100644 index 61a2bceb7..000000000 --- a/examples/all-sky/all_tests.sh +++ /dev/null @@ -1,9 +0,0 @@ -set -eux -./rrtmgp_allsky 24 72 1 rrtmgp-allsky-lw.nc \ - ${RRTMGP_DATA}/rrtmgp-gas-lw-g256.nc ${RRTMGP_DATA}/rrtmgp-clouds-lw.nc ${RRTMGP_DATA}/rrtmgp-aerosols-merra-lw.nc -./rrtmgp_allsky 24 72 1 rrtmgp-allsky-sw.nc \ - ${RRTMGP_DATA}/rrtmgp-gas-sw-g224.nc ${RRTMGP_DATA}/rrtmgp-clouds-sw.nc ${RRTMGP_DATA}/rrtmgp-aerosols-merra-sw.nc -./rrtmgp_allsky 24 72 1 rrtmgp-allsky-lw-no-aerosols.nc \ - ${RRTMGP_DATA}/rrtmgp-gas-lw-g256.nc ${RRTMGP_DATA}/rrtmgp-clouds-lw.nc -./rrtmgp_allsky 24 72 1 rrtmgp-allsky-sw-no-aerosols.nc \ - ${RRTMGP_DATA}/rrtmgp-gas-sw-g224.nc ${RRTMGP_DATA}/rrtmgp-clouds-sw.nc \ No newline at end of file diff --git a/examples/all-sky/make_problem_size_loop.py b/examples/all-sky/make_problem_size_loop.py index 17483e71f..8b4cdee98 100644 --- a/examples/all-sky/make_problem_size_loop.py +++ b/examples/all-sky/make_problem_size_loop.py @@ -60,5 +60,5 @@ for l in args.nlay: for i in args.ncol: print(f"{args.executable} {i:6d} {l:4d} {args.nloops:3d} " + - f"{args.output_file} {args.k_distribution}" + - f"{args.cloud_optics} {args.aerosol_optics}") + f"{args.output_file} {args.k_distribution} " + + f"{args.cloud_optics} {args.aerosol_optics} ") diff --git a/examples/all-sky/rrtmgp_allsky.F90 b/examples/all-sky/rrtmgp_allsky.F90 index e25ec0bf3..f891a194e 100644 --- a/examples/all-sky/rrtmgp_allsky.F90 +++ b/examples/all-sky/rrtmgp_allsky.F90 @@ -89,10 +89,11 @@ program rte_rrtmgp_allsky ! ! Inputs to RRTMGP ! - logical :: top_at_1, is_sw, is_lw + logical :: is_sw, is_lw integer :: nbnd, ngpt integer :: icol, ilay, ibnd, iloop, igas + logical :: top_is_at_1 ! CCE OMP workaround character(len=8) :: char_input integer :: nUserArgs, nloops, ncol, nlay @@ -191,7 +192,6 @@ program rte_rrtmgp_allsky ! nbnd = k_dist%get_nband() ngpt = k_dist%get_ngpt() - top_at_1 = p_lay(1, 1) < p_lay(1, nlay) ! ---------------------------------------------------------------------------- ! LW calculations neglect scattering; SW calculations use the 2-stream approximation @@ -249,9 +249,10 @@ program rte_rrtmgp_allsky !$acc enter data create (t_sfc, emis_sfc) !$omp target enter data map(alloc:t_sfc, emis_sfc) ! Surface temperature + top_is_at_1 = atmos%top_is_at_1() ! CCE OMP workaround !$acc kernels !$omp target - t_sfc = t_lev(1, merge(nlay+1, 1, top_at_1)) + t_sfc = t_lev(1, merge(nlay+1, 1, top_is_at_1)) emis_sfc = 0.98_wp !$acc end kernels !$omp end target @@ -322,9 +323,9 @@ program rte_rrtmgp_allsky tlev = t_lev)) if(do_clouds) call stop_on_err(clouds%increment(atmos)) if(do_aerosols) call stop_on_err(aerosols%increment(atmos)) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - emis_sfc, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + emis_sfc, & fluxes)) !$acc end data !$omp end target data @@ -347,8 +348,7 @@ program rte_rrtmgp_allsky call stop_on_err(aerosols%delta_scale()) call stop_on_err(aerosols%increment(atmos)) end if - call stop_on_err(rte_sw(atmos, top_at_1, & - mu0, toa_flux, & + call stop_on_err(rte_sw(atmos, mu0, toa_flux, & sfc_alb_dir, sfc_alb_dif, & fluxes)) !$acc end data diff --git a/examples/rfmip-clear-sky/CMakeLists.txt b/examples/rfmip-clear-sky/CMakeLists.txt new file mode 100644 index 000000000..37c69177a --- /dev/null +++ b/examples/rfmip-clear-sky/CMakeLists.txt @@ -0,0 +1,83 @@ +add_library(rfmip_clear_utils STATIC mo_rfmip_io.F90) +target_link_libraries(rfmip_clear_utils PUBLIC examples_utils) + +foreach( + test_executable IN + ITEMS # cmake-format: sort + rrtmgp_rfmip_lw + rrtmgp_rfmip_sw +) + add_executable(${test_executable} ${test_executable}.F90) + target_link_libraries(${test_executable} PRIVATE rfmip_clear_utils) +endforeach() + +set(inoutputs + rld_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc + rlu_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc + rsd_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc + rsu_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc +) + +list( + TRANSFORM inoutputs + PREPEND ${RRTMGP_DATA}/examples/rfmip-clear-sky/inputs/ + OUTPUT_VARIABLE inputs +) + +# The tests write to the input files, therefore we copy them: +add_test( + NAME copy_rrtmgp_rfmip_inputs + COMMAND + ${CMAKE_COMMAND} -E copy_if_different ${inputs} + ${CMAKE_CURRENT_BINARY_DIR}/ +) +set_tests_properties( + copy_rrtmgp_rfmip_inputs + PROPERTIES FIXTURES_REQUIRED + fetch_rrtmgp_data + FIXTURES_SETUP + copy_rrtmgp_rfmip_inputs +) + +add_test( + NAME run_rrtmgp_rfmip_lw + COMMAND + rrtmgp_rfmip_lw 8 + ${RRTMGP_DATA}/examples/rfmip-clear-sky/inputs/multiple_input4MIPs_radiation_RFMIP_UColorado-RFMIP-1-2_none.nc + ${RRTMGP_DATA}/rrtmgp-gas-lw-g256.nc +) +set_tests_properties( + run_rrtmgp_rfmip_lw + PROPERTIES FIXTURES_REQUIRED + "fetch_rrtmgp_data;copy_rrtmgp_rfmip_inputs" + FIXTURES_SETUP + run_rrtmgp_rfmip +) + +add_test( + NAME run_rrtmgp_rfmip_sw + COMMAND + rrtmgp_rfmip_sw 8 + ${RRTMGP_DATA}/examples/rfmip-clear-sky/inputs/multiple_input4MIPs_radiation_RFMIP_UColorado-RFMIP-1-2_none.nc + ${RRTMGP_DATA}/rrtmgp-gas-sw-g224.nc +) +set_tests_properties( + run_rrtmgp_rfmip_sw + PROPERTIES FIXTURES_REQUIRED + "fetch_rrtmgp_data;copy_rrtmgp_rfmip_inputs" + FIXTURES_SETUP + run_rrtmgp_rfmip +) + +add_test( + NAME check_rfmip_lw_sw + COMMAND + ${Python3_EXECUTABLE} ${CMAKE_SOURCE_DIR}/examples/compare-to-reference.py + --ref_dir ${RRTMGP_DATA}/examples/rfmip-clear-sky/reference --tst_dir + ${CMAKE_CURRENT_BINARY_DIR} --variables rld rlu rsd rsu --file_names + ${inoutputs} --failure_threshold ${FAILURE_THRESHOLD} +) +set_tests_properties( + check_rfmip_lw_sw + PROPERTIES FIXTURES_REQUIRED "fetch_rrtmgp_data;run_rrtmgp_rfmip" +) diff --git a/examples/rfmip-clear-sky/Makefile b/examples/rfmip-clear-sky/Makefile deleted file mode 100644 index bad1bc80f..000000000 --- a/examples/rfmip-clear-sky/Makefile +++ /dev/null @@ -1,66 +0,0 @@ -#!/usr/bin/env make -# -# Location of RTE+RRTMGP libraries, module files. -# -RRTMGP_BUILD = $(RRTMGP_ROOT)/build -# -# RRTMGP library, module files -# -LDFLAGS += -L$(RRTMGP_BUILD) -LIBS += -lrrtmgp -lrte -FCINCLUDE += -I$(RRTMGP_BUILD) - -# netcdf Fortran module files has to be in the search path or added via environment variable FCINCLUDE e.g. -#FCINCLUDE += -I$(NFHOME)/include - -# netcdf C and Fortran libraries have to be in the search path or added via environment variable LDFLAGS e.g. -#LDFLAGS += -L$(NFHOME)/lib -L$(NCHOME)/lib -LIBS += -lnetcdff -lnetcdf - -VPATH = ../ - -# Compilation rules -%.o: %.F90 - $(FC) $(FCFLAGS) $(FCINCLUDE) -c $< - -%: %.o - $(FC) $(FCFLAGS) -o $@ $^ $(LDFLAGS) $(LIBS) - -# Required netCDF files are in $RRTMGP_DATA -%.nc: - cp ${RRTMGP_DATA}/examples/rfmip-clear-sky/inputs/$@ . - - -# -# Ancillary codes -# -ADDITIONS = mo_simple_netcdf.o mo_rfmip_io.o mo_load_coefficients.o - -all: rrtmgp_rfmip_lw rrtmgp_rfmip_sw - -rrtmgp_rfmip_lw: rrtmgp_rfmip_lw.o $(ADDITIONS) - -rrtmgp_rfmip_lw.o: rrtmgp_rfmip_lw.F90 $(ADDITIONS) - -rrtmgp_rfmip_sw: rrtmgp_rfmip_sw.o $(ADDITIONS) - -rrtmgp_rfmip_sw.o: rrtmgp_rfmip_sw.F90 $(ADDITIONS) - -mo_rfmip_io.o: mo_rfmip_io.F90 mo_simple_netcdf.o - -mo_load_coefficients.o: mo_load_coefficients.F90 mo_simple_netcdf.o - -tests: rrtmgp_rfmip_lw rrtmgp_rfmip_sw \ - multiple_input4MIPs_radiation_RFMIP_UColorado-RFMIP-1-2_none.nc \ - rld_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc rlu_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc \ - rsd_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc rsu_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc - $(RUN_CMD) ./rrtmgp_rfmip_lw 8 multiple_input4MIPs_radiation_RFMIP_UColorado-RFMIP-1-2_none.nc ${RRTMGP_DATA}/rrtmgp-gas-lw-g256.nc - $(RUN_CMD) ./rrtmgp_rfmip_sw 8 multiple_input4MIPs_radiation_RFMIP_UColorado-RFMIP-1-2_none.nc ${RRTMGP_DATA}/rrtmgp-gas-sw-g224.nc - -check: - $${PYTHON-python} ${RRTMGP_ROOT}/examples/compare-to-reference.py \ - --ref_dir ${RRTMGP_DATA}/examples/rfmip-clear-sky/reference --tst_dir ${RRTMGP_ROOT}/examples/rfmip-clear-sky \ - --var rld rlu rsd rsu --file r??_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc - -clean: - -rm rrtmgp_rfmip_sw rrtmgp_rfmip_lw *.o *.mod *.optrpt *.nc diff --git a/examples/rfmip-clear-sky/rrtmgp_rfmip_lw.F90 b/examples/rfmip-clear-sky/rrtmgp_rfmip_lw.F90 index 607e372ef..07fac9c26 100644 --- a/examples/rfmip-clear-sky/rrtmgp_rfmip_lw.F90 +++ b/examples/rfmip-clear-sky/rrtmgp_rfmip_lw.F90 @@ -90,7 +90,6 @@ program rrtmgp_rfmip_lw character(len=132) :: kdist_file = 'coefficients_lw.nc' character(len=132) :: flxdn_file, flxup_file integer :: nargs, ncol, nlay, nbnd, nexp, nblocks, block_size, forcing_index, physics_index, n_quad_angles = 1 - logical :: top_at_1 integer :: b, icol, ibnd character(len=4) :: block_size_char, forcing_index_char = '1', physics_index_char = '1' @@ -168,10 +167,6 @@ program rrtmgp_rfmip_lw ! Allocation on assignment within reading routines ! call read_and_block_pt(rfmip_file, block_size, p_lay, p_lev, t_lay, t_lev) - ! - ! Are the arrays ordered in the vertical with 1 at the top or the bottom of the domain? - ! - top_at_1 = p_lay(1, 1, 1) < p_lay(1, nlay, 1) ! ! Read the gas concentrations and surface properties @@ -193,7 +188,8 @@ program rrtmgp_rfmip_lw ! is set to 10^-3 Pa. Here we pretend the layer is just a bit less deep. ! This introduces an error but shows input sanitizing. ! - if(top_at_1) then + ! Are the arrays ordered in the vertical with 1 at the top or the bottom of the domain? + if(p_lay(1, 1, 1) < p_lay(1, nlay, 1)) then p_lev(:,1,:) = k_dist%get_press_min() + epsilon(k_dist%get_press_min()) else p_lev(:,nlay+1,:) & @@ -256,7 +252,6 @@ program rrtmgp_rfmip_lw ! via ty_fluxes_broadband ! call stop_on_err(rte_lw(optical_props, & - top_at_1, & source, & sfc_emis_spec, & fluxes, n_gauss_angles = n_quad_angles)) diff --git a/examples/rfmip-clear-sky/rrtmgp_rfmip_sw.F90 b/examples/rfmip-clear-sky/rrtmgp_rfmip_sw.F90 index 547d9b4e3..e914f3f93 100644 --- a/examples/rfmip-clear-sky/rrtmgp_rfmip_sw.F90 +++ b/examples/rfmip-clear-sky/rrtmgp_rfmip_sw.F90 @@ -90,7 +90,6 @@ program rrtmgp_rfmip_sw character(len=132) :: kdist_file = 'coefficients_sw.nc' character(len=132) :: flxdn_file, flxup_file integer :: nargs, ncol, nlay, nbnd, ngpt, nexp, nblocks, block_size, forcing_index - logical :: top_at_1 integer :: b, icol, ibnd, igpt character(len=4) :: block_size_char, forcing_index_char = '1' @@ -166,10 +165,6 @@ program rrtmgp_rfmip_sw ! Allocation on assignment within reading routines ! call read_and_block_pt(rfmip_file, block_size, p_lay, p_lev, t_lay, t_lev) - ! - ! Are the arrays ordered in the vertical with 1 at the top or the bottom of the domain? - ! - top_at_1 = p_lay(1, 1, 1) < p_lay(1, nlay, 1) ! ! Read the gas concentrations and surface properties @@ -193,7 +188,10 @@ program rrtmgp_rfmip_sw ! is set to 10^-3 Pa. Here we pretend the layer is just a bit less deep. ! This introduces an error but shows input sanitizing. ! - if(top_at_1) then + ! + ! Are the arrays ordered in the vertical with 1 at the top or the bottom of the domain? + ! + if(p_lay(1, 1, 1) < p_lay(1, nlay, 1)) then p_lev(:,1,:) = k_dist%get_press_min() + epsilon(k_dist%get_press_min()) else p_lev(:,nlay+1,:) & @@ -296,7 +294,6 @@ program rrtmgp_rfmip_sw ! via ty_fluxes_broadband ! call stop_on_err(rte_sw(optical_props, & - top_at_1, & mu0, & toa_flux, & sfc_alb_spec, & diff --git a/extensions/mo_compute_bc.F90 b/extensions/mo_compute_bc.F90 index ef2a50b29..fbec1d974 100644 --- a/extensions/mo_compute_bc.F90 +++ b/extensions/mo_compute_bc.F90 @@ -166,7 +166,6 @@ function compute_bc(k_dist, & ! Compute fluxes ! error_msg = rte_lw(optical_props_1lay, & - top_at_1, & lw_sources_1lay, & lower_bc, fluxes_1lev) else @@ -191,8 +190,7 @@ function compute_bc(k_dist, & optical_props_1lay, & solar_src) error_msg = rte_sw(optical_props_1lay, & - top_at_1, mu0, & - solar_src, & + mu0, solar_src, & lower_bc, lower_bc, fluxes_1lev) endif end function diff --git a/extensions/mo_heating_rates.F90 b/extensions/mo_heating_rates.F90 index dc2ce15c7..9535ae4cd 100644 --- a/extensions/mo_heating_rates.F90 +++ b/extensions/mo_heating_rates.F90 @@ -71,7 +71,6 @@ function compute_heating_rate_solar_varmu0(flux_up, flux_dn, flux_dir, p_lev, mu ! --------- integer :: ncol, nlay, icol, ilay integer :: last_sunlight_layer(size(mu0, 1)) - logical(wl) :: top_at_1 ! --------- error_msg = "" ! diff --git a/extensions/mo_rrtmgp_clr_all_sky.F90 b/extensions/mo_rrtmgp_clr_all_sky.F90 index 755681cbb..eb122d73d 100644 --- a/extensions/mo_rrtmgp_clr_all_sky.F90 +++ b/extensions/mo_rrtmgp_clr_all_sky.F90 @@ -74,7 +74,6 @@ function rte_lw(k_dist, gas_concs, p_lay, t_lay, p_lev, & type(ty_source_func_lw) :: sources integer :: ncol, nlay, ngpt, nband, nstr - logical :: top_at_1 ! -------------------------------- ! Problem sizes ! @@ -85,11 +84,6 @@ function rte_lw(k_dist, gas_concs, p_lay, t_lay, p_lev, & ngpt = k_dist%get_ngpt() nband = k_dist%get_nband() - !$acc kernels copyout(top_at_1) - !$omp target map(from:top_at_1) - top_at_1 = p_lay(1, 1) < p_lay(1, nlay) - !$acc end kernels - !$omp end target ! ------------------------------------------------------------------------------------ ! Error checking @@ -161,8 +155,8 @@ function rte_lw(k_dist, gas_concs, p_lay, t_lay, p_lev, & if(present(aer_props)) error_msg = aer_props%increment(optical_props) if(error_msg /= '') return - error_msg = base_rte_lw(optical_props, top_at_1, sources, & - sfc_emis, clrsky_fluxes, & + error_msg = base_rte_lw(optical_props, sources, & + sfc_emis, clrsky_fluxes, & inc_flux, n_gauss_angles) if(error_msg /= '') return ! ------------------------------------------------------------------------------------ @@ -171,8 +165,8 @@ function rte_lw(k_dist, gas_concs, p_lay, t_lay, p_lev, & error_msg = cloud_props%increment(optical_props) if(error_msg /= '') return - error_msg = base_rte_lw(optical_props, top_at_1, sources, & - sfc_emis, allsky_fluxes, & + error_msg = base_rte_lw(optical_props, sources, & + sfc_emis, allsky_fluxes, & inc_flux, n_gauss_angles) call sources%finalize() @@ -207,7 +201,6 @@ function rte_sw(k_dist, gas_concs, p_lay, t_lay, p_lev, & class(ty_optical_props_arry), allocatable :: optical_props real(wp), dimension(:,:), allocatable :: toa_flux integer :: ncol, nlay, ngpt, nband, nstr - logical :: top_at_1 ! -------------------------------- ! Problem sizes ! @@ -218,12 +211,6 @@ function rte_sw(k_dist, gas_concs, p_lay, t_lay, p_lev, & ngpt = k_dist%get_ngpt() nband = k_dist%get_nband() - !$acc kernels copyout(top_at_1) - !$omp target map(from:top_at_1) - top_at_1 = p_lay(1, 1) < p_lay(1, nlay) - !$acc end kernels - !$omp end target - ! ------------------------------------------------------------------------------------ ! Error checking ! @@ -276,7 +263,7 @@ function rte_sw(k_dist, gas_concs, p_lay, t_lay, p_lev, & ! Gas optical depth -- pressure need to be expressed as Pa ! error_msg = k_dist%gas_optics(p_lay, p_lev, t_lay, gas_concs, & - optical_props, toa_flux, & + optical_props, toa_flux, & col_dry) if (error_msg /= '') return ! @@ -289,9 +276,8 @@ function rte_sw(k_dist, gas_concs, p_lay, t_lay, p_lev, & if(present(aer_props)) error_msg = aer_props%increment(optical_props) if(error_msg /= '') return - error_msg = base_rte_sw(optical_props, top_at_1, & - mu0, toa_flux, & - sfc_alb_dir, sfc_alb_dif, & + error_msg = base_rte_sw(optical_props, mu0, toa_flux, & + sfc_alb_dir, sfc_alb_dif, & clrsky_fluxes) if(error_msg /= '') return @@ -301,9 +287,8 @@ function rte_sw(k_dist, gas_concs, p_lay, t_lay, p_lev, & error_msg = cloud_props%increment(optical_props) if(error_msg /= '') return - error_msg = base_rte_sw(optical_props, top_at_1, & - mu0, toa_flux, & - sfc_alb_dir, sfc_alb_dif, & + error_msg = base_rte_sw(optical_props, mu0, toa_flux, & + sfc_alb_dir, sfc_alb_dif, & allsky_fluxes) call optical_props%finalize() diff --git a/rrtmgp-frontend/CMakeLists.txt b/rrtmgp-frontend/CMakeLists.txt new file mode 100644 index 000000000..386a8756b --- /dev/null +++ b/rrtmgp-frontend/CMakeLists.txt @@ -0,0 +1,27 @@ +set(gas_optics_source_dir ${PROJECT_SOURCE_DIR}/gas-optics) + +add_library( + rrtmgp STATIC # cmake-format: sort + $ + ${gas_optics_source_dir}/mo_gas_concentrations.F90 + ${gas_optics_source_dir}/mo_gas_optics.F90 + ${gas_optics_source_dir}/mo_gas_optics_constants.F90 + ${gas_optics_source_dir}/mo_gas_optics_util_string.F90 + mo_aerosol_optics_rrtmgp_merra.F90 + mo_cloud_optics_rrtmgp.F90 + mo_gas_optics_rrtmgp.F90 +) + +add_library(rte-rrtmgp::rrtmgp ALIAS rrtmgp) + +set_target_properties(rrtmgp PROPERTIES EXPORT_NAME rte-rrtmgp::rrtmgp) + +target_include_directories( + rrtmgp + PUBLIC + $:${CMAKE_Fortran_MODULE_DIRECTORY}>> +) + +target_link_libraries(rrtmgp PRIVATE rte) + +install(TARGETS rrtmgp EXPORT rte-rrtmgp-targets) diff --git a/rrtmgp-frontend/Make.depends b/rrtmgp-frontend/Make.depends deleted file mode 100644 index 1a06c30f5..000000000 --- a/rrtmgp-frontend/Make.depends +++ /dev/null @@ -1,40 +0,0 @@ -RRTMGP_FORTRAN_INTERFACE = \ - mo_gas_optics_util_string.o \ - mo_gas_optics_constants.o \ - mo_gas_concentrations.o \ - mo_gas_optics.o \ - mo_gas_optics_rrtmgp.o - -##### -# RRTMGP: RRTM for GCM Applications - Parallel -# Built on top of RTE, requiring mo_rte_kind.o, mo_rte_util_array.o, mo_rte_util_array_validation.o, mo_optical_props.o -# -# Physical constants -# -mo_gas_optics_constants.o: $(RTE_FORTRAN_INTERFACE) mo_gas_optics_constants.F90 -# -# Utility -# -mo_gas_optics_util_string.o: mo_gas_optics_util_string.F90 - -# -# Gas concentrations - used by gas optics base class -# -mo_gas_concentrations.o: $(RTE_FORTRAN_INTERFACE) mo_gas_concentrations.F90 - -# -# Gas optics base class -# -mo_gas_optics.o: \ - $(RTE_FORTRAN_INTERFACE) mo_gas_concentrations.o \ - mo_gas_optics.F90 - -# -# RRTMGP gas optics -# -mo_gas_optics_rrtmgp.o: \ - $(RTE_FORTRAN_INTERFACE) \ - mo_gas_optics_constants.o mo_gas_optics_util_string.o \ - mo_gas_concentrations.o \ - mo_gas_optics.o \ - mo_gas_optics_rrtmgp_kernels.o mo_gas_optics_rrtmgp.F90 diff --git a/rrtmgp-frontend/mo_cloud_optics_rrtmgp.F90 b/rrtmgp-frontend/mo_cloud_optics_rrtmgp.F90 index 582fe36e0..07102b534 100644 --- a/rrtmgp-frontend/mo_cloud_optics_rrtmgp.F90 +++ b/rrtmgp-frontend/mo_cloud_optics_rrtmgp.F90 @@ -28,10 +28,9 @@ module mo_cloud_optics_rrtmgp ty_optical_props_1scl, & ty_optical_props_2str, & ty_optical_props_nstr + use mo_cloud_optics_rrtmgp_kernels, only: & + compute_cld_from_table, compute_cld_from_pade implicit none - interface pade_eval - module procedure pade_eval_nbnd, pade_eval_1 - end interface pade_eval private ! ----------------------------------------------------------------------------------- type, extends(ty_optical_props), public :: ty_cloud_optics_rrtmgp @@ -486,14 +485,14 @@ function cloud_optics(this, & ! ! Liquid ! - call compute_all_from_table(ncol, nlay, nbnd, liqmsk, clwp, reliq, & + call compute_cld_from_table(ncol, nlay, nbnd, liqmsk, clwp, reliq, & this%liq_nsteps,this%liq_step_size,this%radliq_lwr, & this%lut_extliq, this%lut_ssaliq, this%lut_asyliq, & ltau, ltaussa, ltaussag) ! ! Ice ! - call compute_all_from_table(ncol, nlay, nbnd, icemsk, ciwp, reice, & + call compute_cld_from_table(ncol, nlay, nbnd, icemsk, ciwp, reice, & this%ice_nsteps,this%ice_step_size,this%radice_lwr, & this%lut_extice(:,:,this%icergh), & this%lut_ssaice(:,:,this%icergh), & @@ -505,13 +504,13 @@ function cloud_optics(this, & ! Hard coded assumptions: order of approximants, three size regimes ! nsizereg = size(this%pade_extliq,2) - call compute_all_from_pade(ncol, nlay, nbnd, nsizereg, & + call compute_cld_from_pade(ncol, nlay, nbnd, nsizereg, & liqmsk, clwp, reliq, & 2, 3, this%pade_sizreg_extliq, this%pade_extliq, & 2, 2, this%pade_sizreg_ssaliq, this%pade_ssaliq, & 2, 2, this%pade_sizreg_asyliq, this%pade_asyliq, & ltau, ltaussa, ltaussag) - call compute_all_from_pade(ncol, nlay, nbnd, nsizereg, & + call compute_cld_from_pade(ncol, nlay, nbnd, nsizereg, & icemsk, ciwp, reice, & 2, 3, this%pade_sizreg_extice, this%pade_extice(:,:,:,this%icergh), & 2, 2, this%pade_sizreg_ssaice, this%pade_ssaice(:,:,:,this%icergh), & @@ -619,185 +618,4 @@ function get_max_radius_ice(this) result(r) r = this%radice_upr end function get_max_radius_ice - !-------------------------------------------------------------------------------------------------------------------- - ! - ! Ancillary functions - ! - !-------------------------------------------------------------------------------------------------------------------- - ! - ! Linearly interpolate values from a lookup table with "nsteps" evenly-spaced - ! elements starting at "offset." The table's second dimension is band. - ! Returns 0 where the mask is false. - ! We could also try gather/scatter for efficiency - ! - subroutine compute_all_from_table(ncol, nlay, nbnd, mask, lwp, re, & - nsteps, step_size, offset, & - tau_table, ssa_table, asy_table, & - tau, taussa, taussag) - integer, intent(in) :: ncol, nlay, nbnd, nsteps - logical(wl), dimension(ncol,nlay), intent(in) :: mask - real(wp), dimension(ncol,nlay), intent(in) :: lwp, re - real(wp), intent(in) :: step_size, offset - real(wp), dimension(nsteps, nbnd), intent(in) :: tau_table, ssa_table, asy_table - real(wp), dimension(ncol,nlay,nbnd) :: tau, taussa, taussag - ! --------------------------- - integer :: icol, ilay, ibnd - integer :: index - real(wp) :: fint - real(wp) :: t, ts ! tau, tau*ssa, tau*ssa*g - ! --------------------------- - !$acc parallel loop gang vector default(present) collapse(3) - !$omp target teams distribute parallel do simd collapse(3) - do ibnd = 1, nbnd - do ilay = 1,nlay - do icol = 1, ncol - if(mask(icol,ilay)) then - index = min(floor((re(icol,ilay) - offset)/step_size)+1, nsteps-1) - fint = (re(icol,ilay) - offset)/step_size - (index-1) - t = lwp(icol,ilay) * & - (tau_table(index, ibnd) + fint * (tau_table(index+1,ibnd) - tau_table(index,ibnd))) - ts = t * & - (ssa_table(index, ibnd) + fint * (ssa_table(index+1,ibnd) - ssa_table(index,ibnd))) - taussag(icol,ilay,ibnd) = & - ts * & - (asy_table(index, ibnd) + fint * (asy_table(index+1,ibnd) - asy_table(index,ibnd))) - taussa (icol,ilay,ibnd) = ts - tau (icol,ilay,ibnd) = t - else - tau (icol,ilay,ibnd) = 0._wp - taussa (icol,ilay,ibnd) = 0._wp - taussag(icol,ilay,ibnd) = 0._wp - end if - end do - end do - end do - end subroutine compute_all_from_table - ! - ! Pade functions - ! - !--------------------------------------------------------------------------- - subroutine compute_all_from_pade(ncol, nlay, nbnd, nsizes, & - mask, lwp, re, & - m_ext, n_ext, re_bounds_ext, coeffs_ext, & - m_ssa, n_ssa, re_bounds_ssa, coeffs_ssa, & - m_asy, n_asy, re_bounds_asy, coeffs_asy, & - tau, taussa, taussag) - integer, intent(in) :: ncol, nlay, nbnd, nsizes - logical(wl), & - dimension(ncol,nlay), intent(in) :: mask - real(wp), dimension(ncol,nlay), intent(in) :: lwp, re - real(wp), dimension(nsizes+1), intent(in) :: re_bounds_ext, re_bounds_ssa, re_bounds_asy - integer, intent(in) :: m_ext, n_ext - real(wp), dimension(nbnd,nsizes,0:m_ext+n_ext), & - intent(in) :: coeffs_ext - integer, intent(in) :: m_ssa, n_ssa - real(wp), dimension(nbnd,nsizes,0:m_ssa+n_ssa), & - intent(in) :: coeffs_ssa - integer, intent(in) :: m_asy, n_asy - real(wp), dimension(nbnd,nsizes,0:m_asy+n_asy), & - intent(in) :: coeffs_asy - real(wp), dimension(ncol,nlay,nbnd) :: tau, taussa, taussag - ! --------------------------- - integer :: icol, ilay, ibnd, irad - real(wp) :: t, ts - - !$acc parallel loop gang vector default(present) collapse(3) - !$omp target teams distribute parallel do simd collapse(3) - do ibnd = 1, nbnd - do ilay = 1, nlay - do icol = 1, ncol - if(mask(icol,ilay)) then - ! - ! Finds index into size regime table - ! This works only if there are precisely three size regimes (four bounds) and it's - ! previously guaranteed that size_bounds(1) <= size <= size_bounds(4) - ! - irad = min(floor((re(icol,ilay) - re_bounds_ext(2))/re_bounds_ext(3))+2, 3) - t = lwp(icol,ilay) * & - pade_eval(ibnd, nbnd, nsizes, m_ext, n_ext, irad, re(icol,ilay), coeffs_ext) - - irad = min(floor((re(icol,ilay) - re_bounds_ssa(2))/re_bounds_ssa(3))+2, 3) - ! Pade approximants for co-albedo can sometimes be negative - ts = t * (1._wp - max(0._wp, & - pade_eval(ibnd, nbnd, nsizes, m_ssa, n_ssa, irad, re(icol,ilay), coeffs_ssa))) - - irad = min(floor((re(icol,ilay) - re_bounds_asy(2))/re_bounds_asy(3))+2, 3) - taussag(icol,ilay,ibnd) = & - ts * & - pade_eval(ibnd, nbnd, nsizes, m_asy, n_asy, irad, re(icol,ilay), coeffs_asy) - - taussa (icol,ilay,ibnd) = ts - tau (icol,ilay,ibnd) = t - else - tau (icol,ilay,ibnd) = 0._wp - taussa (icol,ilay,ibnd) = 0._wp - taussag(icol,ilay,ibnd) = 0._wp - end if - end do - end do - end do - - end subroutine compute_all_from_pade - !--------------------------------------------------------------------------- - ! - ! Evaluate Pade approximant of order [m/n] - ! - function pade_eval_nbnd(nbnd, nrads, m, n, irad, re, pade_coeffs) - integer, intent(in) :: nbnd, nrads, m, n, irad - real(wp), dimension(nbnd, nrads, 0:m+n), & - intent(in) :: pade_coeffs - real(wp), intent(in) :: re - real(wp), dimension(nbnd) :: pade_eval_nbnd - - integer :: iband - real(wp) :: numer, denom - integer :: i - - do iband = 1, nbnd - denom = pade_coeffs(iband,irad,n+m) - do i = n-1+m, 1+m, -1 - denom = pade_coeffs(iband,irad,i)+re*denom - end do - denom = 1._wp +re*denom - - numer = pade_coeffs(iband,irad,m) - do i = m-1, 1, -1 - numer = pade_coeffs(iband,irad,i)+re*numer - end do - numer = pade_coeffs(iband,irad,0) +re*numer - - pade_eval_nbnd(iband) = numer/denom - end do - end function pade_eval_nbnd - !--------------------------------------------------------------------------- - ! - ! Evaluate Pade approximant of order [m/n] - ! - function pade_eval_1(iband, nbnd, nrads, m, n, irad, re, pade_coeffs) - !$acc routine seq - !$omp declare target - ! - integer, intent(in) :: iband, nbnd, nrads, m, n, irad - real(wp), dimension(nbnd, nrads, 0:m+n), & - intent(in) :: pade_coeffs - real(wp), intent(in) :: re - real(wp) :: pade_eval_1 - - real(wp) :: numer, denom - integer :: i - - denom = pade_coeffs(iband,irad,n+m) - do i = n-1+m, 1+m, -1 - denom = pade_coeffs(iband,irad,i)+re*denom - end do - denom = 1._wp +re*denom - - numer = pade_coeffs(iband,irad,m) - do i = m-1, 1, -1 - numer = pade_coeffs(iband,irad,i)+re*numer - end do - numer = pade_coeffs(iband,irad,0) +re*numer - - pade_eval_1 = numer/denom - end function pade_eval_1 end module mo_cloud_optics_rrtmgp diff --git a/rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 b/rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 index 906c7cccb..967d656a9 100644 --- a/rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 +++ b/rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 @@ -258,6 +258,10 @@ function gas_optics_int(this, & ngpt = this%get_ngpt() nband = this%get_nband() ! + ! Vertical orientation + ! + call optical_props%set_top_at_1(play(1,1) < play(1, nlay)) + ! ! Gas optics ! !$acc enter data create(jtemp, jpress, tropo, fmajor, jeta) @@ -309,7 +313,7 @@ function gas_optics_int(this, & ! if(present(tlev)) then error_msg = source(this, & - ncol, nlay, nband, ngpt, & + ncol, nlay, nband, ngpt, optical_props%top_is_at_1(), & play, plev, tlay, tsfc, & jtemp, jpress, jeta, tropo, fmajor, & sources, & @@ -318,7 +322,7 @@ function gas_optics_int(this, & !$omp target exit data map(release:tlev) else error_msg = source(this, & - ncol, nlay, nband, ngpt, & + ncol, nlay, nband, ngpt, optical_props%top_is_at_1(), & play, plev, tlay, tsfc, & jtemp, jpress, jeta, tropo, fmajor, & sources) @@ -328,6 +332,7 @@ function gas_optics_int(this, & !$omp target exit data map(release:tsfc) !$acc exit data delete(jtemp, jpress, tropo, fmajor, jeta) !$omp target exit data map(release:jtemp, jpress, tropo, fmajor, jeta) + end function gas_optics_int !------------------------------------------------------------------------------------------ ! @@ -371,6 +376,10 @@ function gas_optics_ext(this, & ngas = this%get_ngas() nflav = get_nflav(this) ! + ! Vertical orientation + ! + call optical_props%set_top_at_1(play(1,1) < play(1, nlay)) + ! ! Gas optics ! !$acc enter data create(jtemp, jpress, tropo, fmajor, jeta) @@ -383,6 +392,7 @@ function gas_optics_ext(this, & col_dry) !$acc exit data delete(jtemp, jpress, tropo, fmajor, jeta) !$omp target exit data map(release:jtemp, jpress, tropo, fmajor, jeta) + if(error_msg /= '') return ! ---------------------------------------------------------- @@ -833,7 +843,7 @@ end function set_tsi ! Compute Planck source functions at layer centers and levels ! function source(this, & - ncol, nlay, nbnd, ngpt, & + ncol, nlay, nbnd, ngpt, top_at_1, & play, plev, tlay, tsfc, & jtemp, jpress, jeta, tropo, fmajor, & sources, & ! Planck sources @@ -842,6 +852,7 @@ function source(this, & ! inputs class(ty_gas_optics_rrtmgp), intent(in ) :: this integer, intent(in ) :: ncol, nlay, nbnd, ngpt + logical, intent(in ) :: top_at_1 real(wp), dimension(ncol,nlay), intent(in ) :: play ! layer pressures [Pa, mb] real(wp), dimension(ncol,nlay+1), intent(in ) :: plev ! level pressures [Pa, mb] real(wp), dimension(ncol,nlay), intent(in ) :: tlay ! layer temperatures [K] @@ -858,7 +869,6 @@ function source(this, & optional, target :: tlev ! level temperatures [K] character(len=128) :: error_msg ! ---------------------------------------------------------- - logical(wl) :: top_at_1 integer :: icol, ilay ! Variables for temperature at layer edges [K] (ncol, nlay+1) real(wp), dimension( ncol,nlay+1), target :: tlev_arr @@ -910,18 +920,12 @@ function source(this, & ! Compute internal (Planck) source functions at layers and levels, ! which depend on mapping from spectral space that creates k-distribution. - !$acc kernels copyout(top_at_1) - !$omp target map(from:top_at_1) - top_at_1 = play(1,1) < play(1, nlay) - !$acc end kernels - !$omp end target - call compute_Planck_source(ncol, nlay, nbnd, ngpt, & get_nflav(this), this%get_neta(), this%get_npres(), this%get_ntemp(), this%get_nPlanckTemp(), & - tlay, tlev_wk, tsfc, merge(nlay, 1, top_at_1), & - fmajor, jeta, tropo, jtemp, jpress, & + tlay, tlev_wk, tsfc, merge(nlay, 1, logical(top_at_1, wl)), & + fmajor, jeta, tropo, jtemp, jpress, & this%get_gpoint_bands(), this%get_band_lims_gpoint(), this%planck_frac, this%temp_ref_min,& - this%totplnk_delta, this%totplnk, this%gpoint_flavor, & + this%totplnk_delta, this%totplnk, this%gpoint_flavor, & sources%sfc_source, sources%lay_source, sources%lev_source, & sources%sfc_source_Jac) !$acc end data diff --git a/rrtmgp-kernels/CMakeLists.txt b/rrtmgp-kernels/CMakeLists.txt new file mode 100644 index 000000000..2c50a17c9 --- /dev/null +++ b/rrtmgp-kernels/CMakeLists.txt @@ -0,0 +1,36 @@ +add_library(rrtmgpkernels OBJECT) + +if(KERNEL_MODE STREQUAL "extern") + target_sources( + rrtmgpkernels + PRIVATE # cmake-format: sort + api/mo_cloud_optics_rrtmgp_kernels.F90 + api/mo_gas_optics_rrtmgp_kernels.F90 + api/rrtmgp_kernels.h + ) +else() + target_sources( + rrtmgpkernels + PRIVATE # cmake-format: sort + mo_cloud_optics_rrtmgp_kernels.F90 + ) + if(KERNEL_MODE STREQUAL "accel") + target_sources( + rrtmgpkernels + PRIVATE # cmake-format: sort + accel/mo_gas_optics_rrtmgp_kernels.F90 + ) + else() + target_sources( + rrtmgpkernels + PRIVATE # cmake-format: sort + mo_gas_optics_rrtmgp_kernels.F90 + ) + endif() +endif() + +target_include_directories( + rrtmgpkernels PRIVATE ${CMAKE_Fortran_MODULE_DIRECTORY} +) + +target_link_libraries(rrtmgpkernels PRIVATE rte) diff --git a/rrtmgp-kernels/Make.depends b/rrtmgp-kernels/Make.depends deleted file mode 100644 index c12e33c2c..000000000 --- a/rrtmgp-kernels/Make.depends +++ /dev/null @@ -1,6 +0,0 @@ -RRTMGP_FORTRAN_KERNELS = mo_gas_optics_rrtmgp_kernels.o - -# -# Gas optics -# -mo_gas_optics_rrtmgp_kernels.o: $(RTE_FORTRAN_KERNELS) mo_gas_optics_rrtmgp_kernels.F90 diff --git a/rrtmgp-kernels/api/mo_cloud_optics_rrtmgp_kernels.F90 b/rrtmgp-kernels/api/mo_cloud_optics_rrtmgp_kernels.F90 new file mode 100644 index 000000000..ed6959fa7 --- /dev/null +++ b/rrtmgp-kernels/api/mo_cloud_optics_rrtmgp_kernels.F90 @@ -0,0 +1,58 @@ +module mo_cloud_optics_rrtmgp_kernels + use mo_rte_kind, only : wp, wl + implicit none + private + public :: compute_cld_from_table, compute_cld_from_pade + interface + !--------------------------------------------------------------------------- + ! + ! Linearly interpolate values from a lookup table with "nsteps" evenly-spaced + ! elements starting at "offset." The table's second dimension is band. + ! Returns 0 where the mask is false. + ! We could also try gather/scatter for efficiency + ! + subroutine compute_cld_from_table(ncol, nlay, nbnd, mask, lwp, re, & + nsteps, step_size, offset, & + tau_table, ssa_table, asy_table, & + tau, taussa, taussag) bind(C, name="rrtmgp_compute_cld_from_table") + use mo_rte_kind, only : wp, wl + integer, intent(in) :: ncol, nlay, nbnd, nsteps + logical(wl), dimension(ncol,nlay), intent(in) :: mask + real(wp), dimension(ncol,nlay), intent(in) :: lwp, re + real(wp), intent(in) :: step_size, offset + real(wp), dimension(nsteps, nbnd), intent(in) :: tau_table, ssa_table, asy_table + real(wp), dimension(ncol,nlay,nbnd) :: tau, taussa, taussag + end subroutine compute_cld_from_table + + !--------------------------------------------------------------------------- + ! + ! Pade functions + ! + !--------------------------------------------------------------------------- + subroutine compute_cld_from_pade(ncol, nlay, nbnd, nsizes, & + mask, lwp, re, & + m_ext, n_ext, re_bounds_ext, coeffs_ext, & + m_ssa, n_ssa, re_bounds_ssa, coeffs_ssa, & + m_asy, n_asy, re_bounds_asy, coeffs_asy, & + tau, taussa, taussag) bind(C, name="rrtmgp_compute_cld_from_pade") + use mo_rte_kind, only : wp, wl + integer, intent(in) :: ncol, nlay, nbnd, nsizes + logical(wl), & + dimension(ncol,nlay), intent(in) :: mask + real(wp), dimension(ncol,nlay), intent(in) :: lwp, re + real(wp), dimension(nsizes+1), intent(in) :: re_bounds_ext, re_bounds_ssa, re_bounds_asy + integer, intent(in) :: m_ext, n_ext + real(wp), dimension(nbnd,nsizes,0:m_ext+n_ext), & + intent(in) :: coeffs_ext + integer, intent(in) :: m_ssa, n_ssa + real(wp), dimension(nbnd,nsizes,0:m_ssa+n_ssa), & + intent(in) :: coeffs_ssa + integer, intent(in) :: m_asy, n_asy + real(wp), dimension(nbnd,nsizes,0:m_asy+n_asy), & + intent(in) :: coeffs_asy + real(wp), dimension(ncol,nlay,nbnd) :: tau, taussa, taussag + end subroutine compute_cld_from_pade + end interface + !--------------------------------------------------------------------------- +end module mo_cloud_optics_rrtmgp_kernels + diff --git a/rrtmgp-kernels/api/rrtmgp_kernels.h b/rrtmgp-kernels/api/rrtmgp_kernels.h index 92c36957c..6d8976c3b 100644 --- a/rrtmgp-kernels/api/rrtmgp_kernels.h +++ b/rrtmgp-kernels/api/rrtmgp_kernels.h @@ -19,6 +19,7 @@ This header files defines the C bindings for the kernels used in RRTMGP extern "C" { + /* Gas optics kernels */ void rrtmgp_interpolation( const int& ncol, const int& nlay, const int& ngas, const int& nflav, const int& neta, @@ -120,4 +121,55 @@ extern "C" Float* lev_src, // [out] (ncol,nlay+1,ngpt) Float* sfc_src_jac // [out] (ncol, ngpt) ); + + /* Cloud optics kernels */ + void rrtmgp_compute_tau_rayleigh( + const int& ncol, const int& nlay, const int& nband, const int& ngpt, + const int& ngas, const int& nflav, const int& neta, const int& npres, const int& ntemp, + const int* gpoint_flavor, // (2,ngpt) + const int* band_lims_gpt, // (2,nbnd) + const Float* krayl, // (ntemp,neta,ngpt,2) + const int& idx_h2o, + const Float* col_dry, // (ncol,nlay) + const Float* col_gas, // (ncol,nlay,ngas+1) + const Float* fminor, // (2,2,ncol,nlay,nflav) + const int* jeta, // (2, ncol,nlay,nflav) + const Bool* tropo, // (ncol,nlay) + const int* jtemp, // (ncol,nlay) + Float* tau_rayleigh // [inout] (ncol,nlay.ngpt) + ); + + void rrtmgp_compute_cld_from_table( + const int& ncol, int& nlay, int& nbnd, int& nsteps, + const Bool* mask, // (ncol,nlay) + const Float* lwp, // (ncol,nlay) + const Float* re, // (ncol,nlay) + const Float& step_size, + const Float& offset, + const Float* tau_table, // (nsteps, nbnd) + const Float* ssa_table, // (nsteps, nbnd) + const Float* asy_table, // (nsteps, nbnd) + Float* tau, // (ncol,nlay,nbnd) + Float* taussa, // (ncol,nlay,nbnd) + Float* taussag // (ncol,nlay,nbnd) + ); + + void rrtmgp_compute_cld_from_pade( + const int& ncol, int& nlay, int& nbnd, int& nsizes, + const Bool* mask, // (ncol,nlay) + const Float* lwp, // (ncol,nlay) + const Float* re, // (ncol,nlay) + const Float* re_bounds_ext, // (nsizes+1) + const Float* re_bounds_ssa, // (nsizes+1) + const Float* re_bounds_asy, // (nsizes+1) + const int& m_ext, int& n_ext, + const Float* coeffs_ext, // (nbnd,nsizes,0:m_ext+n_ext) + const int& m_ssa, int& n_ssa, + const Float* coeffs_ssa, // (nbnd,nsizes,0:m_ssa+n_ssa) + const int& m_asy, int& n_asy, + const Float* coeffs_asy, // (nbnd,nsizes,0:m_asy+n_asy) + Float* tau, // (ncol,nlay,nbnd) + Float* taussa, // (ncol,nlay,nbnd) + Float* taussag // (ncol,nlay,nbnd) + ); } diff --git a/rrtmgp-kernels/mo_cloud_optics_rrtmgp_kernels.F90 b/rrtmgp-kernels/mo_cloud_optics_rrtmgp_kernels.F90 new file mode 100644 index 000000000..a5c9e295e --- /dev/null +++ b/rrtmgp-kernels/mo_cloud_optics_rrtmgp_kernels.F90 @@ -0,0 +1,198 @@ +! This code is part of +! RRTM for GCM Applications - Parallel (RRTMGP) +! +! Copyright 2024-, Atmospheric and Environmental Research, +! Trustees of Columbia University. 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 +! +module mo_cloud_optics_rrtmgp_kernels + use mo_rte_kind, only : wp, wl + implicit none + private + public :: compute_cld_from_table, compute_cld_from_pade + interface pade_eval + module procedure pade_eval_nbnd, pade_eval_1 + end interface pade_eval +contains + !--------------------------------------------------------------------------- + ! + ! Linearly interpolate values from a lookup table with "nsteps" evenly-spaced + ! elements starting at "offset." The table's second dimension is band. + ! Returns 0 where the mask is false. + ! We could also try gather/scatter for efficiency + ! + subroutine compute_cld_from_table(ncol, nlay, nbnd, mask, lwp, re, & + nsteps, step_size, offset, & + tau_table, ssa_table, asy_table, & + tau, taussa, taussag) bind(C, name="rrtmgp_compute_cld_from_table") + integer, intent(in) :: ncol, nlay, nbnd, nsteps + logical(wl), dimension(ncol,nlay), intent(in) :: mask + real(wp), dimension(ncol,nlay), intent(in) :: lwp, re + real(wp), intent(in) :: step_size, offset + real(wp), dimension(nsteps, nbnd), intent(in) :: tau_table, ssa_table, asy_table + real(wp), dimension(ncol,nlay,nbnd) :: tau, taussa, taussag + ! --------------------------- + integer :: icol, ilay, ibnd + integer :: index + real(wp) :: fint + real(wp) :: t, ts ! tau, tau*ssa, tau*ssa*g + ! --------------------------- + !$acc parallel loop gang vector default(present) collapse(3) + !$omp target teams distribute parallel do simd collapse(3) + do ibnd = 1, nbnd + do ilay = 1,nlay + do icol = 1, ncol + if(mask(icol,ilay)) then + index = min(floor((re(icol,ilay) - offset)/step_size)+1, nsteps-1) + fint = (re(icol,ilay) - offset)/step_size - (index-1) + t = lwp(icol,ilay) * & + (tau_table(index, ibnd) + fint * (tau_table(index+1,ibnd) - tau_table(index,ibnd))) + ts = t * & + (ssa_table(index, ibnd) + fint * (ssa_table(index+1,ibnd) - ssa_table(index,ibnd))) + taussag(icol,ilay,ibnd) = & + ts * & + (asy_table(index, ibnd) + fint * (asy_table(index+1,ibnd) - asy_table(index,ibnd))) + taussa (icol,ilay,ibnd) = ts + tau (icol,ilay,ibnd) = t + else + tau (icol,ilay,ibnd) = 0._wp + taussa (icol,ilay,ibnd) = 0._wp + taussag(icol,ilay,ibnd) = 0._wp + end if + end do + end do + end do + end subroutine compute_cld_from_table + !--------------------------------------------------------------------------- + ! + ! Pade functions + ! + !--------------------------------------------------------------------------- + subroutine compute_cld_from_pade(ncol, nlay, nbnd, nsizes, & + mask, lwp, re, & + m_ext, n_ext, re_bounds_ext, coeffs_ext, & + m_ssa, n_ssa, re_bounds_ssa, coeffs_ssa, & + m_asy, n_asy, re_bounds_asy, coeffs_asy, & + tau, taussa, taussag) bind(C, name="rrtmgp_compute_cld_from_pade") + integer, intent(in) :: ncol, nlay, nbnd, nsizes + logical(wl), & + dimension(ncol,nlay), intent(in) :: mask + real(wp), dimension(ncol,nlay), intent(in) :: lwp, re + real(wp), dimension(nsizes+1), intent(in) :: re_bounds_ext, re_bounds_ssa, re_bounds_asy + integer, intent(in) :: m_ext, n_ext + real(wp), dimension(nbnd,nsizes,0:m_ext+n_ext), & + intent(in) :: coeffs_ext + integer, intent(in) :: m_ssa, n_ssa + real(wp), dimension(nbnd,nsizes,0:m_ssa+n_ssa), & + intent(in) :: coeffs_ssa + integer, intent(in) :: m_asy, n_asy + real(wp), dimension(nbnd,nsizes,0:m_asy+n_asy), & + intent(in) :: coeffs_asy + real(wp), dimension(ncol,nlay,nbnd) :: tau, taussa, taussag + ! --------------------------- + integer :: icol, ilay, ibnd, irad + real(wp) :: t, ts + + !$acc parallel loop gang vector default(present) collapse(3) + !$omp target teams distribute parallel do simd collapse(3) + do ibnd = 1, nbnd + do ilay = 1, nlay + do icol = 1, ncol + if(mask(icol,ilay)) then + ! + ! Finds index into size regime table + ! This works only if there are precisely three size regimes (four bounds) and it's + ! previously guaranteed that size_bounds(1) <= size <= size_bounds(4) + ! + irad = min(floor((re(icol,ilay) - re_bounds_ext(2))/re_bounds_ext(3))+2, 3) + t = lwp(icol,ilay) * & + pade_eval(ibnd, nbnd, nsizes, m_ext, n_ext, irad, re(icol,ilay), coeffs_ext) + + irad = min(floor((re(icol,ilay) - re_bounds_ssa(2))/re_bounds_ssa(3))+2, 3) + ! Pade approximants for co-albedo can sometimes be negative + ts = t * (1._wp - max(0._wp, & + pade_eval(ibnd, nbnd, nsizes, m_ssa, n_ssa, irad, re(icol,ilay), coeffs_ssa))) + + irad = min(floor((re(icol,ilay) - re_bounds_asy(2))/re_bounds_asy(3))+2, 3) + taussag(icol,ilay,ibnd) = & + ts * & + pade_eval(ibnd, nbnd, nsizes, m_asy, n_asy, irad, re(icol,ilay), coeffs_asy) + + taussa (icol,ilay,ibnd) = ts + tau (icol,ilay,ibnd) = t + else + tau (icol,ilay,ibnd) = 0._wp + taussa (icol,ilay,ibnd) = 0._wp + taussag(icol,ilay,ibnd) = 0._wp + end if + end do + end do + end do + + end subroutine compute_cld_from_pade + !--------------------------------------------------------------------------- + ! + ! Evaluate Pade approximant of order [m/n] + ! + function pade_eval_nbnd(nbnd, nrads, m, n, irad, re, pade_coeffs) + integer, intent(in) :: nbnd, nrads, m, n, irad + real(wp), dimension(nbnd, nrads, 0:m+n), & + intent(in) :: pade_coeffs + real(wp), intent(in) :: re + real(wp), dimension(nbnd) :: pade_eval_nbnd + + integer :: iband + real(wp) :: numer, denom + integer :: i + + do iband = 1, nbnd + denom = pade_coeffs(iband,irad,n+m) + do i = n-1+m, 1+m, -1 + denom = pade_coeffs(iband,irad,i)+re*denom + end do + denom = 1._wp +re*denom + + numer = pade_coeffs(iband,irad,m) + do i = m-1, 1, -1 + numer = pade_coeffs(iband,irad,i)+re*numer + end do + numer = pade_coeffs(iband,irad,0) +re*numer + + pade_eval_nbnd(iband) = numer/denom + end do + end function pade_eval_nbnd + !--------------------------------------------------------------------------- + ! + ! Evaluate Pade approximant of order [m/n] + ! + function pade_eval_1(iband, nbnd, nrads, m, n, irad, re, pade_coeffs) + !$acc routine seq + !$omp declare target + ! + integer, intent(in) :: iband, nbnd, nrads, m, n, irad + real(wp), dimension(nbnd, nrads, 0:m+n), & + intent(in) :: pade_coeffs + real(wp), intent(in) :: re + real(wp) :: pade_eval_1 + + real(wp) :: numer, denom + integer :: i + + denom = pade_coeffs(iband,irad,n+m) + do i = n-1+m, 1+m, -1 + denom = pade_coeffs(iband,irad,i)+re*denom + end do + denom = 1._wp +re*denom + + numer = pade_coeffs(iband,irad,m) + do i = m-1, 1, -1 + numer = pade_coeffs(iband,irad,i)+re*numer + end do + numer = pade_coeffs(iband,irad,0) +re*numer + + pade_eval_1 = numer/denom + end function pade_eval_1 + +end module mo_cloud_optics_rrtmgp_kernels diff --git a/rrtmgp-kernels/mo_gas_optics_rrtmgp_kernels.F90 b/rrtmgp-kernels/mo_gas_optics_rrtmgp_kernels.F90 index 7a050fb84..df656859b 100644 --- a/rrtmgp-kernels/mo_gas_optics_rrtmgp_kernels.F90 +++ b/rrtmgp-kernels/mo_gas_optics_rrtmgp_kernels.F90 @@ -71,22 +71,9 @@ subroutine interpolation( & logical(wl), dimension(ncol,nlay), intent(out) :: tropo !! use lower (or upper) atmosphere tables integer, dimension(2, ncol,nlay,nflav), intent(out) :: jeta - !! Index for binary species interpolation -#if !defined(__INTEL_LLVM_COMPILER) && __INTEL_COMPILER >= 1910 - ! A performance-hitting workaround for the vectorization problem reported in - ! https://github.com/earth-system-radiation/rte-rrtmgp/issues/159 - ! The known affected compilers are Intel Fortran Compiler Classic - ! 2021.4, 2021.5 and 2022.1. We do not limit the workaround to these - ! versions because it is not clear when the compiler bug will be fixed, see - ! https://community.intel.com/t5/Intel-Fortran-Compiler/Compiler-vectorization-bug/m-p/1362591. - ! We, however, limit the workaround to the Classic versions only since the - ! problem is not confirmed for the Intel Fortran Compiler oneAPI (a.k.a - ! 'ifx'), which does not mean there is none though. - real(wp), dimension(:, :, :, :), intent(out) :: col_mix -#else + !! Index for binary species interpolation real(wp), dimension(2, ncol,nlay,nflav), intent(out) :: col_mix !! combination of major species's column amounts (first index is strat/trop) -#endif real(wp), dimension(2,2,2,ncol,nlay,nflav), intent(out) :: fmajor !! Interpolation weights in pressure, eta, strat/trop real(wp), dimension(2,2, ncol,nlay,nflav), intent(out) :: fminor @@ -100,29 +87,40 @@ subroutine interpolation( & real(wp) :: eta, feta ! binary_species_parameter, interpolation variable for eta real(wp) :: loceta ! needed to find location in eta grid real(wp) :: ftemp_term + real(wp) :: press_ref_trop + real(wp) :: temp_ref_delta_inv + real(wp) :: press_ref_log_1 + real(wp) :: press_ref_log_delta_inv + real(wp) :: jpress_aint ! ----------------- ! local indexes - integer :: icol, ilay, iflav, igases(2), itropo, itemp + integer :: icol, ilay, iflav, igas_1, igas_2, itropo, itemp, jtemp_ + press_ref_trop = exp(press_ref_trop_log) + temp_ref_delta_inv = 1.0_wp / temp_ref_delta + press_ref_log_1 = press_ref_log(1) + press_ref_log_delta_inv = 1.0_wp / press_ref_log_delta do ilay = 1, nlay do icol = 1, ncol ! index and factor for temperature interpolation - jtemp(icol,ilay) = int((tlay(icol,ilay) - (temp_ref_min - temp_ref_delta)) / temp_ref_delta) - jtemp(icol,ilay) = min(ntemp - 1, max(1, jtemp(icol,ilay))) ! limit the index range - ftemp(icol,ilay) = (tlay(icol,ilay) - temp_ref(jtemp(icol,ilay))) / temp_ref_delta + jtemp_ = INT((tlay(icol,ilay) - (temp_ref_min - temp_ref_delta)) * temp_ref_delta_inv) + jtemp(icol,ilay) = min(ntemp - 1, max(1, jtemp_)) ! limit the index range + ftemp(icol,ilay) = (tlay(icol,ilay) - temp_ref(jtemp_)) * temp_ref_delta_inv ! index and factor for pressure interpolation - locpress = 1._wp + (log(play(icol,ilay)) - press_ref_log(1)) / press_ref_log_delta - jpress(icol,ilay) = min(npres-1, max(1, int(locpress))) - fpress(icol,ilay) = locpress - float(jpress(icol,ilay)) + locpress = 1._wp + (log(play(icol,ilay)) - press_ref_log_1) * press_ref_log_delta_inv + jpress_aint = min(real(npres-1, wp), max(1.0_wp, aint(locpress))) + jpress(icol,ilay) = int(jpress_aint) + fpress(icol,ilay) = locpress - jpress_aint ! determine if in lower or upper part of atmosphere - tropo(icol,ilay) = log(play(icol,ilay)) > press_ref_trop_log + tropo(icol,ilay) = play(icol,ilay) > press_ref_trop end do end do do iflav = 1, nflav - igases(:) = flavor(:,iflav) + igas_1 = flavor(1,iflav) + igas_2 = flavor(2,iflav) do ilay = 1, nlay do icol = 1, ncol ! itropo = 1 lower atmosphere; itropo = 2 upper atmosphere @@ -132,9 +130,9 @@ subroutine interpolation( & ! compute interpolation fractions needed for lower, then upper reference temperature level ! compute binary species parameter (eta) for flavor and temperature and ! associated interpolation index and factors - ratio_eta_half = vmr_ref(itropo,igases(1),(jtemp(icol,ilay)+itemp-1)) / & - vmr_ref(itropo,igases(2),(jtemp(icol,ilay)+itemp-1)) - col_mix(itemp,icol,ilay,iflav) = col_gas(icol,ilay,igases(1)) + ratio_eta_half * col_gas(icol,ilay,igases(2)) + ratio_eta_half = vmr_ref(itropo,igas_1,(jtemp(icol,ilay)+itemp-1)) / & + vmr_ref(itropo,igas_2,(jtemp(icol,ilay)+itemp-1)) + col_mix(itemp,icol,ilay,iflav) = col_gas(icol,ilay,igas_1) + ratio_eta_half * col_gas(icol,ilay,igas_2) ! Keep this commented lines. Fortran does allow for ! substantial optimizations and in this merge cases may ! happen that all expressions are evaluated and so create @@ -142,18 +140,18 @@ subroutine interpolation( & ! save. Merge is the way to do it in general inside of ! loops, but sometimes it may not work. ! - ! eta = merge(col_gas(icol,ilay,igases(1)) / col_mix(itemp,icol,ilay,iflav), 0.5_wp, & + ! eta = merge(col_gas(icol,ilay,igas_1) / col_mix(itemp,icol,ilay,iflav), 0.5_wp, & ! col_mix(itemp,icol,ilay,iflav) > 2._wp * tiny(col_mix)) ! ! In essence: do not turn it back to merge(...)! if (col_mix(itemp,icol,ilay,iflav) > 2._wp * tiny(col_mix)) then - eta = col_gas(icol,ilay,igases(1)) / col_mix(itemp,icol,ilay,iflav) + eta = col_gas(icol,ilay,igas_1) / col_mix(itemp,icol,ilay,iflav) else eta = 0.5_wp endif - loceta = eta * float(neta-1) + loceta = eta * real(neta-1, wp) jeta(itemp,icol,ilay,iflav) = min(int(loceta)+1, neta-1) - feta = mod(loceta, 1.0_wp) + feta = loceta - aint(loceta) ! compute interpolation fractions needed for minor species ! ftemp_term = (1._wp-ftemp(icol,ilay)) for itemp = 1, ftemp(icol,ilay) for itemp=2 ftemp_term = (real(2-itemp, wp) + real(2*itemp-3, wp) * ftemp(icol,ilay)) diff --git a/rte-frontend/CMakeLists.txt b/rte-frontend/CMakeLists.txt new file mode 100644 index 000000000..9a84452a0 --- /dev/null +++ b/rte-frontend/CMakeLists.txt @@ -0,0 +1,23 @@ +add_library( + rte STATIC # cmake-format: sort + $ + mo_fluxes.F90 + mo_optical_props.F90 + mo_rte_config.F90 + mo_rte_lw.F90 + mo_rte_sw.F90 + mo_rte_util_array_validation.F90 + mo_source_functions.F90 +) + +add_library(rte-rrtmgp::rte ALIAS rte) + +set_target_properties(rte PROPERTIES EXPORT_NAME rte-rrtmgp::rte) + +target_include_directories( + rte + PUBLIC + $:${CMAKE_Fortran_MODULE_DIRECTORY}>> +) + +install(TARGETS rte EXPORT rte-rrtmgp-targets) diff --git a/rte-frontend/Make.depends b/rte-frontend/Make.depends deleted file mode 100644 index 5b5b098a5..000000000 --- a/rte-frontend/Make.depends +++ /dev/null @@ -1,51 +0,0 @@ -RTE_FORTRAN_INTERFACE = \ - mo_rte_kind.o \ - mo_rte_config.o \ - mo_rte_util_array_validation.o \ - mo_optical_props.o \ - mo_source_functions.o \ - mo_fluxes.o \ - mo_rte_lw.o \ - mo_rte_sw.o - -################################## -# RTE - Radiative transfer for energetics -################################## -# -# -mo_rte_config.o: mo_rte_config.F90 mo_rte_kind.o -# -# -mo_rte_util_array_validation.o: mo_rte_util_array_validation.F90 mo_rte_kind.o -# -# Optical properties -# -mo_optical_props.o: mo_rte_kind.o mo_rte_util_array_validation.o mo_optical_props_kernels.o mo_optical_props.F90 -# -# Source functions -# -mo_source_functions.o: mo_rte_kind.o mo_optical_props.o mo_source_functions.F90 -# -# Flux reduction -# -mo_fluxes.o: mo_rte_kind.o mo_fluxes_broadband_kernels.o mo_rte_config.o mo_optical_props.o mo_rte_util_array_validation.o mo_fluxes.F90 - -mo_rte_lw.o: mo_rte_kind.o \ - mo_rte_config.o \ - mo_rte_util_array.o \ - mo_rte_util_array_validation.o \ - mo_optical_props.o \ - mo_source_functions.o \ - mo_fluxes.o \ - mo_rte_solver_kernels.o \ - mo_rte_lw.F90 - -mo_rte_sw.o: mo_rte_kind.o \ - mo_rte_config.o \ - mo_rte_util_array.o \ - mo_rte_util_array_validation.o \ - mo_optical_props.o \ - mo_source_functions.o \ - mo_fluxes.o \ - mo_rte_solver_kernels.o \ - mo_rte_sw.F90 diff --git a/rte-frontend/mo_optical_props.F90 b/rte-frontend/mo_optical_props.F90 index c49fd5dc8..5f216a52e 100644 --- a/rte-frontend/mo_optical_props.F90 +++ b/rte-frontend/mo_optical_props.F90 @@ -33,6 +33,7 @@ !> !> These classes must be allocated before use. Initialization and allocation can be combined. !> The classes have a validate() function that checks all arrays for valid values (e.g. tau > 0.) +!> The vertical orientation can be specified via this%set_top_at_1() or obtained via this%top_at_1(). !> !> Optical properties can be delta-scaled (though this is currently implemented only for two-stream arrays) !> @@ -51,7 +52,7 @@ !@endnote !> ------------------------------------------------------------------------------------------------- module mo_optical_props - use mo_rte_kind, only: wp + use mo_rte_kind, only: wp, wl use mo_rte_config, only: check_extents, check_values use mo_rte_util_array_validation, & only: any_vals_less_than, any_vals_outside, extents_are @@ -109,6 +110,7 @@ module mo_optical_props ! ------------------------------------------------------------------------------------------------- type, extends(ty_optical_props), abstract, public :: ty_optical_props_arry real(wp), dimension(:,:,:), allocatable :: tau !! optical depth (ncol, nlay, ngpt) + logical(wl), private :: top_at_1 ! No default - maybe uninitialized values will get caught? contains procedure, public :: get_ncol procedure, public :: get_nlay @@ -116,6 +118,11 @@ module mo_optical_props !> Increment another set of values !> procedure, public :: increment + !> + !> + !> + procedure, public :: top_is_at_1 + procedure, public :: set_top_at_1 !> !> Deferred procedures -- each must be implemented in each child class with @@ -757,6 +764,7 @@ function subset_1scl_range(full, start, n, subset) result(err_message) subset%p(:,1:n,:,:) = 0._wp end select call extract_subset(ncol, nlay, ngpt, full%tau, start, start+n-1, subset%tau) + call subset%set_top_at_1(full%top_is_at_1()) end function subset_1scl_range ! ------------------------------------------------------------------------------------------ @@ -810,6 +818,7 @@ function subset_2str_range(full, start, n, subset) result(err_message) subset%p(1,1:n,:,:) = full%g (start:start+n-1,:,:) subset%p(2:,:, :,:) = 0._wp end select + call subset%set_top_at_1(full%top_is_at_1()) end function subset_2str_range ! ------------------------------------------------------------------------------------------ @@ -858,6 +867,7 @@ function subset_nstr_range(full, start, n, subset) result(err_message) call extract_subset( ncol, nlay, ngpt, full%ssa, start, start+n-1, subset%ssa) call extract_subset(nmom, ncol, nlay, ngpt, full%p , start, start+n-1, subset%p ) end select + call subset%set_top_at_1(full%top_is_at_1()) end function subset_nstr_range !> ------------------------------------------------------------------------------------------ @@ -1057,6 +1067,30 @@ pure function get_nmom(this) get_nmom = 0 end if end function get_nmom + !> ----------------------------------------------------------------------------------------------- + !> + !> Routines for array classes: vertical orientation + !> + ! ------------------------------------------------------------------------------------------ + pure function top_is_at_1(this) + ! + ! Vertical orientation - .true. if array index 1 is top of atmosphere + ! + class(ty_optical_props_arry), intent(in) :: this + logical :: top_is_at_1 + + top_is_at_1 = this%top_at_1 + end function top_is_at_1 + ! ------------------------------------------------------------------------------------------ + subroutine set_top_at_1(this, top_at_1) + ! + !> Set vertical orientation of class - .true. if array index 1 is top of atmosphere + ! + class(ty_optical_props_arry), intent(inout) :: this + logical, intent(in ) :: top_at_1 + + this%top_at_1 = top_at_1 + end subroutine set_top_at_1 ! ----------------------------------------------------------------------------------------------- ! ! Routines for base class: spectral discretization diff --git a/rte-frontend/mo_rte_lw.F90 b/rte-frontend/mo_rte_lw.F90 index 47a11ebe6..d03eaf1e9 100644 --- a/rte-frontend/mo_rte_lw.F90 +++ b/rte-frontend/mo_rte_lw.F90 @@ -67,15 +67,13 @@ module mo_rte_lw ! Interface using only optical properties and source functions as inputs; fluxes as outputs. ! ! -------------------------------------------------- - function rte_lw(optical_props, top_at_1, & - sources, sfc_emis, & - fluxes, & + function rte_lw(optical_props, & + sources, sfc_emis, & + fluxes, & inc_flux, n_gauss_angles, use_2stream, & lw_Ds, flux_up_Jac) result(error_msg) class(ty_optical_props_arry), intent(in ) :: optical_props !! Set of optical properties as one or more arrays - logical, intent(in ) :: top_at_1 - !! Is the top of the domain at index 1? (if not, ordering is bottom-to-top) type(ty_source_func_lw), intent(in ) :: sources !! Derived type with Planck source functions real(wp), dimension(:,:), intent(in ) :: sfc_emis @@ -353,8 +351,8 @@ function rte_lw(optical_props, top_at_1, & end do end if call lw_solver_noscat(ncol, nlay, ngpt, & - logical(top_at_1, wl), n_quad_angs, & - secants, gauss_wts(1:n_quad_angs,n_quad_angs), & + logical(optical_props%top_is_at_1(), wl), & + n_quad_angs, secants, gauss_wts(1:n_quad_angs,n_quad_angs), & optical_props%tau, & sources%lay_source, & sources%lev_source, & @@ -373,7 +371,8 @@ function rte_lw(optical_props, top_at_1, & ! ! two-stream calculation with scattering ! - call lw_solver_2stream(ncol, nlay, ngpt, logical(top_at_1, wl), & + call lw_solver_2stream(ncol, nlay, ngpt, & + logical(optical_props%top_is_at_1(), wl), & optical_props%tau, optical_props%ssa, optical_props%g, & sources%lay_source, sources%lev_source, & sfc_emis_gpt, sources%sfc_source, & @@ -396,7 +395,8 @@ function rte_lw(optical_props, top_at_1, & ! Re-scaled solution to account for scattering ! call lw_solver_noscat(ncol, nlay, ngpt, & - logical(top_at_1, wl), n_quad_angs, & + logical(optical_props%top_is_at_1(), wl), & + n_quad_angs, & secants, gauss_wts(1:n_quad_angs,n_quad_angs), & optical_props%tau, & sources%lay_source, & @@ -438,7 +438,7 @@ function rte_lw(optical_props, top_at_1, & ! ! ...or reduce spectral fluxes to desired output quantities ! - error_msg = fluxes%reduce(gpt_flux_up, gpt_flux_dn, optical_props, top_at_1) + error_msg = fluxes%reduce(gpt_flux_up, gpt_flux_dn, optical_props, optical_props%top_is_at_1()) end select end if ! no error message from validation !$acc end data diff --git a/rte-frontend/mo_rte_sw.F90 b/rte-frontend/mo_rte_sw.F90 index 22cdf3aeb..7a4b359c0 100644 --- a/rte-frontend/mo_rte_sw.F90 +++ b/rte-frontend/mo_rte_sw.F90 @@ -53,14 +53,12 @@ module mo_rte_sw contains ! ------------------------------------------------------------------------------------------------- - function rte_sw_mu0_bycol(atmos, top_at_1, & + function rte_sw_mu0_bycol(atmos, & mu0, inc_flux, & sfc_alb_dir, sfc_alb_dif, & fluxes, inc_flux_dif) result(error_msg) class(ty_optical_props_arry), intent(in ) :: atmos !! Optical properties provided as arrays - logical, intent(in ) :: top_at_1 - !! Is the top of the domain at index 1? (if not, ordering is bottom-to-top) real(wp), dimension(:), intent(in ) :: mu0 !! cosine of solar zenith angle (ncol) - will be assumed constant with height real(wp), dimension(:,:), intent(in ) :: inc_flux @@ -94,22 +92,19 @@ function rte_sw_mu0_bycol(atmos, top_at_1, & end do end do - error_msg = rte_sw_mu0_full(atmos, top_at_1, & - mu0_bylay, inc_flux, & - sfc_alb_dir, sfc_alb_dif, & + error_msg = rte_sw_mu0_full(atmos, & + mu0_bylay, inc_flux, & + sfc_alb_dir, sfc_alb_dif, & fluxes, inc_flux_dif) !$acc end data !$omp end target data end function rte_sw_mu0_bycol ! ------------------------------------------------------------------------------------------------- - function rte_sw_mu0_full(atmos, top_at_1, & - mu0, inc_flux, & - sfc_alb_dir, sfc_alb_dif, & - fluxes, inc_flux_dif) result(error_msg) + function rte_sw_mu0_full(atmos, mu0, inc_flux, & + sfc_alb_dir, sfc_alb_dif, & + fluxes, inc_flux_dif) result(error_msg) class(ty_optical_props_arry), intent(in ) :: atmos !! Optical properties provided as arrays - logical, intent(in ) :: top_at_1 - !! Is the top of the domain at index 1? (if not, ordering is bottom-to-top) real(wp), dimension(:,:), intent(in ) :: mu0 !! cosine of solar zenith angle (ncol, nlay) real(wp), dimension(:,:), intent(in ) :: inc_flux @@ -293,8 +288,9 @@ function rte_sw_mu0_full(atmos, top_at_1, & ! ! Direct beam only - for completeness, unlikely to be used in practice ! - call sw_solver_noscat(ncol, nlay, ngpt, logical(top_at_1, wl), & - atmos%tau, mu0, inc_flux, & + call sw_solver_noscat(ncol, nlay, ngpt, & + logical(atmos%top_is_at_1(), wl), & + atmos%tau, mu0, inc_flux, & gpt_flux_dir) call zero_array(ncol, nlay+1, ngpt, gpt_flux_up) ! @@ -308,7 +304,8 @@ function rte_sw_mu0_full(atmos, top_at_1, & ! ! two-stream calculation with scattering ! - call sw_solver_2stream(ncol, nlay, ngpt, logical(top_at_1, wl), & + call sw_solver_2stream(ncol, nlay, ngpt, & + logical(atmos%top_is_at_1(), wl), & atmos%tau, atmos%ssa, atmos%g, mu0, & sfc_alb_dir_gpt, sfc_alb_dif_gpt, & inc_flux, & @@ -348,7 +345,7 @@ function rte_sw_mu0_full(atmos, top_at_1, & ! ! ...or reduce spectral fluxes to desired output quantities ! - error_msg = fluxes%reduce(gpt_flux_up, gpt_flux_dn, atmos, top_at_1, gpt_flux_dir) + error_msg = fluxes%reduce(gpt_flux_up, gpt_flux_dn, atmos, atmos%top_is_at_1(), gpt_flux_dir) end select end if ! In case of an error we exit here diff --git a/rte-kernels/CMakeLists.txt b/rte-kernels/CMakeLists.txt new file mode 100644 index 000000000..056feb047 --- /dev/null +++ b/rte-kernels/CMakeLists.txt @@ -0,0 +1,47 @@ +add_library(rtekernels OBJECT mo_rte_kind.F90) + +if(KERNEL_MODE STREQUAL "extern") + target_sources( + rtekernels + PRIVATE # cmake-format: sort + api/mo_fluxes_broadband_kernels.F90 + api/mo_optical_props_kernels.F90 + api/mo_rte_solver_kernels.F90 + api/mo_rte_util_array.F90 + api/rte_kernels.h + api/rte_types.h + ) +else() + target_sources( + rtekernels + PRIVATE # cmake-format: sort + mo_fluxes_broadband_kernels.F90 + mo_rte_util_array.F90 + ) + if(KERNEL_MODE STREQUAL "accel") + target_sources( + rtekernels + PRIVATE # cmake-format: sort + accel/mo_optical_props_kernels.F90 + accel/mo_rte_solver_kernels.F90 + ) + else() + target_sources( + rtekernels + PRIVATE # cmake-format: sort + mo_optical_props_kernels.F90 + mo_rte_solver_kernels.F90 + ) + endif() +endif() + +include(CheckFortranNeedsCBool) +check_fortran_needs_cbool(RTE_USE_C_BOOL) + +target_compile_definitions( + rtekernels + PRIVATE $<$:RTE_USE_SP> + $<$:RTE_USE_CBOOL> +) + +target_include_directories(rtekernels PRIVATE ${CMAKE_Fortran_MODULE_DIRECTORY}) diff --git a/rte-kernels/Make.depends b/rte-kernels/Make.depends deleted file mode 100644 index 6a8911e09..000000000 --- a/rte-kernels/Make.depends +++ /dev/null @@ -1,27 +0,0 @@ -RTE_FORTRAN_KERNELS = \ - mo_rte_kind.o \ - mo_rte_util_array.o \ - mo_rte_solver_kernels.o \ - mo_optical_props_kernels.o \ - mo_fluxes_broadband_kernels.o \ - -# -# Array utilities -# -mo_rte_util_array.o: mo_rte_kind.o mo_rte_util_array.F90 - -# -# Optical properties -# -mo_optical_props_kernels.o: mo_rte_kind.o mo_optical_props_kernels.F90 - -# -# Flux reduction -# -mo_fluxes_broadband_kernels.o : mo_rte_kind.o mo_fluxes_broadband_kernels.F90 - -# -# Radiative transfer -# -mo_rte_solver_kernels.o: mo_rte_kind.o mo_rte_util_array.o mo_rte_solver_kernels.F90 - diff --git a/rte-frontend/mo_rte_kind.F90 b/rte-kernels/mo_rte_kind.F90 similarity index 100% rename from rte-frontend/mo_rte_kind.F90 rename to rte-kernels/mo_rte_kind.F90 diff --git a/setup.sh b/setup.sh new file mode 100755 index 000000000..66bb3a651 --- /dev/null +++ b/setup.sh @@ -0,0 +1,23 @@ +rm -rf build +# conda --version +# conda env create -f environment-dev.yml + +FC=gfortran +FFLAGS='-ffree-line-length-none -m64 -std=f2008 -march=native -fbounds-check -fmodule-private -fimplicit-none -finit-real=nan' +RTE_ENABLE_SP=OFF +KERNEL_MODE=default +FAILURE_THRESHOLD='7.e-4' + +cmake -S . -B build -G "Ninja" \ + -DCMAKE_Fortran_COMPILER=$FC \ + -DCMAKE_Fortran_FLAGS="$FFLAGS" \ + -DRTE_ENABLE_SP=$RTE_ENABLE_SP \ + -DKERNEL_MODE=$KERNEL_MODE \ + -DBUILD_TESTING=ON \ + -DFAILURE_THRESHOLD=$FAILURE_THRESHOLD \ + -DCMAKE_BUILD_TYPE=Release + +# cmake --build build --target all --parallel + +# The --test-dir option is available only starting CMake 3.20: +# ctest -V --test-dir build diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt new file mode 100644 index 000000000..643374dcd --- /dev/null +++ b/tests/CMakeLists.txt @@ -0,0 +1,110 @@ +set(CMAKE_Fortran_MODULE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/modules) + +set(extensions_source_dir ${PROJECT_SOURCE_DIR}/extensions) + +add_library( + test_utils STATIC # cmake-format: sort + ${extensions_source_dir}/mo_compute_bc.F90 + ${extensions_source_dir}/mo_heating_rates.F90 + ${extensions_source_dir}/mo_rrtmgp_clr_all_sky.F90 + ${extensions_source_dir}/mo_zenith_angle_spherical_correction.F90 + ${extensions_source_dir}/solar_variability/mo_solar_variability.F90 + mo_gas_optics_defs_rrtmgp.F90 + mo_rcemip_profiles.F90 + mo_testing_utils.F90 +) + +target_include_directories( + test_utils + PUBLIC + $:${CMAKE_Fortran_MODULE_DIRECTORY}>> + ${NetCDF_Fortran_INCLUDE_DIR} +) + +target_link_libraries(test_utils PUBLIC rfmip_clear_utils) + +foreach( + test_executable IN + ITEMS # cmake-format: sort + check_equivalence + check_variants + rte_lw_solver_unit_tests + rte_optic_prop_unit_tests + rte_sw_solver_unit_tests + test_zenith_angle_spherical_correction +) + add_executable(${test_executable} ${test_executable}.F90) + target_link_libraries(${test_executable} PRIVATE test_utils) +endforeach() + +foreach( + test_executable IN + ITEMS # cmake-format: sort + rte_lw_solver_unit_tests + rte_optic_prop_unit_tests + rte_sw_solver_unit_tests +) + add_test(NAME ${test_executable} COMMAND ${test_executable}) +endforeach() + +foreach(g_value IN ITEMS 256 128) + add_test( + NAME check_equivalence_lw_g${g_value} + COMMAND + check_equivalence + ${RRTMGP_DATA}/examples/rfmip-clear-sky/inputs/multiple_input4MIPs_radiation_RFMIP_UColorado-RFMIP-1-2_none.nc + ${RRTMGP_DATA}/rrtmgp-gas-lw-g${g_value}.nc + ) + set_tests_properties( + check_equivalence_lw_g${g_value} + PROPERTIES FIXTURES_REQUIRED fetch_rrtmgp_data + ) +endforeach() + +foreach(g_value IN ITEMS 224 112) + add_test( + NAME check_equivalence_sw_g${g_value} + COMMAND + check_equivalence + ${RRTMGP_DATA}/examples/rfmip-clear-sky/inputs/multiple_input4MIPs_radiation_RFMIP_UColorado-RFMIP-1-2_none.nc + ${RRTMGP_DATA}/rrtmgp-gas-sw-g${g_value}.nc + ) + set_tests_properties( + check_equivalence_sw_g${g_value} + PROPERTIES FIXTURES_REQUIRED fetch_rrtmgp_data + ) +endforeach() + +add_test( + NAME check_variants_lw + COMMAND + check_variants + ${RRTMGP_DATA}/examples/rfmip-clear-sky/inputs/multiple_input4MIPs_radiation_RFMIP_UColorado-RFMIP-1-2_none.nc + lw_flux_variants.nc ${RRTMGP_DATA}/rrtmgp-gas-lw-g256.nc + ${RRTMGP_DATA}/rrtmgp-gas-lw-g128.nc +) +set_tests_properties( + check_variants_lw PROPERTIES FIXTURES_REQUIRED fetch_rrtmgp_data +) + +add_test( + NAME check_variants_sw + COMMAND + check_variants + ${RRTMGP_DATA}/examples/rfmip-clear-sky/inputs/multiple_input4MIPs_radiation_RFMIP_UColorado-RFMIP-1-2_none.nc + sw_flux_variants.nc ${RRTMGP_DATA}/rrtmgp-gas-sw-g224.nc + ${RRTMGP_DATA}/rrtmgp-gas-sw-g112.nc +) +set_tests_properties( + check_variants_sw PROPERTIES FIXTURES_REQUIRED fetch_rrtmgp_data +) + +add_custom_target( + validation-plots + COMMAND + ${Python3_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/validation-plots.py + --state_file + ${RRTMGP_DATA}/examples/rfmip-clear-sky/inputs/multiple_input4MIPs_radiation_RFMIP_UColorado-RFMIP-1-2_none.nc + --lw_vars_file lw_flux_variants.nc --sw_vars_file sw_flux_variants.nc + COMMENT "Generating validation plots" +) diff --git a/tests/Makefile b/tests/Makefile deleted file mode 100644 index a5060c623..000000000 --- a/tests/Makefile +++ /dev/null @@ -1,91 +0,0 @@ -#!/usr/bin/env make -# -# Location of RTE+RRTMGP libraries, module files. -# -RRTMGP_BUILD = $(RRTMGP_ROOT)/build -# -# RRTMGP library, module files -# -LDFLAGS += -L$(RRTMGP_BUILD) -LIBS += -lrrtmgp -lrte -FCINCLUDE += -I$(RRTMGP_BUILD) - -# netcdf Fortran module files has to be in the search path or added via environment variable FCINCLUDE e.g. -#FCINCLUDE += -I$(NFHOME)/include - -# netcdf C and Fortran libraries have to be in the search path or added via environment variable LDFLAGS e.g. -#LDFLAGS += -L$(NFHOME)/lib -L$(NCHOME)/lib -LIBS += -lnetcdff -lnetcdf - -VPATH = $(RRTMGP_ROOT)/examples:$(RRTMGP_ROOT)/examples/rfmip-clear-sky:$(RRTMGP_ROOT)/examples/all-sky -VPATH += $(RRTMGP_ROOT)/rrtmgp-frontend:$(RRTMGP_ROOT)/extensions:$(RRTMGP_ROOT)/:$(RRTMGP_ROOT)/extensions/solar_variability - -# Compilation rules -%.o: %.F90 - $(FC) $(FCFLAGS) $(FCINCLUDE) -c $< -%: %.o - $(FC) $(FCFLAGS) -o $@ $^ $(LDFLAGS) $(LIBS) - - -# -# Extra sources -- extensions to RRTMGP classes, shared infrastructure, local sources -# -ADDITIONS = mo_heating_rates.o mo_compute_bc.o mo_rrtmgp_clr_all_sky.o -ADDITIONS += mo_gas_optics_defs_rrtmgp.o -# File I/O -ADDITIONS += mo_simple_netcdf.o mo_rfmip_io.o -ADDITIONS += mo_testing_io.o mo_testing_utils.o -# Cloud optics -CLOUDS += mo_cloud_sampling.o mo_cloud_optics_rrtmgp.o mo_load_cloud_coefficients.o mo_garand_atmos_io.o -# Solar variability -ADDITIONS += mo_solar_variability.o - -# Many codes will need to be updated if the library changes -# LIB_DEPS = $(RRTMGP_BUILD)/librte.a $(RRTMGP_BUILD)/librrtmgp.a -# -# Targets -# -all: check_variants check_equivalence test_zenith_angle_spherical_correction rte_sw_solver_unit_tests rte_optic_prop_unit_tests rte_lw_solver_unit_tests - -check_equivalence: $(ADDITIONS) $(LIB_DEPS) check_equivalence.o -check_equivalence.o: $(ADDITIONS) $(LIB_DEPS) check_equivalence.F90 - -check_variants: $(ADDITIONS) $(LIB_DEPS) check_variants.o -check_variants.o: $(ADDITIONS) $(LIB_DEPS) check_variants.F90 - -test_zenith_angle_spherical_correction: mo_zenith_angle_spherical_correction.o mo_rcemip_profiles.o $(ADDITIONS) $(LIB_DEPS) test_zenith_angle_spherical_correction.o -test_zenith_angle_spherical_correction.o: mo_zenith_angle_spherical_correction.o mo_rcemip_profiles.o $(ADDITIONS) $(LIB_DEPS) test_zenith_angle_spherical_correction.F90 - -mo_testing_io.o: $(LIB_DEPS) mo_simple_netcdf.o mo_testing_io.F90 - -mo_cloud_optics_rrtmgp.o: $(LIB_DEPS) mo_cloud_optics_rrtmgp.F90 -mo_load_cloud_coefficients.o: $(LIB_DEPS) mo_simple_netcdf.o mo_cloud_optics_rrtmgp.o mo_load_cloud_coefficients.F90 -mo_cloud_sampling.o: $(LIB_DEPS) mo_cloud_sampling.F90 - -mo_gas_optics_defs_rrtmgp.o: $(LIB_DEPS) mo_testing_utils.o mo_simple_netcdf.o mo_gas_optics_defs_rrtmgp.F90 - -mo_load_coefficients.o: $(LIB_DEPS) mo_simple_netcdf.o mo_load_coefficients.F90 -mo_rfmip_io.o: $(LIB_DEPS) mo_simple_netcdf.o mo_rfmip_io.F90 -mo_simple_netcdf.o: $(LIB_DEPS) mo_simple_netcdf.F90 - -rte_optic_prop_unit_tests.o: $(LIB_DEPS) mo_testing_utils.o rte_optic_prop_unit_tests.F90 -rte_optic_prop_unit_tests : $(LIB_DEPS) mo_testing_utils.o rte_optic_prop_unit_tests.o - -rte_lw_solver_unit_tests.o: $(LIB_DEPS) mo_testing_utils.o rte_lw_solver_unit_tests.F90 -rte_lw_solver_unit_tests : $(LIB_DEPS) mo_testing_utils.o rte_lw_solver_unit_tests.o - -rte_sw_solver_unit_tests.o: $(LIB_DEPS) mo_testing_utils.o rte_sw_solver_unit_tests.F90 -rte_sw_solver_unit_tests : $(LIB_DEPS) mo_testing_utils.o rte_sw_solver_unit_tests.o - - -.PHONY: tests -tests: check_variants check_equivalence test_zenith_angle_spherical_correction rte_sw_solver_unit_tests rte_optic_prop_unit_tests rte_lw_solver_unit_tests - cp ${RRTMGP_DATA}/examples/rfmip-clear-sky/inputs/multiple_input4MIPs_radiation_RFMIP_UColorado-RFMIP-1-2_none.nc ./test_atmospheres.nc - $(RUN_CMD) bash all_tests.sh - - -check: - echo "Nothing to check in tests/" - -clean: - -rm clear_sky_regression *.o *.optrpt *.mod diff --git a/tests/all_tests.sh b/tests/all_tests.sh deleted file mode 100644 index 7752f4064..000000000 --- a/tests/all_tests.sh +++ /dev/null @@ -1,10 +0,0 @@ -set -eux -./rte_optic_prop_unit_tests -./rte_lw_solver_unit_tests -./rte_sw_solver_unit_tests -./check_equivalence test_atmospheres.nc ${RRTMGP_DATA}/rrtmgp-gas-lw-g256.nc -./check_equivalence test_atmospheres.nc ${RRTMGP_DATA}/rrtmgp-gas-lw-g128.nc -./check_equivalence test_atmospheres.nc ${RRTMGP_DATA}/rrtmgp-gas-sw-g224.nc -./check_equivalence test_atmospheres.nc ${RRTMGP_DATA}/rrtmgp-gas-sw-g112.nc -./check_variants test_atmospheres.nc ${RRTMGP_DATA}/rrtmgp-gas-lw-g256.nc ${RRTMGP_DATA}/rrtmgp-gas-lw-g128.nc -./check_variants test_atmospheres.nc ${RRTMGP_DATA}/rrtmgp-gas-sw-g224.nc ${RRTMGP_DATA}/rrtmgp-gas-sw-g112.nc \ No newline at end of file diff --git a/tests/check_equivalence.F90 b/tests/check_equivalence.F90 index f08ca8a26..86d064232 100644 --- a/tests/check_equivalence.F90 +++ b/tests/check_equivalence.F90 @@ -99,7 +99,7 @@ program rte_check_equivalence ! ! Inputs to RRTMGP ! - logical :: top_at_1, is_sw, is_lw + logical :: is_sw, is_lw integer :: ncol, nlay, nbnd, ngpt, nexp integer :: icol, ilay, ibnd, iloop, igas @@ -168,7 +168,6 @@ program rte_check_equivalence ! nbnd = gas_optics%get_nband() ngpt = gas_optics%get_ngpt() - top_at_1 = p_lay(1, 1) < p_lay(1, nlay) ! ---------------------------------------------------------------------------- ! ! Boundary conditions @@ -236,9 +235,9 @@ program rte_check_equivalence atmos, & lw_sources, & tlev = t_lev)) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes)) print *, " Default calculation" ! @@ -253,17 +252,17 @@ program rte_check_equivalence nullify(fluxes%flux_up) nullify(fluxes%flux_dn) allocate(fluxes%flux_net(ncol,nlay+1)) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes)) if(.not. allclose(fluxes%flux_net, ref_flux_dn-ref_flux_up) ) & call stop_on_err("Net fluxes don't match when computed alone") fluxes%flux_up => tst_flux_up(:,:) fluxes%flux_dn => tst_flux_dn(:,:) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes)) if(.not. allclose(fluxes%flux_net, ref_flux_dn-ref_flux_up) ) & call report_err("Net fluxes don't match when computed in tandem") @@ -293,36 +292,36 @@ program rte_check_equivalence ! atmos%tau(:,:,:) = 0.5_wp * atmos%tau(:,:,:) call stop_on_err(atmos%increment(atmos)) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes)) if(.not. allclose(tst_flux_up, ref_flux_up) .or. & .not. allclose(tst_flux_dn, ref_flux_dn) ) & call report_err(" halving/doubling fails") call increment_with_1scl(atmos) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes)) if(.not. allclose(tst_flux_up, ref_flux_up) .or. & .not. allclose(tst_flux_dn, ref_flux_dn) ) & call report_err(" Incrementing with 1scl fails") call increment_with_2str(atmos) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes)) if(.not. allclose(tst_flux_up, ref_flux_up) .or. & .not. allclose(tst_flux_dn, ref_flux_dn) ) & call report_err(" Incrementing with 2str fails") call increment_with_nstr(atmos) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes)) if(.not. allclose(tst_flux_up, ref_flux_up) .or. & .not. allclose(tst_flux_dn, ref_flux_dn) ) & @@ -332,10 +331,10 @@ program rte_check_equivalence ! ! Computing Jacobian shouldn't change net fluxes ! - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & - fluxes, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & + fluxes, & flux_up_Jac = jFluxUp)) if(.not. allclose(tst_flux_up, ref_flux_up) .or. & .not. allclose(tst_flux_dn, ref_flux_dn) ) & @@ -349,9 +348,9 @@ program rte_check_equivalence atmos, & lw_sources, & tlev = t_lev)) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes)) ! ! Comparision of fluxes with increased surface T aren't expected to match @@ -385,8 +384,7 @@ program rte_check_equivalence gas_concs, & atmos, & toa_flux)) - call stop_on_err(rte_sw(atmos, top_at_1, & - mu0, toa_flux, & + call stop_on_err(rte_sw(atmos, mu0, toa_flux, & sfc_alb_dir, sfc_alb_dif, & fluxes)) print *, " Default calculation" @@ -427,13 +425,12 @@ program rte_check_equivalence toa_flux)) atmos%tau(:,:,:) = 0.5_wp * atmos%tau(:,:,:) call stop_on_err(atmos%increment(atmos)) - call stop_on_err(rte_sw(atmos, top_at_1, & - mu0, toa_flux, & + call stop_on_err(rte_sw(atmos, mu0, toa_flux, & sfc_alb_dir, sfc_alb_dif, & fluxes)) if(.not. allclose(tst_flux_up, ref_flux_up, tol = 8._wp) .or. & - .not. allclose(tst_flux_dn, ref_flux_dn, tol = 10._wp) .or. & - .not. allclose(tst_flux_dir,ref_flux_dir,tol = 10._wp)) & + .not. allclose(tst_flux_dn, ref_flux_dn, tol = 12._wp) .or. & + .not. allclose(tst_flux_dir,ref_flux_dir,tol = 12._wp)) & call report_err(" halving/doubling fails") ! @@ -445,13 +442,13 @@ program rte_check_equivalence atmos, & toa_flux)) call increment_with_1scl(atmos) - call stop_on_err(rte_sw(atmos, top_at_1, & + call stop_on_err(rte_sw(atmos, & mu0, toa_flux, & sfc_alb_dir, sfc_alb_dif, & fluxes)) if(.not. allclose(tst_flux_up, ref_flux_up, tol = 8._wp) .or. & - .not. allclose(tst_flux_dn, ref_flux_dn, tol = 10._wp) .or. & - .not. allclose(tst_flux_dir,ref_flux_dir,tol = 10._wp)) & + .not. allclose(tst_flux_dn, ref_flux_dn, tol = 12._wp) .or. & + .not. allclose(tst_flux_dir,ref_flux_dir,tol = 12._wp)) & call report_err(" Incrementing with 1scl fails") call stop_on_err(gas_optics%gas_optics(p_lay, p_lev, & @@ -460,13 +457,13 @@ program rte_check_equivalence atmos, & toa_flux)) call increment_with_2str(atmos) - call stop_on_err(rte_sw(atmos, top_at_1, & + call stop_on_err(rte_sw(atmos, & mu0, toa_flux, & sfc_alb_dir, sfc_alb_dif, & fluxes)) if(.not. allclose(tst_flux_up, ref_flux_up, tol = 8._wp) .or. & - .not. allclose(tst_flux_dn, ref_flux_dn, tol = 10._wp) .or. & - .not. allclose(tst_flux_dir,ref_flux_dir,tol = 10._wp)) & + .not. allclose(tst_flux_dn, ref_flux_dn, tol = 12._wp) .or. & + .not. allclose(tst_flux_dir,ref_flux_dir,tol = 12._wp)) & call report_err(" Incrementing with 2str fails") call stop_on_err(gas_optics%gas_optics(p_lay, p_lev, & @@ -475,13 +472,13 @@ program rte_check_equivalence atmos, & toa_flux)) call increment_with_nstr(atmos) - call stop_on_err(rte_sw(atmos, top_at_1, & + call stop_on_err(rte_sw(atmos, & mu0, toa_flux, & sfc_alb_dir, sfc_alb_dif, & fluxes)) if(.not. allclose(tst_flux_up, ref_flux_up, tol = 8._wp) .or. & - .not. allclose(tst_flux_dn, ref_flux_dn, tol = 10._wp) .or. & - .not. allclose(tst_flux_dir,ref_flux_dir,tol = 10._wp)) & + .not. allclose(tst_flux_dn, ref_flux_dn, tol = 12._wp) .or. & + .not. allclose(tst_flux_dir,ref_flux_dir,tol = 12._wp)) & call report_err(" Incrementing with nstr fails") print *, " Incrementing" end if @@ -507,7 +504,6 @@ subroutine lw_clear_sky_vr t_lay (:,:) = t_lay (:, nlay :1:-1) p_lev (:,:) = p_lev (:,(nlay+1):1:-1) t_lev (:,:) = t_lev (:,(nlay+1):1:-1) - top_at_1 = .not. top_at_1 ! ! No direct access to gas concentrations so use the classes ! This also tests otherwise uncovered routines for ty_gas_concs @@ -526,9 +522,9 @@ subroutine lw_clear_sky_vr atmos, & lw_sources, & tlev=t_lev)) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes)) tst_flux_up(:,:) = tst_flux_up(:,(nlay+1):1:-1) tst_flux_dn(:,:) = tst_flux_dn(:,(nlay+1):1:-1) @@ -536,7 +532,6 @@ subroutine lw_clear_sky_vr t_lay (:,:) = t_lay (:, nlay :1:-1) p_lev (:,:) = p_lev (:,(nlay+1):1:-1) t_lev (:,:) = t_lev (:,(nlay+1):1:-1) - top_at_1 = .not. top_at_1 end subroutine lw_clear_sky_vr ! ---------------------------------------------------------------------------- ! @@ -566,8 +561,8 @@ subroutine lw_clear_sky_subset colE = i * ncol/2 call stop_on_err(atmos%get_subset (colS, ncol/2, atmos_subset)) call stop_on_err(lw_sources%get_subset(colS, ncol/2, sources_subset)) - call stop_on_err(rte_lw(atmos_subset, top_at_1, & - sources_subset, & + call stop_on_err(rte_lw(atmos_subset, & + sources_subset, & sfc_emis(:, colS:colE), & fluxes)) tst_flux_up(colS:colE,:) = up @@ -594,7 +589,6 @@ subroutine sw_clear_sky_vr p_lay (:,:) = p_lay (:, nlay :1:-1) t_lay (:,:) = t_lay (:, nlay :1:-1) p_lev (:,:) = p_lev (:,(nlay+1):1:-1) - top_at_1 = .not. top_at_1 ! ! No direct access to gas concentrations so use the classes ! This also tests otherwise uncovered routines for ty_gas_concs @@ -612,8 +606,7 @@ subroutine sw_clear_sky_vr gas_concs_vr, & atmos, & toa_flux)) - call stop_on_err(rte_sw(atmos, top_at_1, & - mu0, toa_flux, & + call stop_on_err(rte_sw(atmos, mu0, toa_flux, & sfc_alb_dir, sfc_alb_dif, & fluxes)) ! @@ -626,7 +619,6 @@ subroutine sw_clear_sky_vr p_lay (:,:) = p_lay (:, nlay :1:-1) t_lay (:,:) = t_lay (:, nlay :1:-1) p_lev (:,:) = p_lev (:,(nlay+1):1:-1) - top_at_1 = .not. top_at_1 end subroutine sw_clear_sky_vr ! ---------------------------------------------------------------------------- ! @@ -645,8 +637,7 @@ subroutine sw_clear_sky_tsi gas_concs, & atmos, & toa_flux)) - call stop_on_err(rte_sw(atmos, top_at_1, & - mu0, toa_flux, & + call stop_on_err(rte_sw(atmos, mu0, toa_flux, & sfc_alb_dir, sfc_alb_dif, & fluxes)) tst_flux_up (:,:) = tst_flux_up (:,:) / tsi_scale diff --git a/tests/check_variants.F90 b/tests/check_variants.F90 index a57b9ec84..385ca04a8 100644 --- a/tests/check_variants.F90 +++ b/tests/check_variants.F90 @@ -47,7 +47,8 @@ program rte_clear_sky_regression read_and_block_lw_bc, read_and_block_sw_bc, determine_gas_names use mo_simple_netcdf, only: get_dim_size, read_field use mo_heating_rates, only: compute_heating_rate - use mo_testing_io, only: write_broadband_field + use netcdf + use mo_simple_netcdf implicit none ! ---------------------------------------------------------------------------------- ! Variables @@ -100,7 +101,7 @@ program rte_clear_sky_regression ! ! Inputs to RRTMGP ! - logical :: top_at_1, is_sw, is_lw + logical :: is_sw, is_lw integer :: ncol, nlay, nbnd, ngpt, nexp integer :: icol, ilay, ibnd, iloop, igas @@ -110,7 +111,8 @@ program rte_clear_sky_regression character(len=32 ), & dimension(:), allocatable :: kdist_gas_names, rfmip_gas_games - character(len=256) :: input_file = "", gas_optics_file = "", gas_optics_file_2 = "" + character(len=256) :: input_file = "", output_file = "", gas_optics_file = "", gas_optics_file_2 = "" + integer :: ncid, dimid ! ---------------------------------------------------------------------------------- ! Code ! ---------------------------------------------------------------------------------- @@ -118,11 +120,12 @@ program rte_clear_sky_regression ! Parse command line for any file names, block size ! nUserArgs = command_argument_count() - if (nUserArgs < 2) call stop_on_err("Need to supply input_file gas_optics_file [gas_optics_file_2]") + if (nUserArgs < 3) call stop_on_err("Need to supply input_file output_file gas_optics_file [gas_optics_file_2]") if (nUserArgs >= 1) call get_command_argument(1,input_file) - if (nUserArgs >= 2) call get_command_argument(2,gas_optics_file) - if (nUserArgs >= 3) call get_command_argument(3,gas_optics_file_2) - if (nUserArgs > 4) print *, "Ignoring command line arguments beyond the first four..." + if (nUserArgs >= 2) call get_command_argument(2,output_file) + if (nUserArgs >= 3) call get_command_argument(3,gas_optics_file) + if (nUserArgs >= 4) call get_command_argument(4,gas_optics_file_2) + if (nUserArgs > 5) print *, "Ignoring command line arguments beyond the first four..." if(trim(input_file) == '-h' .or. trim(input_file) == "--help") then call stop_on_err("clear_sky_regression input_file absorption_coefficients_file") end if @@ -166,7 +169,6 @@ program rte_clear_sky_regression ! nbnd = gas_optics%get_nband() ngpt = gas_optics%get_ngpt() - top_at_1 = p_lay(1, 1) < p_lay(1, nlay) ! ---------------------------------------------------------------------------- ! ! Boundary conditions @@ -202,6 +204,21 @@ program rte_clear_sky_regression ! Fluxes ! allocate(flux_up(ncol,nlay+1), flux_dn(ncol,nlay+1), flux_dir(ncol,nlay+1)) + ! ---------------------------------------------------------------------------- + ! + ! Create output file and site, level, layer dimensions + ! + if(nf90_create(trim(output_file), NF90_CLOBBER, ncid) /= NF90_NOERR) & + call stop_on_err("write_fluxes: can't create file " // trim(output_file)) + if(nf90_def_dim(ncid, "site", ncol, dimid) /= NF90_NOERR) & + call stop_on_err("fail to create 'site' dimension") + if(nf90_def_dim(ncid, "level", nlay+1, dimid) /= NF90_NOERR) & + call stop_on_err("fail to create 'level' dimension") + if(nf90_def_dim(ncid, "layer", nlay, dimid) /= NF90_NOERR) & + call stop_on_err("fail to create 'layer' dimension") + if(nf90_enddef(ncid) /= NF90_NOERR) & + call stop_on_err("fail to to end redefinition??") + ! ---------------------------------------------------------------------------- ! ! Solvers @@ -245,6 +262,7 @@ program rte_clear_sky_regression call sw_clear_sky_alt end if end if + ncid = nf90_close(ncid) contains ! ---------------------------------------------------------------------------- ! @@ -261,27 +279,27 @@ subroutine lw_clear_sky_default atmos, & lw_sources, & tlev = t_lev)) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes)) - call write_broadband_field(input_file, flux_up, "lw_flux_up", "LW flux up") - call write_broadband_field(input_file, flux_dn, "lw_flux_dn", "LW flux dn") - call write_broadband_field(input_file, flux_net, "lw_flux_net", "LW flux net") + call write_broadband_field(flux_up, "lw_flux_up", "LW flux up") + call write_broadband_field(flux_dn, "lw_flux_dn", "LW flux dn") + call write_broadband_field(flux_net, "lw_flux_net", "LW flux net") call stop_on_err(compute_heating_rate(flux_up, flux_dn, p_lev, heating_rate)) - call write_broadband_field(input_file, heating_rate, & + call write_broadband_field(heating_rate, & "lw_flux_hr_default", "LW heating rate", vert_dim_name = "layer") ! ! Test for mo_fluxes_broadband for computing only net flux ! nullify(fluxes%flux_up) nullify(fluxes%flux_dn) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes)) - call write_broadband_field(input_file, flux_net, "lw_flux_net_2", "LW flux net, direct") + call write_broadband_field(flux_net, "lw_flux_net_2", "LW flux net, direct") fluxes%flux_up => flux_up fluxes%flux_dn => flux_dn nullify(fluxes%flux_net) @@ -297,12 +315,12 @@ subroutine lw_clear_sky_notlev gas_concs, & atmos, & lw_sources)) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes)) - call write_broadband_field(input_file, flux_up, "lw_flux_up_notlev", "LW flux up, no level temperatures") - call write_broadband_field(input_file, flux_dn, "lw_flux_dn_notlev", "LW flux dn, no level temperatures") + call write_broadband_field(flux_up, "lw_flux_up_notlev", "LW flux up, no level temperatures") + call write_broadband_field(flux_dn, "lw_flux_dn_notlev", "LW flux dn, no level temperatures") end subroutine lw_clear_sky_notlev ! ---------------------------------------------------------------------------- ! @@ -315,12 +333,12 @@ subroutine lw_clear_sky_3ang atmos, & lw_sources, & tlev = t_lev)) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes, n_gauss_angles=3)) - call write_broadband_field(input_file, flux_up, "lw_flux_up_3ang", "LW flux up, three quadrature angles") - call write_broadband_field(input_file, flux_dn, "lw_flux_dn_3ang", "LW flux dn, three quadrature angles") + call write_broadband_field(flux_up, "lw_flux_up_3ang", "LW flux up, three quadrature angles") + call write_broadband_field(flux_dn, "lw_flux_dn_3ang", "LW flux dn, three quadrature angles") end subroutine lw_clear_sky_3ang ! ---------------------------------------------------------------------------- ! @@ -335,12 +353,12 @@ subroutine lw_clear_sky_optangle lw_sources, & tlev = t_lev)) call stop_on_err(gas_optics%compute_optimal_angles(atmos, lw_Ds)) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes, lw_Ds=lw_Ds)) - call write_broadband_field(input_file, flux_up, "lw_flux_up_optang", "LW flux up, single optimal angles") - call write_broadband_field(input_file, flux_dn, "lw_flux_dn_optang", "LW flux dn, single optimal angles") + call write_broadband_field(flux_up, "lw_flux_up_optang", "LW flux up, single optimal angles") + call write_broadband_field(flux_dn, "lw_flux_dn_optang", "LW flux dn, single optimal angles") end subroutine lw_clear_sky_optangle ! ---------------------------------------------------------------------------- ! @@ -355,14 +373,14 @@ subroutine lw_clear_sky_jaco atmos, & lw_sources, & tlev=t_lev)) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & - fluxes, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & + fluxes, & flux_up_Jac = jFluxUp)) - call write_broadband_field(input_file, flux_up, "lw_flux_up_jaco", "LW flux up, computing Jaobians") - call write_broadband_field(input_file, flux_dn, "lw_flux_dn_jaco", "LW flux dn, computing Jaobians") - call write_broadband_field(input_file, jFluxUp, "lw_jaco_up" , "Jacobian of LW flux up to surface temperature") + call write_broadband_field(flux_up, "lw_flux_up_jaco", "LW flux up, computing Jaobians") + call write_broadband_field(flux_dn, "lw_flux_dn_jaco", "LW flux dn, computing Jaobians") + call write_broadband_field(jFluxUp, "lw_jaco_up" , "Jacobian of LW flux up to surface temperature") call stop_on_err(gas_optics%gas_optics(p_lay, p_lev, & t_lay, sfc_t + 1._wp, & @@ -370,12 +388,12 @@ subroutine lw_clear_sky_jaco atmos, & lw_sources, & tlev=t_lev)) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes)) - call write_broadband_field(input_file, flux_up, "lw_flux_up_stp1", "LW flux up, surface T+1K") - call write_broadband_field(input_file, flux_dn, "lw_flux_dn_stp1", "LW flux dn, surface T+1K") + call write_broadband_field(flux_up, "lw_flux_up_stp1", "LW flux up, surface T+1K") + call write_broadband_field(flux_dn, "lw_flux_dn_stp1", "LW flux dn, surface T+1K") end subroutine lw_clear_sky_jaco ! ---------------------------------------------------------------------------- ! @@ -391,19 +409,19 @@ subroutine lw_clear_sky_2str atmos, & lw_sources, & tlev=t_lev)) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes)) - call write_broadband_field(input_file, flux_up, "lw_flux_up_1rescl", "LW flux up, clear-sky _1rescl") - call write_broadband_field(input_file, flux_dn, "lw_flux_dn_1rescl", "LW flux dn, clear-sky _1rescl") + call write_broadband_field(flux_up, "lw_flux_up_1rescl", "LW flux up, clear-sky _1rescl") + call write_broadband_field(flux_dn, "lw_flux_dn_1rescl", "LW flux dn, clear-sky _1rescl") - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes, use_2stream=.true.)) - call write_broadband_field(input_file, flux_up, "lw_flux_up_2str", "LW flux up, clear-sky _2str") - call write_broadband_field(input_file, flux_dn, "lw_flux_dn_2str", "LW flux dn, clear-sky _2str") + call write_broadband_field(flux_up, "lw_flux_up_2str", "LW flux up, clear-sky _2str") + call write_broadband_field(flux_dn, "lw_flux_dn_2str", "LW flux dn, clear-sky _2str") end subroutine lw_clear_sky_2str ! ---------------------------------------------------------------------------- @@ -419,28 +437,28 @@ subroutine lw_clear_sky_alt atmos, & lw_sources, & tlev = t_lev)) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes)) - call write_broadband_field(input_file, flux_up, "lw_flux_up_alt", "LW flux up, fewer g-points") - call write_broadband_field(input_file, flux_dn, "lw_flux_dn_alt", "LW flux dn, fewer g-points") - call write_broadband_field(input_file, flux_net, "lw_flux_net_alt", "LW flux ne, fewer g-pointst") + call write_broadband_field(flux_up, "lw_flux_up_alt", "LW flux up, fewer g-points") + call write_broadband_field(flux_dn, "lw_flux_dn_alt", "LW flux dn, fewer g-points") + call write_broadband_field(flux_net, "lw_flux_net_alt", "LW flux ne, fewer g-pointst") call stop_on_err(compute_heating_rate(flux_up, flux_dn, p_lev, heating_rate)) - call write_broadband_field(input_file, heating_rate, & + call write_broadband_field(heating_rate, & "lw_flux_hr_alt", "LW heating rate, fewer g-points", & vert_dim_name = "layer") call stop_on_err(gas_optics%compute_optimal_angles(atmos, lw_Ds)) - call stop_on_err(rte_lw(atmos, top_at_1, & - lw_sources, & - sfc_emis, & + call stop_on_err(rte_lw(atmos, & + lw_sources, & + sfc_emis, & fluxes, lw_Ds=lw_Ds)) - call write_broadband_field(input_file, flux_up, "lw_flux_up_alt_oa", "LW flux up, fewer g-points, opt. angle") - call write_broadband_field(input_file, flux_dn, "lw_flux_dn_alt_oa", "LW flux dn, fewer g-points, opt. angle") - call write_broadband_field(input_file, flux_net, "lw_flux_net_alt_oa", "LW flux ne, fewer g-points, opt. angle") + call write_broadband_field(flux_up, "lw_flux_up_alt_oa", "LW flux up, fewer g-points, opt. angle") + call write_broadband_field(flux_dn, "lw_flux_dn_alt_oa", "LW flux dn, fewer g-points, opt. angle") + call write_broadband_field(flux_net, "lw_flux_net_alt_oa", "LW flux ne, fewer g-points, opt. angle") call stop_on_err(compute_heating_rate(flux_up, flux_dn, p_lev, heating_rate)) - call write_broadband_field(input_file, heating_rate, & + call write_broadband_field(heating_rate, & "lw_flux_hr_alt_oa", "LW heating rate, fewer g-points, opt. angle", & vert_dim_name = "layer") call gas_optics%finalize() @@ -470,8 +488,7 @@ subroutine sw_clear_sky_default rfmip_tsi_scale(:,:) = spread(tsi_3d(:,1)/rrtmgp_tsi, dim=2, ncopies=ngpt) toa_flux(:,:) = toa_flux(:,:) * rfmip_tsi_scale(:,:) - call stop_on_err(rte_sw(atmos, top_at_1, & - mu0, toa_flux, & + call stop_on_err(rte_sw(atmos, mu0, toa_flux, & sfc_alb_dir, sfc_alb_dif, & fluxes)) ! @@ -481,8 +498,8 @@ subroutine sw_clear_sky_default flux_up = 0._wp flux_dn = 0._wp end where - call write_broadband_field(input_file, flux_up, "sw_flux_up", "SW flux up") - call write_broadband_field(input_file, flux_dn, "sw_flux_dn", "SW flux dn") + call write_broadband_field(flux_up, "sw_flux_up", "SW flux up") + call write_broadband_field(flux_dn, "sw_flux_dn", "SW flux dn") end subroutine sw_clear_sky_default ! ---------------------------------------------------------------------------- subroutine sw_clear_sky_alt @@ -501,8 +518,7 @@ subroutine sw_clear_sky_alt rfmip_tsi_scale(:,:) = spread(tsi_3d(:,1)/rrtmgp_tsi, dim=2, ncopies=gas_optics%get_ngpt()) toa_flux(:,:) = toa_flux(:,:) * rfmip_tsi_scale(:,:) - call stop_on_err(rte_sw(atmos, top_at_1, & - mu0, toa_flux, & + call stop_on_err(rte_sw(atmos, mu0, toa_flux, & sfc_alb_dir, sfc_alb_dif, & fluxes)) ! @@ -512,8 +528,8 @@ subroutine sw_clear_sky_alt flux_up = 0._wp flux_dn = 0._wp end where - call write_broadband_field(input_file, flux_up, "sw_flux_up_alt", "SW flux up, fewer g-points") - call write_broadband_field(input_file, flux_dn, "sw_flux_dn_alt", "SW flux dn, fewer g-points") + call write_broadband_field(flux_up, "sw_flux_up_alt", "SW flux up, fewer g-points") + call write_broadband_field(flux_dn, "sw_flux_dn_alt", "SW flux dn, fewer g-points") end subroutine sw_clear_sky_alt ! ---------------------------------------------------------------------------- subroutine make_optical_props_1scl(gas_optics) @@ -555,4 +571,44 @@ subroutine make_optical_props_2str(gas_optics) end select end subroutine make_optical_props_2str ! ---------------------------------------------------------------------------- -end program rte_clear_sky_regression + subroutine write_broadband_field(field, field_name, field_description, col_dim_name, vert_dim_name) + ! + ! Write a field defined by column and some vertical dimension (lev or lay)) + ! + real(wp), dimension(:,:), intent(in) :: field + character(len=*), intent(in) :: field_name, field_description + character(len=*), optional, & + intent(in) ::col_dim_name, vert_dim_name + ! ------------------- + integer :: varid, ncol, nlev + ! + ! Names of column (first) and vertical (second) dimension. + ! Because they are used in an array constuctor the need to have the same number of characters + ! + character(len=32) :: cdim, vdim + ! ------------------- + cdim = "site " + vdim = "level" + if(present(col_dim_name)) cdim = col_dim_name + if(present(vert_dim_name)) vdim = vert_dim_name + + ncol = size(field, dim=1) + nlev = size(field, dim=2) + + call create_var(ncid, trim(field_name), [cdim, vdim], [ncol, nlev]) + call stop_on_err(write_field(ncid, trim(field_name), field)) + ! + ! Adding descriptive text as an attribute means knowing the varid + ! + if(nf90_inq_varid(ncid, trim(field_name), varid) /= NF90_NOERR) & + call stop_on_err("Can't find variable " // trim(field_name)) + if(nf90_redef(ncid) /= NF90_NOERR) & + call stop_on_err("write_broadband_field: can't put file into redefine mode") + if(nf90_put_att(ncid, varid, "description", trim(field_description)) /= NF90_NOERR) & + call stop_on_err("Can't write 'description' attribute to variable " // trim(field_name)) + if(nf90_enddef(ncid) /= NF90_NOERR) & + call stop_on_err("write_broadband_field: fail to to end redefinition??") + + end subroutine write_broadband_field + ! ---------------------------------------------------------------------------- +end program diff --git a/tests/intel-codecov.sh b/tests/intel-codecov.sh deleted file mode 100644 index 90422fda2..000000000 --- a/tests/intel-codecov.sh +++ /dev/null @@ -1,75 +0,0 @@ -#!/bin/bash -ulimit -s hard -export FC=ifort -export FCFLAGS="-m64 -prof-gen=srcpos -g -traceback -heap-arrays -assume realloc_lhs -extend-source 132" -# -# Intel specific - where will the profiling files be generated? -# -export RRTMGP_ROOT=`cd ..;pwd` -export RRTMGP_BUILD=${RRTMGP_ROOT}/build -export PROF_DIR=$PWD -# -# Environment variables for netCDF Fortran and C installations -# -export NFHOME=${HOME}/Applications/${FC} -export NCHOME=/opt/local - -# -# An Anaconda environent with modules needed for other python scripts -# (xarray, matplotlib, ...) -# -source activate pangeo -cd ${PROF_DIR} -rm -rf *.dyn pgopti.* CODE_COVERAGE.html CodeCoverage/ -# -# Build RTE+RRTMGP librarues -# -make -C ${RRTMGP_BUILD} clean -make -C ${RRTMGP_BUILD} -j 4 || exit 1 - -# -# Build and run RFMIP examples -# -cd ${RRTMGP_ROOT}/examples/rfmip-clear-sky || exit 1 -export FCFLAGS+=" -I${RRTMGP_BUILD} -I${NFHOME}/include" -export LDFLAGS+=" -L${RRTMGP_BUILD} -L${NFHOME}/lib -L${NCHOME}/lib -lrte -lrrtmgp -lnetcdff -lnetcdf" -make clean || exit 1 -make -j 4 || exit 1 -python ./stage_files.py -python ./run-rfmip-examples.py -python ./compare-to-reference.py --fail=7.e-4 -make clean -# -# Build and run all-sky examples -# -cd ${RRTMGP_ROOT}/examples/all-sky || exit 1 -make clean || exit 1 -make -j 4 || exit 1 -python ./run-allsky-example.py -python ./compare-to-reference.py -make clean -# -# Build and run regression tests -# -cd ${PROF_DIR} || exit 1 -make clean || exit 1 -make -j 4 || exit 1 -cp ${RRTMGP_DATA}/examples/rfmip-clear-sky/inputs/multiple_input4MIPs_radiation_RFMIP_UColorado-RFMIP-1-2_none.nc test_atmospheres.nc -./clear_sky_regression test_atmospheres.nc ${RRTMGP_DATA}/rrtmgp-gas-lw-g256.nc -./clear_sky_regression test_atmospheres.nc ${RRTMGP_DATA}/rrtmgp-gas-sw-g224.nc -# Need to repeat for openacc-kernels - -# -# Merge -# -cd ${PROF_DIR} -profmerge -a -echo " -mo_ -~mo_simple_netcdf -~mo_rfmip_io -~mo_testing_io -~mo_garand_atmos_io -~mo_load" > intel_codecov_filter.txt -codecov -prj rte-rrtmgp -comp intel_codecov_filter.txt -rm intel_codecov_filter.txt diff --git a/tests/mo_testing_io.F90 b/tests/mo_testing_io.F90 deleted file mode 100644 index 8444cac56..000000000 --- a/tests/mo_testing_io.F90 +++ /dev/null @@ -1,60 +0,0 @@ -!-------------------------------------------------------------------------------------------------------------------- -module mo_testing_io - ! - ! RTE+RRTMGP modules - ! - use mo_rte_kind, only: wp - ! - ! NetCDF I/O routines, shared with other RTE+RRTMGP examples - ! - use mo_simple_netcdf, only: write_field, create_dim, create_var - use netcdf - implicit none - private - public :: write_broadband_field -contains - !-------------------------------------------------------------------------------------------------------------------- - subroutine write_broadband_field(fileName, field, field_name, field_description, col_dim_name, vert_dim_name) - ! - ! Write a field defined by column and some vertical dimension (lev or lay)) - ! - character(len=*), intent(in) :: fileName - real(wp), dimension(:,:), intent(in) :: field - character(len=*), intent(in) :: field_name, field_description - character(len=*), optional, & - intent(in) ::col_dim_name, vert_dim_name - ! ------------------- - integer :: ncid, varid, ncol, nlev - ! - ! Names of column (first) and vertical (second) dimension. - ! Because they are used in an array constuctor the need to have the same number of characters - ! - character(len=32) :: cdim, vdim - ! ------------------- - cdim = "site " - vdim = "level" - if(present(col_dim_name)) cdim = col_dim_name - if(present(vert_dim_name)) vdim = vert_dim_name - - if(nf90_open(trim(fileName), NF90_WRITE, ncid) /= NF90_NOERR) & - call stop_on_err("write_fluxes: can't open file " // trim(fileName)) - - ncol = size(field, dim=1) - nlev = size(field, dim=2) - - call create_dim(ncid, trim(cdim), ncol) - call create_var(ncid, trim(field_name), [cdim, vdim], [ncol, nlev]) - call stop_on_err(write_field(ncid, trim(field_name), field)) - ! - ! Adding descriptive text as an attribute means knowing the varid - ! - if(nf90_inq_varid(ncid, trim(field_name), varid) /= NF90_NOERR) & - call stop_on_err("Can't find variable " // trim(field_name)) - if(nf90_put_att(ncid, varid, "description", trim(field_description)) /= NF90_NOERR) & - call stop_on_err("Can't write 'description' attribute to variable " // trim(field_name)) - - ncid = nf90_close(ncid) - - end subroutine write_broadband_field - !-------------------------------------------------------------------------------------------------------------------- -end module mo_testing_io diff --git a/tests/mo_testing_utils.F90 b/tests/mo_testing_utils.F90 index 07adffa39..95dce46c0 100644 --- a/tests/mo_testing_utils.F90 +++ b/tests/mo_testing_utils.F90 @@ -256,6 +256,7 @@ subroutine vr(atmos, sources) ! ----------------------- nlay = atmos%get_nlay() + call atmos%set_top_at_1(.not. atmos%top_is_at_1()) atmos%tau(:,:,:) = atmos%tau(:,nlay:1:-1,:) select type (atmos) diff --git a/tests/rte_lw_solver_unit_tests.F90 b/tests/rte_lw_solver_unit_tests.F90 index baf08fd3f..ed42440e4 100644 --- a/tests/rte_lw_solver_unit_tests.F90 +++ b/tests/rte_lw_solver_unit_tests.F90 @@ -95,7 +95,7 @@ program rte_lw_solver_unit_tests fluxes%flux_up => ref_flux_up (:,:) fluxes%flux_dn => ref_flux_dn (:,:) fluxes%flux_net => ref_flux_net(:,:) - call stop_on_err(rte_lw(lw_atmos, top_at_1, & + call stop_on_err(rte_lw(lw_atmos, & lw_sources, & sfc_emis, & fluxes)) @@ -116,7 +116,7 @@ program rte_lw_solver_unit_tests ! nullify(fluxes%flux_up) nullify(fluxes%flux_dn) - call stop_on_err(rte_lw(lw_atmos, top_at_1, & + call stop_on_err(rte_lw(lw_atmos, & lw_sources, sfc_emis,& fluxes)) call check_fluxes(ref_flux_net, ref_flux_dn-ref_flux_up, & @@ -126,7 +126,7 @@ program rte_lw_solver_unit_tests ! fluxes%flux_up => tst_flux_up (:,:) fluxes%flux_dn => tst_flux_dn (:,:) - call stop_on_err(rte_lw(lw_atmos, top_at_1, & + call stop_on_err(rte_lw(lw_atmos, & lw_sources, sfc_emis, & fluxes)) call check_fluxes(ref_flux_net, tst_flux_dn-tst_flux_up, & @@ -149,7 +149,7 @@ program rte_lw_solver_unit_tests print *, " Vertical orientation invariance" call gray_rad_equil(sfc_t, total_tau, nlay, top_at_1, lw_atmos, lw_sources) call vr(lw_atmos, lw_sources) - call stop_on_err(rte_lw(lw_atmos, .not. top_at_1, & + call stop_on_err(rte_lw(lw_atmos, & lw_sources, sfc_emis, & fluxes)) ! @@ -168,7 +168,7 @@ program rte_lw_solver_unit_tests ! print *, " Jacobian" call gray_rad_equil(sfc_t, total_tau, nlay, top_at_1, lw_atmos, lw_sources) - call stop_on_err(rte_lw(lw_atmos, top_at_1, & + call stop_on_err(rte_lw(lw_atmos, & lw_sources, & sfc_emis, & fluxes, & @@ -180,7 +180,7 @@ program rte_lw_solver_unit_tests ! lw_sources%sfc_source (:,1) = sigma/pi * (sfc_t + 1._wp)**4 lw_sources%sfc_source_Jac(:,1) = 4._wp * sigma/pi * (sfc_t + 1._wp)**3 - call stop_on_err(rte_lw(lw_atmos, top_at_1, & + call stop_on_err(rte_lw(lw_atmos, & lw_sources, & sfc_emis, & fluxes)) @@ -200,8 +200,10 @@ program rte_lw_solver_unit_tests sw_atmos%tau = lw_atmos%tau sw_atmos%ssa = 0._wp sw_atmos%g = 0._wp + call sw_atmos%set_top_at_1(lw_atmos%top_is_at_1()) - call stop_on_err(rte_lw(sw_atmos, top_at_1, & + + call stop_on_err(rte_lw(sw_atmos, & lw_sources, & sfc_emis, & fluxes, & @@ -214,7 +216,7 @@ program rte_lw_solver_unit_tests ! Specifying diffusivity angle ! print *, " Specified transport angle" - call stop_on_err(rte_lw(lw_atmos, top_at_1, & + call stop_on_err(rte_lw(lw_atmos, & lw_sources, & sfc_emis, & fluxes, & @@ -258,7 +260,7 @@ subroutine gray_rad_equil(sfc_t, total_tau, nlay, top_at_1, atmos, sources) ! Divide optical depth evenly among layers ! atmos%tau(1:ncol,1:nlay,1) = spread(total_tau(1:ncol)/real(nlay, wp), dim=2, ncopies=nlay) - + call atmos%set_top_at_1(top_at_1) ! ! Longwave sources - for broadband these are sigma/pi T^4 ! (isotropic radiation) @@ -271,17 +273,17 @@ subroutine gray_rad_equil(sfc_t, total_tau, nlay, top_at_1, atmos, sources) ! ! Calculation with top_at_1 ! - ilay = 1 - sources%lev_source(:,ilay, 1) = 0.5_wp/pi * olr(:) - do ilay = 2, nlay+1 - sources%lev_source(:,ilay, 1) = 0.5_wp/pi * olr(:) * & - (1._wp + D * sum(atmos%tau(:,:ilay-1,1),dim=2)) - ! - ! The source is linear in optical depth so layer source is average of edges - ! - sources%lay_source(:,ilay-1,1) = 0.5_wp * (sources%lev_source(:,ilay, 1) + & - sources%lev_source(:,ilay-1,1)) - end do + ilay = 1 + sources%lev_source(:,ilay, 1) = 0.5_wp/pi * olr(:) + do ilay = 2, nlay+1 + sources%lev_source(:,ilay, 1) = 0.5_wp/pi * olr(:) * & + (1._wp + D * sum(atmos%tau(:,:ilay-1,1),dim=2)) + ! + ! The source is linear in optical depth so layer source is average of edges + ! + sources%lay_source(:,ilay-1,1) = 0.5_wp * (sources%lev_source(:,ilay, 1) + & + sources%lev_source(:,ilay-1,1)) + end do if (.not. top_at_1) then ! ! Reverse vertical ordering of source functions @@ -372,8 +374,8 @@ subroutine clear_sky_subset(atmos, sources, sfc_emis, flux_up, flux_dn) colE = i * ncol/2 call stop_on_err(atmos%get_subset (colS, ncol/2, atmos_subset)) call stop_on_err(sources%get_subset(colS, ncol/2, sources_subset)) - call stop_on_err(rte_lw(atmos_subset, top_at_1, & - sources_subset, & + call stop_on_err(rte_lw(atmos_subset, & + sources_subset, & sfc_emis(:,colS:colE), & fluxes)) flux_up(colS:colE,:) = up diff --git a/tests/rte_sw_solver_unit_tests.F90 b/tests/rte_sw_solver_unit_tests.F90 index 196fd1f4a..610e20d3a 100644 --- a/tests/rte_sw_solver_unit_tests.F90 +++ b/tests/rte_sw_solver_unit_tests.F90 @@ -69,7 +69,7 @@ program rte_sw_solver_unit_tests type(ty_optical_props_2str) :: atmos type(ty_fluxes_broadband) :: fluxes - logical :: top_at_1 + logical, parameter :: top_at_1 = .true. real(wp), dimension(ncol,nlay+1), target :: & ref_flux_up, ref_flux_dn, ref_flux_dir, ref_flux_net, & tst_flux_up, tst_flux_dn, tst_flux_dir, tst_flux_net @@ -81,8 +81,6 @@ program rte_sw_solver_unit_tests logical :: passed - ! ------------------------------------------------------------------------------------------------------ - top_at_1 = .true. ! ------------------------------------------------------------------------------------ ! ! Shortwave tests - thin atmospheres @@ -95,6 +93,7 @@ program rte_sw_solver_unit_tests band_lims_gpt = reshape([1, 1 ], shape = [2, 1]), & name = "Gray atmosphere")) call stop_on_err(atmos%alloc_2str(ncol, nlay)) + call atmos%set_top_at_1(top_at_1) call thin_scattering(tau, ssa, g, nlay, atmos) do imu0 = 1, nmu0 @@ -104,7 +103,7 @@ program rte_sw_solver_unit_tests fluxes%flux_dn => ref_flux_dn (:,:) fluxes%flux_dn_dir => ref_flux_dir(:,:) fluxes%flux_net => ref_flux_net(:,:) - call stop_on_err(rte_sw(atmos, top_at_1, & + call stop_on_err(rte_sw(atmos, & mu0_arr, & toa_flux, & sfc_albedo, sfc_albedo, & @@ -120,7 +119,7 @@ program rte_sw_solver_unit_tests ! ! Check direct beam for correctness with Beer-Lambert-Bouguier ! - if(top_at_1) then + if(atmos%top_is_at_1()) then sfc => ref_flux_dir(:,nlay+1) else sfc => ref_flux_dir(:, 1) @@ -143,7 +142,7 @@ program rte_sw_solver_unit_tests ! nullify(fluxes%flux_up) nullify(fluxes%flux_dn) - call stop_on_err(rte_sw(atmos, top_at_1, & + call stop_on_err(rte_sw(atmos, & mu0_arr, & toa_flux, & sfc_albedo, sfc_albedo, & @@ -155,7 +154,7 @@ program rte_sw_solver_unit_tests ! fluxes%flux_up => tst_flux_up (:,:) fluxes%flux_dn => tst_flux_dn (:,:) - call stop_on_err(rte_sw(atmos, top_at_1, & + call stop_on_err(rte_sw(atmos, & mu0_arr, & toa_flux, & sfc_albedo, sfc_albedo, & @@ -178,7 +177,7 @@ program rte_sw_solver_unit_tests print *, " Vertical orientation invariance" call thin_scattering(tau, ssa, g, nlay, atmos) call vr(atmos) - call stop_on_err(rte_sw(atmos, .not. top_at_1, & + call stop_on_err(rte_sw(atmos, & mu0_arr, & toa_flux, & sfc_albedo, sfc_albedo, & @@ -201,7 +200,7 @@ program rte_sw_solver_unit_tests ! print *, " Linear in TOA flux" call thin_scattering(tau, ssa, g, nlay, atmos) - call stop_on_err(rte_sw(atmos, top_at_1, & + call stop_on_err(rte_sw(atmos, & mu0_arr, & toa_flux * factor, & sfc_albedo, sfc_albedo, & @@ -332,7 +331,7 @@ subroutine clear_sky_subset(atmos, mu0, toa_flux, sfc_albedo, flux_up, flux_dn) colS = ((i-1) * ncol/2) + 1 colE = i * ncol/2 call stop_on_err(atmos%get_subset(colS, ncol/2, atmos_subset)) - call stop_on_err(rte_sw(atmos_subset, top_at_1, & + call stop_on_err(rte_sw(atmos_subset, & mu0(colS:colE), & toa_flux(colS:colE,:), & sfc_albedo(:,colS:colE), sfc_albedo(:,colS:colE), & diff --git a/tests/test_zenith_angle_spherical_correction.F90 b/tests/test_zenith_angle_spherical_correction.F90 index 26ad615e0..917349441 100644 --- a/tests/test_zenith_angle_spherical_correction.F90 +++ b/tests/test_zenith_angle_spherical_correction.F90 @@ -38,7 +38,6 @@ program test_solar_zenith_angle p_lev(ncol, nlay+1) real(wp), dimension(ncol, nlay+1), target & :: flux_up, flux_dn, flux_dir - logical :: top_at_1 character(len=128) :: k_dist_file = "rrtmgp-gas-sw-g112.nc" real(wp), dimension(:,:), allocatable & :: toa_flux, sfc_alb_dir, sfc_alb_dif @@ -95,9 +94,7 @@ program test_solar_zenith_angle ! Set small solar zenith angles, varying by column; compute fluxes and heating rates ! call stop_on_err(zenith_angle_with_height(z_lay(:,1), 3 * [0.01_wp, 0.02_wp, 0.03_wp, 0.04_wp], z_lay, mu0)) - top_at_1 = p_lay(1, 1) < p_lay(1,nlay) - call stop_on_err(rte_sw(atmos, top_at_1, & - mu0, toa_flux, & + call stop_on_err(rte_sw(atmos, mu0, toa_flux, & sfc_alb_dir, sfc_alb_dif, & fluxes)) call stop_on_err(compute_heating_rate(flux_up, flux_dn, flux_dir, p_lev, mu0, heating_rate)) diff --git a/tests/validation-plots.py b/tests/validation-plots.py old mode 100644 new mode 100755 index 4dcc9434b..9564e9216 --- a/tests/validation-plots.py +++ b/tests/validation-plots.py @@ -1,3 +1,5 @@ +#! /usr/bin/env python + import colorcet as cc import matplotlib as mpl import matplotlib.pyplot as plt @@ -7,6 +9,7 @@ import xarray as xr from matplotlib.backends.backend_pdf import PdfPages +import argparse def mae(diff, col_dim): # @@ -62,6 +65,36 @@ def construct_lbl_esgf_root(var, esgf_node="llnl"): ######################################################################## def main(): + parser = argparse.ArgumentParser() + parser.add_argument( + "--state_file", + help="Path to the state information NetCDF file.", + default="multiple_input4MIPs_radiation_RFMIP_UColorado-RFMIP-1-2_none.nc" + ) + parser.add_argument( + "--lw_vars_file", + help="Path to the LW results file", + default="lw_flux_variants.nc" + + ) + parser.add_argument( + "--sw_vars_file", + help="Path to the SW results file", + default="sw_flux_variants.nc" + + ) + parser.add_argument( + "--output_pdf", + help="Path to the output PDF file for validation plots.", + default="validation-figures.pdf" + ) + args = parser.parse_args() + + state_file = args.state_file + lw_vars_file = args.lw_vars_file + sw_vars_file = args.sw_vars_file + output_pdf = args.output_pdf + warnings.simplefilter("ignore", xr.SerializationWarning) # # Reference values from LBLRTM - download locally, since OpenDAP access is @@ -87,7 +120,7 @@ def main(): # # Open the test results # - gp = xr.open_dataset("test_atmospheres.nc") + gp = xr.open_mfdataset([state_file, lw_vars_file, sw_vars_file]) # # Does the flux plus the Jacobian equal a calculation with perturbed surface # temperature? @@ -120,7 +153,7 @@ def main(): plev.load() gpi.load() lbli.load() - with PdfPages('validation-figures.pdf') as pdf: + with PdfPages(output_pdf) as pdf: ######################################################################## # Longwave ########################################################################