From 6616f4f797af5b719f5be933313e0591d40e4932 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Thu, 14 Dec 2023 21:10:14 +0000 Subject: [PATCH 01/20] Bump actions/upload-artifact from 3 to 4 Bumps [actions/upload-artifact](https://github.com/actions/upload-artifact) from 3 to 4. - [Release notes](https://github.com/actions/upload-artifact/releases) - [Commits](https://github.com/actions/upload-artifact/compare/v3...v4) --- updated-dependencies: - dependency-name: actions/upload-artifact dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] --- .github/workflows/containerized-ci.yml | 2 +- .github/workflows/doc-deployment.yml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/containerized-ci.yml b/.github/workflows/containerized-ci.yml index ba5382a57..e4167a473 100644 --- a/.github/workflows/containerized-ci.yml +++ b/.github/workflows/containerized-ci.yml @@ -110,7 +110,7 @@ jobs: # - name: Upload validation plots if: matrix.fortran-compiler == 'ifort' && matrix.rte-kernels == 'default' && matrix.fpmodel == 'DP' - uses: actions/upload-artifact@v3 + uses: actions/upload-artifact@v4 with: name: valdiation-plot path: tests/validation-figures.pdf diff --git a/.github/workflows/doc-deployment.yml b/.github/workflows/doc-deployment.yml index 1d67ab9e8..2ee1de23b 100644 --- a/.github/workflows/doc-deployment.yml +++ b/.github/workflows/doc-deployment.yml @@ -51,7 +51,7 @@ jobs: # Upload documentation # - name: Upload Documentation - uses: actions/upload-artifact@v3 + uses: actions/upload-artifact@v4 with: name: documentation path: public/ From 5ac4d0a8e9fe4e8a2769634f92e3d4e4c35f4e98 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Mon, 18 Dec 2023 11:28:21 -0500 Subject: [PATCH 02/20] Temperature and source function refinements (#253) Remove a workaround for PGI Fortran 19 and `present` status; interpolate level temperatures on device. --- rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 | 65 +++++++++++------------- 1 file changed, 31 insertions(+), 34 deletions(-) diff --git a/rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 b/rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 index 7583fe579..fcf414455 100644 --- a/rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 +++ b/rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 @@ -304,30 +304,21 @@ function gas_optics_int(this, & ! ! Interpolate source function - ! - if(present(tlev)) then - ! - ! present status of optional argument should be passed to source() - ! but isn't with PGI 19.10 - ! - error_msg = source(this, & - ncol, nlay, nband, ngpt, & - play, plev, tlay, tsfc, & - jtemp, jpress, jeta, tropo, fmajor, & - sources, & - tlev) - !$acc exit data delete(tlev) + ! present status of optional argument is passed to source() + ! + error_msg = source(this, & + ncol, nlay, nband, ngpt, & + play, plev, tlay, tsfc, & + jtemp, jpress, jeta, tropo, fmajor, & + sources, & + tlev) + if(present(tlev)) then + !$acc exit data delete(tlev) !$omp target exit data map(release:tlev) - else - error_msg = source(this, & - ncol, nlay, nband, ngpt, & - play, plev, tlay, tsfc, & - jtemp, jpress, jeta, tropo, fmajor, & - sources) - end if - !$acc exit data delete(tsfc) + end if + !$acc exit data delete(tsfc) !$omp target exit data map(release:tsfc) - !$acc exit data delete(jtemp, jpress, tropo, fmajor, jeta) + !$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 !------------------------------------------------------------------------------------------ @@ -858,7 +849,15 @@ function source(this, & error_msg = "" ! ! Source function needs temperature at interfaces/levels and at layer centers + ! Allocate small local array for tlev unconditionally ! + !$acc data copyin(sources) copyout( sources%lay_source, sources%lev_source) & + !$acc copyout( sources%sfc_source, sources%sfc_source_Jac) & + !$acc create(tlev_arr) + !$omp target data map(from:sources%lay_source, sources%lev_source) & + !$omp map(from:sources%sfc_source, sources%sfc_source_Jac) & + !$omp map(alloc:tlev_arr) + if (present(tlev)) then ! Users might have provided these tlev_wk => tlev @@ -868,32 +867,30 @@ function source(this, & ! Interpolate temperature to levels if not provided ! Interpolation and extrapolation at boundaries is weighted by pressure ! + !$acc parallel loop gang vector + !$omp target teams distribute parallel do simd do icol = 1, ncol - tlev_arr(icol,1) = tlay(icol,1) & + tlev_arr(icol,1) = tlay(icol,1) & + (plev(icol,1)-play(icol,1))*(tlay(icol,2)-tlay(icol,1)) & - & / (play(icol,2)-play(icol,1)) + / (play(icol,2)-play(icol,1)) + tlev_arr(icol,nlay+1) = tlay(icol,nlay) & + + (plev(icol,nlay+1)-play(icol,nlay))*(tlay(icol,nlay)-tlay(icol,nlay-1)) & + / (play(icol,nlay)-play(icol,nlay-1)) end do - do ilay = 2, nlay + !$acc parallel loop gang vector collapse(2) + !$omp target teams distribute parallel do simd collapse(2) + do ilay = 2, nlay do icol = 1, ncol tlev_arr(icol,ilay) = (play(icol,ilay-1)*tlay(icol,ilay-1)*(plev(icol,ilay )-play(icol,ilay)) & + play(icol,ilay )*tlay(icol,ilay )*(play(icol,ilay-1)-plev(icol,ilay))) / & (plev(icol,ilay)*(play(icol,ilay-1) - play(icol,ilay))) end do end do - do icol = 1, ncol - tlev_arr(icol,nlay+1) = tlay(icol,nlay) & - + (plev(icol,nlay+1)-play(icol,nlay))*(tlay(icol,nlay)-tlay(icol,nlay-1)) & - / (play(icol,nlay)-play(icol,nlay-1)) - end do end if !------------------------------------------------------------------- ! Compute internal (Planck) source functions at layers and levels, ! which depend on mapping from spectral space that creates k-distribution. - !$acc data copyin(sources) copyout( sources%lay_source, sources%lev_source) & - !$acc copyout( sources%sfc_source, sources%sfc_source_Jac) - !$omp target data map(from:sources%lay_source, sources%lev_source) & - !$omp map(from:sources%sfc_source, sources%sfc_source_Jac) !$acc kernels copyout(top_at_1) !$omp target map(from:top_at_1) From 5c55f2c0bdaee8b2c6f0cc53ebfd6cc97d8feadc Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Wed, 17 Jan 2024 15:30:13 -0500 Subject: [PATCH 03/20] Bump actions/cache from 3 to 4 (#254) Bumps [actions/cache](https://github.com/actions/cache) from 3 to 4. - [Release notes](https://github.com/actions/cache/releases) - [Changelog](https://github.com/actions/cache/blob/main/RELEASES.md) - [Commits](https://github.com/actions/cache/compare/v3...v4) --- updated-dependencies: - dependency-name: actions/cache dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/continuous-integration.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/continuous-integration.yml b/.github/workflows/continuous-integration.yml index 31d6f7328..3e19910a2 100644 --- a/.github/workflows/continuous-integration.yml +++ b/.github/workflows/continuous-integration.yml @@ -64,7 +64,7 @@ jobs: # Cache Conda packages # - name: Cache Conda packages - uses: actions/cache@v3 + uses: actions/cache@v4 with: path: ~/conda_pkgs_dir key: conda-pkgs From 44e4352d295c285214ca640e35f00193d3eddeff Mon Sep 17 00:00:00 2001 From: Sergey Kosukhin Date: Wed, 24 Jan 2024 18:23:05 +0100 Subject: [PATCH 04/20] Run CI jobs on Levante via GitLab (#255) Add CI jobs that run the tests on the GPU partition of Levante. The compiler version and flags are what we currently use to build ICON by default. For security reasons, the jobs are skipped for PRs from the forks. --- .github/workflows/gitlab-ci.yml | 89 +++++++++++++++++++++++++++++++++ .gitlab/levante.yml | 73 +++++++++++++++++++++++++++ 2 files changed, 162 insertions(+) create mode 100644 .github/workflows/gitlab-ci.yml create mode 100644 .gitlab/levante.yml diff --git a/.github/workflows/gitlab-ci.yml b/.github/workflows/gitlab-ci.yml new file mode 100644 index 000000000..0567d97b5 --- /dev/null +++ b/.github/workflows/gitlab-ci.yml @@ -0,0 +1,89 @@ +name: GitLab CI +on: + push: + branches-ignore: + - documentation + pull_request: + branches-ignore: + - documentation + +defaults: + run: + shell: bash + +jobs: + # + # Deferred GitLab pipelines on Levante at DKRZ (see .gitlab/levante.yml): + # + levante-init: + if: | + github.repository_owner == 'earth-system-radiation' && + ( github.event_name != 'pull_request' || + github.event.pull_request.head.repo.owner.login == github.repository_owner ) + runs-on: ubuntu-latest + outputs: + ref-name: ${{ steps.g-push-rev.outputs.ref-name }} + pipeline-id: ${{ steps.gl-trigger-pipeline.outputs.pipeline-id }} + steps: + # + # Check out GitHub repository + # + - name: Check out GitHub repository + uses: actions/checkout@v3 + with: + fetch-depth: 0 + # + # Push to GitLab repository + # + - name: Push to GitLab repository + id: g-push-rev + uses: "skosukhin/git-ci-hub-lab/g-push-rev@v1" + with: + remote-url: ${{ vars.DKRZ_GITLAB_SERVER }}/${{ vars.DKRZ_GITLAB_PROJECT }}.git + password: ${{ secrets.DKRZ_GITLAB_TOKEN }} + ref-type: tag + ref-message: ${{ github.server_url }}/${{ github.repository }}/actions/runs/${{ github.run_id }} + force-push: true + # + # Trigger GitLab CI/CD Pipeline + # + - name: Trigger GitLab CI/CD Pipeline + id: gl-trigger-pipeline + uses: "skosukhin/git-ci-hub-lab/gl-trigger-pipeline@v1" + with: + server-url: ${{ vars.DKRZ_GITLAB_SERVER }} + project-name: ${{ vars.DKRZ_GITLAB_PROJECT }} + token: ${{ secrets.DKRZ_GITLAB_TRIGGER_TOKEN }} + ref-name: ${{ steps.g-push-rev.outputs.ref-name }} + expected-sha: ${{ github.sha }} + levante: + runs-on: ubuntu-latest + needs: levante-init + strategy: + fail-fast: false + matrix: + config-name: [nvhpc-gpu-openacc-DP, nvhpc-gpu-openacc-SP] + steps: + # + # Build, run and check (fetch the log) + # + - name: Build, run and check (fetch the log) + uses: "skosukhin/git-ci-hub-lab/gl-attach-job@v1" + with: + server-url: ${{ vars.DKRZ_GITLAB_SERVER }} + project-name: ${{ vars.DKRZ_GITLAB_PROJECT }} + token: ${{ secrets.DKRZ_GITLAB_TOKEN }} + pipeline-id: ${{ needs.levante-init.outputs.pipeline-id }} + job-name: ${{ matrix.config-name }} + levante-cleanup: + runs-on: ubuntu-latest + needs: [levante-init, levante] + if: always() && needs.levante-init.result != 'skipped' + continue-on-error: true + steps: + - uses: "skosukhin/git-ci-hub-lab/g-delete-ref@v1" + with: + remote-url: ${{ vars.DKRZ_GITLAB_SERVER }}/${{ vars.DKRZ_GITLAB_PROJECT }}.git + password: ${{ secrets.DKRZ_GITLAB_TOKEN }} + ref-type: tag + ref-name: ${{ needs.levante-init.outputs.ref-name }} diff --git a/.gitlab/levante.yml b/.gitlab/levante.yml new file mode 100644 index 000000000..93d01fa87 --- /dev/null +++ b/.gitlab/levante.yml @@ -0,0 +1,73 @@ +workflow: + rules: + - if: $CI_PIPELINE_SOURCE == "trigger" + +include: + - project: 'anw_dienste/ci-templates' + file: '.slurm-ci.yml' + +variables: + SCHEDULER_PARAMETERS: >- + --account=mh0287 + --partition=gpu + --gpus=1 + --time=05:00 + +.build-common: + extends: .default + variables: + # Core variables: + FC: /sw/spack-levante/nvhpc-22.5-v4oky3/Linux_x86_64/22.5/compilers/bin/nvfortran + # Production flags for ICON model: + FCFLAGS: -g -O2 -Mrecursive -Mallocatable=03 -Mstack_arrays -Minfo=accel,inline -acc=gpu,verystrict -gpu=cc80,cuda11.7 -DRTE_USE_${FPMODEL} + # Convenience variables: + NFHOME: /sw/spack-levante/netcdf-fortran-4.5.4-syv4qr + NCHOME: /sw/spack-levante/netcdf-c-4.9.0-gc7kgj + 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 + RRTMGP_ROOT: ${CI_PROJECT_DIR} + RRTMGP_DATA: ${CI_PROJECT_DIR}/rrtmgp-data + RTE_KERNELS: accel + 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-}" + # The -Mstack_arrays compiler flag requires a large stack: + - ulimit -s unlimited + script: + # + # Build libraries, examples and tests + # + - ${FC} --version + - make libs + - make -C build separate-libs + # + # Check out data + # + - git clone --depth 1 https://github.com/earth-system-radiation/rrtmgp-data.git "${RRTMGP_DATA}" + # + # Run examples and tests + # + - make tests + # + # Compare the results + # + - make check + +nvhpc-gpu-openacc-DP: + extends: .build-common + variables: + FPMODEL: DP + FAILURE_THRESHOLD: "5.8e-2" + +nvhpc-gpu-openacc-SP: + extends: .build-common + variables: + FPMODEL: SP + FAILURE_THRESHOLD: "3.5e-1" From 592c30a072b977f032d32e828188f4213ed64c5c Mon Sep 17 00:00:00 2001 From: Sergey Kosukhin Date: Fri, 26 Jan 2024 15:14:05 +0100 Subject: [PATCH 05/20] Extend Levante CI with NAG (#259) This adds CI tests with the NAG compiler on Levante with compiler version and flags currently used to build ICON by default. Both DP and SP floating models are tested, as are default and accelerator kernels, but not all possible combinations. Accelerator kernels fail with a run-time error and are marked as experimental. --- .github/workflows/gitlab-ci.yml | 19 +++++- .gitlab/levante.yml | 113 +++++++++++++++++++++++++++----- 2 files changed, 113 insertions(+), 19 deletions(-) diff --git a/.github/workflows/gitlab-ci.yml b/.github/workflows/gitlab-ci.yml index 0567d97b5..941e85796 100644 --- a/.github/workflows/gitlab-ci.yml +++ b/.github/workflows/gitlab-ci.yml @@ -29,7 +29,7 @@ jobs: # Check out GitHub repository # - name: Check out GitHub repository - uses: actions/checkout@v3 + uses: actions/checkout@v4 with: fetch-depth: 0 # @@ -59,10 +59,24 @@ jobs: levante: runs-on: ubuntu-latest needs: levante-init + continue-on-error: ${{ matrix.experimental }} strategy: fail-fast: false matrix: - config-name: [nvhpc-gpu-openacc-DP, nvhpc-gpu-openacc-SP] + config-name: + - nvhpc-gpu-openacc-DP + - nvhpc-gpu-openacc-SP + #- nag-cpu-default-DP + - nag-cpu-default-SP + - nag-cpu-accel-DP + #- nag-cpu-accel-SP + include: + # The tests are not experimental by default: + - experimental: false + - config-name: nag-cpu-accel-DP + experimental: true + #- config-name: nag-cpu-accel-SP + # experimental: true steps: # # Build, run and check (fetch the log) @@ -87,3 +101,4 @@ jobs: password: ${{ secrets.DKRZ_GITLAB_TOKEN }} ref-type: tag ref-name: ${{ needs.levante-init.outputs.ref-name }} + force: true diff --git a/.gitlab/levante.yml b/.gitlab/levante.yml index 93d01fa87..8104acd52 100644 --- a/.gitlab/levante.yml +++ b/.gitlab/levante.yml @@ -9,42 +9,74 @@ include: variables: SCHEDULER_PARAMETERS: >- --account=mh0287 - --partition=gpu - --gpus=1 --time=05:00 + ${EXTRA_SCHEDULER_PARAMETERS} + EXTRA_SCHEDULER_PARAMETERS: -.build-common: +.gpu: extends: .default + variables: + EXTRA_SCHEDULER_PARAMETERS: >- + --partition=gpu + --gpus=1 + +.cpu: + extends: .default + variables: + EXTRA_SCHEDULER_PARAMETERS: >- + --partition=shared + +.nvhpc: variables: # Core variables: FC: /sw/spack-levante/nvhpc-22.5-v4oky3/Linux_x86_64/22.5/compilers/bin/nvfortran - # Production flags for ICON model: - FCFLAGS: -g -O2 -Mrecursive -Mallocatable=03 -Mstack_arrays -Minfo=accel,inline -acc=gpu,verystrict -gpu=cc80,cuda11.7 -DRTE_USE_${FPMODEL} # 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 + +.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 + +.dp: + variables: + FPMODEL: DP + FAILURE_THRESHOLD: "5.8e-2" + +.sp: + variables: + FPMODEL: SP + FAILURE_THRESHOLD: "3.5e-1" + +.common: + variables: PYHOME: /sw/spack-levante/mambaforge-22.9.0-2-Linux-x86_64-kptncg - # Suppress an irrelevant but annoying error message: + # 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 RRTMGP_ROOT: ${CI_PROJECT_DIR} RRTMGP_DATA: ${CI_PROJECT_DIR}/rrtmgp-data - RTE_KERNELS: accel 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-}" - # The -Mstack_arrays compiler flag requires a large stack: + # Some tests require a large stack: - ulimit -s unlimited script: # # Build libraries, examples and tests # - - ${FC} --version + - ${FC} ${VERSION_FCFLAGS} - make libs - make -C build separate-libs # @@ -60,14 +92,61 @@ variables: # - make check -nvhpc-gpu-openacc-DP: - extends: .build-common +.nvhpc-gpu-openacc: + extends: + - .gpu + - .nvhpc + - .common variables: - FPMODEL: DP - FAILURE_THRESHOLD: "5.8e-2" + # 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 -nvhpc-gpu-openacc-SP: - extends: .build-common +.nag-cpu: + extends: + - .cpu + - .nag + - .common variables: - FPMODEL: SP - FAILURE_THRESHOLD: "3.5e-1" + # 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} + +.nag-cpu-default: + extends: .nag-cpu + variables: + RTE_KERNELS: default + +.nag-cpu-accel: + extends: .nag-cpu + variables: + RTE_KERNELS: accel + +nvhpc-gpu-openacc-DP: + extends: + - .dp + - .nvhpc-gpu-openacc + +nvhpc-gpu-openacc-SP: + extends: + - .sp + - .nvhpc-gpu-openacc + +#nag-cpu-default-DP: +# extends: +# - .dp +# - .nag-cpu-default + +nag-cpu-default-SP: + extends: + - .sp + - .nag-cpu-default + +nag-cpu-accel-DP: + extends: + - .dp + - .nag-cpu-accel + +#nag-cpu-accel-SP: +# extends: +# - .sp +# - .nag-cpu-accel From 0a17084bf89ef061713789d875b94aede22624c3 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Fri, 26 Jan 2024 13:21:28 -0500 Subject: [PATCH 06/20] Adjust tolerances, fix error in accelerator source function kernel (#258) Scripts fails when individual tests fail, small tolerance adjustments to tests pass, small correction to LW source calculation in accel kernels, revert change in 5ac4d0a: nvfortran fails to pass present status of args to lower-level procedures --- examples/all-sky/all_tests.sh | 1 + rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 | 24 ++++++++++++------ .../accel/mo_gas_optics_rrtmgp_kernels.F90 | 14 ++++++----- tests/all_tests.sh | 1 + tests/check_equivalence.F90 | 25 ++++++++++--------- 5 files changed, 39 insertions(+), 26 deletions(-) diff --git a/examples/all-sky/all_tests.sh b/examples/all-sky/all_tests.sh index bbc272ef6..61a2bceb7 100644 --- a/examples/all-sky/all_tests.sh +++ b/examples/all-sky/all_tests.sh @@ -1,3 +1,4 @@ +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 \ diff --git a/rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 b/rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 index fcf414455..9aeb476f5 100644 --- a/rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 +++ b/rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 @@ -304,17 +304,25 @@ function gas_optics_int(this, & ! ! Interpolate source function - ! present status of optional argument is passed to source() - ! - error_msg = source(this, & - ncol, nlay, nband, ngpt, & - play, plev, tlay, tsfc, & - jtemp, jpress, jeta, tropo, fmajor, & - sources, & - tlev) + ! present status of optional argument should be passed to source() + ! but nvfortran (and PGI Fortran before it) do not do so + ! if(present(tlev)) then + error_msg = source(this, & + ncol, nlay, nband, ngpt, & + play, plev, tlay, tsfc, & + jtemp, jpress, jeta, tropo, fmajor, & + sources, & + tlev) !$acc exit data delete(tlev) !$omp target exit data map(release:tlev) + else + error_msg = source(this, & + ncol, nlay, nband, ngpt, & + play, plev, tlay, tsfc, & + jtemp, jpress, jeta, tropo, fmajor, & + sources) + end if !$acc exit data delete(tsfc) !$omp target exit data map(release:tsfc) diff --git a/rrtmgp-kernels/accel/mo_gas_optics_rrtmgp_kernels.F90 b/rrtmgp-kernels/accel/mo_gas_optics_rrtmgp_kernels.F90 index 44ed0a08c..103850776 100644 --- a/rrtmgp-kernels/accel/mo_gas_optics_rrtmgp_kernels.F90 +++ b/rrtmgp-kernels/accel/mo_gas_optics_rrtmgp_kernels.F90 @@ -624,7 +624,7 @@ subroutine compute_Planck_source( & ibnd = gpoint_bands(igpt) ! itropo = 1 lower atmosphere; itropo = 2 upper atmosphere itropo = merge(1,2,tropo(icol,ilay)) !WS moved itropo inside loop for GPU - iflav = gpoint_flavor(itropo, igpt) !eta interpolation depends on band's flavor + iflav = gpoint_flavor(itropo, igpt) !eta interpolation depends on band's flavor ! interpolation in temperature, pressure, and eta pfrac = & interpolate3D(one, fmajor(:,:,:,icol,ilay,iflav), pfracin, & @@ -634,20 +634,22 @@ subroutine compute_Planck_source( & planck_function_1 = interpolate1D(tlay(icol,ilay), temp_ref_min, totplnk_delta, totplnk(:,ibnd)) lay_src (icol,ilay,igpt) = pfrac * planck_function_1 - ! Compute layer source irradiance for g-point, equals band irradiance x fraction for g-point + ! Compute level source irradiance for g-point planck_function_1 = interpolate1D(tlev(icol,ilay), temp_ref_min, totplnk_delta, totplnk(:,ibnd)) if (ilay == 1) then lev_src(icol,ilay, igpt) = pfrac * planck_function_1 - else if (ilay == nlay) then - lev_src(icol,ilay, igpt) = pfrac * planck_function_1 - planck_function_2 = interpolate1D(tlev(icol,nlay+1), temp_ref_min, totplnk_delta, totplnk(:,ibnd)) - lev_src(icol,nlay+1,igpt) = pfrac * planck_function_2 else + itropo = merge(1,2,tropo(icol,ilay-1)) !WS moved itropo inside loop for GPU + iflav = gpoint_flavor(itropo, igpt) !eta interpolation depends on band's flavor pfrac_m1 = & interpolate3D(one, fmajor(:,:,:,icol,ilay-1,iflav), pfracin, & igpt, jeta(:,icol,ilay-1,iflav), jtemp(icol,ilay-1),jpress(icol,ilay-1)+itropo) lev_src(icol,ilay, igpt) = sqrt(pfrac * pfrac_m1) * planck_function_1 end if + if (ilay == nlay) then + planck_function_1 = interpolate1D(tlev(icol,nlay+1), temp_ref_min, totplnk_delta, totplnk(:,ibnd)) + lev_src(icol,nlay+1,igpt) = pfrac * planck_function_1 + end if if (ilay == sfc_lay) then planck_function_1 = interpolate1D(tsfc(icol) , temp_ref_min, totplnk_delta, totplnk(:,ibnd)) planck_function_2 = interpolate1D(tsfc(icol) + delta_Tsurf, temp_ref_min, totplnk_delta, totplnk(:,ibnd)) diff --git a/tests/all_tests.sh b/tests/all_tests.sh index dba254b6d..9af6c7021 100644 --- a/tests/all_tests.sh +++ b/tests/all_tests.sh @@ -1,3 +1,4 @@ +set -eux ./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 diff --git a/tests/check_equivalence.F90 b/tests/check_equivalence.F90 index 8ac814965..40270178e 100644 --- a/tests/check_equivalence.F90 +++ b/tests/check_equivalence.F90 @@ -275,8 +275,8 @@ program rte_check_equivalence ! Orientation invariance ! call lw_clear_sky_vr - if(.not. allclose(tst_flux_up, ref_flux_up) .or. & - .not. allclose(tst_flux_dn, ref_flux_dn) ) & + if(.not. allclose(tst_flux_up, ref_flux_up, tol=4._wp) .or. & + .not. allclose(tst_flux_dn, ref_flux_dn, tol=4._wp) ) & call report_err(" Vertical invariance failure") print *, " Vertical orientation invariance" ! ------------------------------------------------------- @@ -411,9 +411,9 @@ program rte_check_equivalence ! Threshold of 4x spacing() works on CPUs but 8x is needed for GPUs ! call sw_clear_sky_tsi - if(.not. allclose(tst_flux_up, ref_flux_up, tol = 8._wp) .or. & - .not. allclose(tst_flux_dn, ref_flux_dn, tol = 8._wp) .or. & - .not. allclose(tst_flux_dir,ref_flux_dir,tol = 8._wp)) & + if(.not. allclose(tst_flux_up, ref_flux_up, tol = 10._wp) .or. & + .not. allclose(tst_flux_dn, ref_flux_dn, tol = 8._wp) .or. & + .not. allclose(tst_flux_dir,ref_flux_dir,tol = 8._wp)) & call report_err(" Changing TSI fails") print *, " TSI invariance" ! ------------------------------------------------------- @@ -432,9 +432,9 @@ program rte_check_equivalence mu0, toa_flux, & sfc_alb_dir, sfc_alb_dif, & fluxes)) - if(.not. allclose(tst_flux_up, ref_flux_up, tol = 6._wp) .or. & + if(.not. allclose(tst_flux_up, ref_flux_up, tol = 8._wp) .or. & .not. allclose(tst_flux_dn, ref_flux_dn, tol = 6._wp) .or. & - .not. allclose(tst_flux_dir,ref_flux_dir,tol = 6._wp)) & + .not. allclose(tst_flux_dir,ref_flux_dir,tol = 8._wp)) & call report_err(" halving/doubling fails") call increment_with_1scl @@ -442,7 +442,7 @@ program rte_check_equivalence mu0, toa_flux, & sfc_alb_dir, sfc_alb_dif, & fluxes)) - if(.not. allclose(tst_flux_up, ref_flux_up, tol = 6._wp) .or. & + if(.not. allclose(tst_flux_up, ref_flux_up, tol = 8._wp) .or. & .not. allclose(tst_flux_dn, ref_flux_dn, tol = 6._wp) .or. & .not. allclose(tst_flux_dir,ref_flux_dir,tol = 6._wp)) & call report_err(" Incrementing with 1scl fails") @@ -453,9 +453,9 @@ program rte_check_equivalence atmos, & toa_flux)) call increment_with_2str - if(.not. allclose(tst_flux_up, ref_flux_up, tol = 6._wp) .or. & - .not. allclose(tst_flux_dn, ref_flux_dn, tol = 6._wp) .or. & - .not. allclose(tst_flux_dir,ref_flux_dir,tol = 6._wp)) & + if(.not. allclose(tst_flux_up, ref_flux_up, tol = 8._wp) .or. & + .not. allclose(tst_flux_dn, ref_flux_dn, tol = 6._wp) .or. & + .not. allclose(tst_flux_dir,ref_flux_dir,tol = 6._wp)) & call report_err(" Incrementing with 2str fails") call stop_on_err(k_dist%gas_optics(p_lay, p_lev, & @@ -468,7 +468,7 @@ program rte_check_equivalence mu0, toa_flux, & sfc_alb_dir, sfc_alb_dif, & fluxes)) - if(.not. allclose(tst_flux_up, ref_flux_up, tol = 6._wp) .or. & + if(.not. allclose(tst_flux_up, ref_flux_up, tol = 8._wp) .or. & .not. allclose(tst_flux_dn, ref_flux_dn, tol = 6._wp) .or. & .not. allclose(tst_flux_dir,ref_flux_dir,tol = 6._wp)) & call report_err(" Incrementing with nstr fails") @@ -489,6 +489,7 @@ subroutine lw_clear_sky_vr character(32), & dimension(gas_concs%get_num_gases()) & :: gc_gas_names + ! ! Reverse the orientation of the problem ! From 9dcc9ba6346887710b74ce7d2b2255c994333762 Mon Sep 17 00:00:00 2001 From: Sergey Kosukhin Date: Tue, 30 Jan 2024 19:03:41 +0100 Subject: [PATCH 07/20] Fix intent in subroutine lw_transport_noscat_up (#261) Fix subroutine intent, mark corresponding CI as not-experimental --- .github/workflows/gitlab-ci.yml | 4 ---- rte-kernels/accel/mo_rte_solver_kernels.F90 | 2 +- 2 files changed, 1 insertion(+), 5 deletions(-) diff --git a/.github/workflows/gitlab-ci.yml b/.github/workflows/gitlab-ci.yml index 941e85796..392d4ba8f 100644 --- a/.github/workflows/gitlab-ci.yml +++ b/.github/workflows/gitlab-ci.yml @@ -73,10 +73,6 @@ jobs: include: # The tests are not experimental by default: - experimental: false - - config-name: nag-cpu-accel-DP - experimental: true - #- config-name: nag-cpu-accel-SP - # experimental: true steps: # # Build, run and check (fetch the log) diff --git a/rte-kernels/accel/mo_rte_solver_kernels.F90 b/rte-kernels/accel/mo_rte_solver_kernels.F90 index ba1175957..42d234abd 100644 --- a/rte-kernels/accel/mo_rte_solver_kernels.F90 +++ b/rte-kernels/accel/mo_rte_solver_kernels.F90 @@ -790,7 +790,7 @@ subroutine lw_transport_noscat_up(ncol, nlay, ngpt, & logical(wl), intent(in ) :: top_at_1 ! real(wp), dimension(ncol,nlay ,ngpt), intent(in ) :: trans ! transmissivity = exp(-tau) real(wp), dimension(ncol,nlay ,ngpt), intent(in ) :: source_up ! Diffuse radiation emitted by the layer - real(wp), dimension(ncol,nlay+1,ngpt), intent( out) :: radn_up ! Radiances [W/m2-str] + real(wp), dimension(ncol,nlay+1,ngpt), intent(inout) :: radn_up ! Radiances [W/m2-str] logical(wl), intent(in ) :: do_Jacobians real(wp), dimension(ncol,nlay+1,ngpt), intent(inout) :: radn_upJac ! surface temperature Jacobian of Radiances [W/m2-str / K] ! Local variables From a0ed5bfcb493c20c412f6e5185e65637afdabb67 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Fri, 9 Feb 2024 12:54:36 -0800 Subject: [PATCH 08/20] Reduce use of CI (#264) Enable manual runs of CI , don't run CI on push except for main, develop --- .github/workflows/containerized-ci.yml | 6 ++++-- .github/workflows/continuous-integration.yml | 6 ++++-- .github/workflows/gitlab-ci.yml | 6 ++++-- .github/workflows/self-hosted-ci.yml | 6 ++++-- 4 files changed, 16 insertions(+), 8 deletions(-) diff --git a/.github/workflows/containerized-ci.yml b/.github/workflows/containerized-ci.yml index e4167a473..bf48d3265 100644 --- a/.github/workflows/containerized-ci.yml +++ b/.github/workflows/containerized-ci.yml @@ -1,11 +1,13 @@ name: Continuous integration in a box on: push: - branches-ignore: - - documentation + branches: + - main + - develop pull_request: branches-ignore: - documentation + workflow_dispatch: jobs: Containerized-CI: diff --git a/.github/workflows/continuous-integration.yml b/.github/workflows/continuous-integration.yml index 3e19910a2..2c5a931bb 100644 --- a/.github/workflows/continuous-integration.yml +++ b/.github/workflows/continuous-integration.yml @@ -1,11 +1,13 @@ name: Continuous Integration on: push: - branches-ignore: - - documentation + branches: + - main + - develop pull_request: branches-ignore: - documentation + workflow_dispatch: defaults: run: diff --git a/.github/workflows/gitlab-ci.yml b/.github/workflows/gitlab-ci.yml index 392d4ba8f..fce24d1f8 100644 --- a/.github/workflows/gitlab-ci.yml +++ b/.github/workflows/gitlab-ci.yml @@ -1,11 +1,13 @@ name: GitLab CI on: push: - branches-ignore: - - documentation + branches: + - main + - develop pull_request: branches-ignore: - documentation + workflow_dispatch: defaults: run: diff --git a/.github/workflows/self-hosted-ci.yml b/.github/workflows/self-hosted-ci.yml index a5409d80a..64fc77adb 100644 --- a/.github/workflows/self-hosted-ci.yml +++ b/.github/workflows/self-hosted-ci.yml @@ -1,11 +1,13 @@ name: Self-hosted CI on: push: - branches-ignore: - - documentation + branches: + - main + - develop pull_request: branches-ignore: - documentation + workflow_dispatch: defaults: run: From ebb1b8caaa977aae5bf34334298755b9fcf5a670 Mon Sep 17 00:00:00 2001 From: Sergey Kosukhin Date: Fri, 16 Feb 2024 18:06:40 +0100 Subject: [PATCH 09/20] Run CI jobs on Lumi via GitLab (#265) --- .github/workflows/gitlab-ci.yml | 103 +++++++++++++++++++++++++ .gitlab/lumi.yml | 130 ++++++++++++++++++++++++++++++++ 2 files changed, 233 insertions(+) create mode 100644 .gitlab/lumi.yml diff --git a/.github/workflows/gitlab-ci.yml b/.github/workflows/gitlab-ci.yml index fce24d1f8..ce17ebc5c 100644 --- a/.github/workflows/gitlab-ci.yml +++ b/.github/workflows/gitlab-ci.yml @@ -100,3 +100,106 @@ jobs: ref-type: tag ref-name: ${{ needs.levante-init.outputs.ref-name }} force: true + # + # Deferred GitLab pipelines on Lumi at CSC (see .gitlab/lumi.yml): + # + lumi-init: + if: | + github.repository_owner == 'earth-system-radiation' && + ( github.event_name != 'pull_request' || + github.event.pull_request.head.repo.owner.login == github.repository_owner ) + runs-on: ubuntu-latest + outputs: + ref-name: ${{ steps.g-push-rev.outputs.ref-name }} + pipeline-id: ${{ steps.gl-create-pipeline.outputs.pipeline-id }} + steps: + # + # Check out GitHub repository + # + - name: Check out GitHub repository + uses: actions/checkout@v4 + with: + fetch-depth: 0 + # + # Push to GitLab repository + # + - name: Push to GitLab repository + id: g-push-rev + uses: "skosukhin/git-ci-hub-lab/g-push-rev@v1" + with: + remote-url: ${{ vars.GITLAB_SERVER }}/${{ vars.GITLAB_PROJECT }}.git + password: ${{ secrets.GITLAB_TOKEN }} + 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 }} + force-push: true + # + # Create GitLab CI/CD Pipeline + # + - name: Create GitLab CI/CD Pipeline + id: gl-create-pipeline + uses: "skosukhin/git-ci-hub-lab/gl-create-pipeline@v1" + with: + server-url: ${{ vars.GITLAB_SERVER }} + project-name: ${{ vars.GITLAB_PROJECT }} + token: ${{ secrets.GITLAB_TOKEN }} + ref-name: ${{ steps.g-push-rev.outputs.ref-name }} + expected-sha: ${{ steps.g-push-rev.outputs.ref-commit }} + # + # Set up Python virtual environment (fetch the log) + # + - name: Set up Python virtual environment (fetch the log) + uses: "skosukhin/git-ci-hub-lab/gl-attach-job@v1" + with: + server-url: ${{ vars.GITLAB_SERVER }} + project-name: ${{ vars.GITLAB_PROJECT }} + token: ${{ secrets.GITLAB_TOKEN }} + pipeline-id: ${{ steps.gl-create-pipeline.outputs.pipeline-id }} + job-name: setup-python + lumi: + runs-on: ubuntu-latest + needs: lumi-init + continue-on-error: ${{ matrix.experimental }} + strategy: + fail-fast: false + matrix: + config-name: + - cce-gpu-openacc-DP + - cce-gpu-openacc-SP + include: + # The tests are not experimental by default: + - experimental: false + steps: + # + # Build, run and check (fetch the log) + # + - name: Build, run and check (fetch the log) + uses: "skosukhin/git-ci-hub-lab/gl-attach-job@v1" + with: + server-url: ${{ vars.GITLAB_SERVER }} + project-name: ${{ vars.GITLAB_PROJECT }} + token: ${{ secrets.GITLAB_TOKEN }} + pipeline-id: ${{ needs.lumi-init.outputs.pipeline-id }} + job-name: ${{ matrix.config-name }} + lumi-cleanup: + runs-on: ubuntu-latest + needs: [lumi-init, lumi] + if: always() && needs.lumi-init.result != 'skipped' + continue-on-error: true + steps: + - uses: "skosukhin/git-ci-hub-lab/gl-cancel-pipeline@v1" + with: + server-url: ${{ vars.GITLAB_SERVER }} + project-name: ${{ vars.GITLAB_PROJECT }} + token: ${{ secrets.GITLAB_TOKEN }} + pipeline-id: ${{ needs.lumi-init.outputs.pipeline-id }} + - 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-name: ${{ needs.lumi-init.outputs.ref-name }} + force: true diff --git a/.gitlab/lumi.yml b/.gitlab/lumi.yml new file mode 100644 index 000000000..6899ac95d --- /dev/null +++ b/.gitlab/lumi.yml @@ -0,0 +1,130 @@ +workflow: + rules: + - if: $CI_PIPELINE_SOURCE == "api" + +default: + tags: + - lumi + +variables: + SCHEDULER_PARAMETERS: >- + --account=project_465000454 + --nodes=1 + --ntasks=1 + --cpus-per-task=4 + --mem-per-cpu=1G + --time=05:00 + ${EXTRA_SCHEDULER_PARAMETERS} + EXTRA_SCHEDULER_PARAMETERS: + +.gpu: + variables: + EXTRA_SCHEDULER_PARAMETERS: >- + --partition=dev-g + --gpus=1 + +.cpu: + variables: + EXTRA_SCHEDULER_PARAMETERS: >- + --partition=debug + +.cce: + variables: + # Core variables: + FC: ftn + # Convenience variables: + VERSION_FCFLAGS: -V + COMPILER_MODULES: PrgEnv-cray cce/16.0.1 craype-x86-milan + +.dp: + variables: + FPMODEL: DP + FAILURE_THRESHOLD: "5.8e-2" + +.sp: + variables: + FPMODEL: SP + FAILURE_THRESHOLD: "3.5e-1" + +# +# Set up Python virtual environment +# +.python-common: + variables: + PYHOME: ${CI_PROJECT_DIR}/python-venv + FF_USE_FASTZIP: 1 + +setup-python: + extends: + - .cpu + - .python-common + script: + - test ! -d "${PYHOME}" || exit 0 + - module load cray-python + - python -m venv ${PYHOME} + - ${PYHOME}/bin/python -m pip install --upgrade pip + - ${PYHOME}/bin/python -m pip install dask[array] netCDF4 numpy xarray + cache: + # Update the key to regenerate the virtual environment: + key: python-venv-version-1 + paths: + - ${PYHOME} + artifacts: + paths: + - ${PYHOME} + expire_in: 60 minutes + +.common: + extends: .python-common + needs: + - setup-python + variables: + # Make variables: + RRTMGP_ROOT: ${CI_PROJECT_DIR} + RRTMGP_DATA: ${CI_PROJECT_DIR}/rrtmgp-data + before_script: + - module --force purge + - module load ${COMPILER_MODULES} ${EXTRA_COMPILER_MODULES} cray-hdf5 cray-netcdf + # Extend the existing environment variables: + - export PATH="${PYHOME}/bin:${PATH}" + script: + # + # Build libraries, examples and tests + # + - ${FC} ${VERSION_FCFLAGS} + - make libs + - make -C build separate-libs + # + # Check out data + # + - git clone --depth 1 https://github.com/earth-system-radiation/rrtmgp-data.git "${RRTMGP_DATA}" + # + # Run examples and tests + # + - make tests + # + # Compare the results + # + - make check + +.cce-gpu-openacc: + extends: + - .gpu + - .cce + - .common + 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 + # Convenience variables: + EXTRA_COMPILER_MODULES: craype-accel-amd-gfx90a rocm + +cce-gpu-openacc-DP: + extends: + - .dp + - .cce-gpu-openacc + +cce-gpu-openacc-SP: + extends: + - .dp + - .cce-gpu-openacc From 6cda8687d7ccdfc21b299db59814486870530596 Mon Sep 17 00:00:00 2001 From: Sergey Kosukhin Date: Mon, 19 Feb 2024 17:31:32 +0100 Subject: [PATCH 10/20] Align success/failure of CI with expectations (#267) CI now reports success if a) the code compiles but fails at runtime with `cce-gpu-openmp` in self-hosted CI, and b) fails with an internal compiler error for `ifx-gpu-openmp`. See the PR for more details. --- .github/workflows/containerized-ci.yml | 33 ++++++++++++++++++++------ .github/workflows/gitlab-ci.yml | 8 ------- .github/workflows/self-hosted-ci.yml | 23 ++++++++++++------ 3 files changed, 42 insertions(+), 22 deletions(-) diff --git a/.github/workflows/containerized-ci.yml b/.github/workflows/containerized-ci.yml index bf48d3265..01343d395 100644 --- a/.github/workflows/containerized-ci.yml +++ b/.github/workflows/containerized-ci.yml @@ -12,7 +12,6 @@ on: jobs: Containerized-CI: runs-on: ubuntu-22.04 - continue-on-error: ${{ matrix.experimental }} strategy: fail-fast: false matrix: @@ -20,8 +19,6 @@ jobs: rte-kernels: [default, accel] fpmodel: [DP, SP] include: - # The tests are not experimental by default: - - experimental: false # 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 @@ -32,7 +29,6 @@ jobs: - fortran-compiler: ifx rte-kernels: accel fcflags: -debug -traceback -O0 -heap-arrays -assume realloc_lhs -extend-source 132 -stand f08 -fiopenmp -fopenmp-targets=spir64 - experimental: true # Set flags for NVIDIA Fortran compiler - fortran-compiler: nvfortran rte-kernels: default @@ -77,28 +73,51 @@ jobs: repository: earth-system-radiation/rrtmgp-data path: rrtmgp-data # - # Build libraries, examples and tests + # Build libraries, examples and tests (expect success) # - - name: Build libraries, examples and tests + - name: Build libraries, examples and tests (expect success) + id: build-success + if: matrix.fortran-compiler != 'ifx' || matrix.rte-kernels != 'accel' run: | $FC --version make libs make -C build separate-libs # + # Build libraries, examples and tests (expect failure) + # + - name: Build libraries, examples and tests (expect failure) + if: steps.build-success.outcome == 'skipped' + shell: bash + run: | + $FC --version + make 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 + } + } + # # Run examples and tests # - name: Run examples and tests + if: steps.build-success.outcome != 'skipped' run: make tests # # Relax failure thresholds for single precision # - name: Relax failure threshold for single precision - if: matrix.fpmodel == 'SP' + 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 check # # Generate validation plots diff --git a/.github/workflows/gitlab-ci.yml b/.github/workflows/gitlab-ci.yml index ce17ebc5c..0c987d844 100644 --- a/.github/workflows/gitlab-ci.yml +++ b/.github/workflows/gitlab-ci.yml @@ -61,7 +61,6 @@ jobs: levante: runs-on: ubuntu-latest needs: levante-init - continue-on-error: ${{ matrix.experimental }} strategy: fail-fast: false matrix: @@ -72,9 +71,6 @@ jobs: - nag-cpu-default-SP - nag-cpu-accel-DP #- nag-cpu-accel-SP - include: - # The tests are not experimental by default: - - experimental: false steps: # # Build, run and check (fetch the log) @@ -161,16 +157,12 @@ jobs: lumi: runs-on: ubuntu-latest needs: lumi-init - continue-on-error: ${{ matrix.experimental }} strategy: fail-fast: false matrix: config-name: - cce-gpu-openacc-DP - cce-gpu-openacc-SP - include: - # The tests are not experimental by default: - - experimental: false steps: # # Build, run and check (fetch the log) diff --git a/.github/workflows/self-hosted-ci.yml b/.github/workflows/self-hosted-ci.yml index 64fc77adb..87dd7e32d 100644 --- a/.github/workflows/self-hosted-ci.yml +++ b/.github/workflows/self-hosted-ci.yml @@ -18,15 +18,12 @@ jobs: if: github.repository == 'earth-system-radiation/rte-rrtmgp' runs-on: labels: cscs-ci - continue-on-error: ${{ matrix.experimental }} strategy: fail-fast: false matrix: config-name: [nvidia-gpu-openacc, cce-cpu-icon-production, cce-gpu-openmp] fpmodel: [DP, SP] include: - # The tests are not experimental by default: - - experimental: false - config-name: nvidia-gpu-openacc rte-kernels: accel compiler-modules: "PrgEnv-nvidia nvidia craype-accel-nvidia60 cdt-cuda/21.09 !cray-libsci_acc" @@ -42,7 +39,6 @@ jobs: 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" - experimental: true env: # Core variables: FC: ftn @@ -103,18 +99,31 @@ jobs: make libs make -C build separate-libs # - # Run examples and tests + # Run examples and tests (expect success) # - - name: Run examples and tests + - name: Run examples and tests (expect success) + id: run-success + if: matrix.config-name != 'cce-gpu-openmp' run: make tests # + # Run examples and tests (expect failure) + # + - name: Run examples and tests (expect failure) + if: steps.run-success.outcome == 'skipped' + run: | + make 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' + 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 check From 43adc58fbd5cbd373f3bc1c031b81861e03a70aa Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Mon, 19 Feb 2024 08:36:39 -0800 Subject: [PATCH 11/20] Add initial unit tests (solvers) (#263) Adds stand-alone unit tests for solvers alone (in tests/rte[sl]w_unit_tests.F90). Longwave solvers are tested for correctness in radiative equilibrium and for invariance to choices (vertical orientation, problem size...). Shortwave solvers are tested only for invariance. --- CITATION.cff | 4 +- Makefile | 8 +- rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 | 4 +- rte-frontend/mo_rte_sw.F90 | 2 +- tests/Makefile | 21 +- tests/all_tests.sh | 11 +- tests/check_equivalence.F90 | 119 ++---- tests/check_variants.F90 | 117 +++--- tests/mo_gas_optics_defs_rrtmgp.F90 | 250 ++++++++++++ tests/mo_testing_utils.F90 | 276 +++++++++++++ tests/rte_lw_solver_unit_tests.F90 | 383 ++++++++++++++++++ tests/rte_optic_prop_unit_tests.F90 | 231 +++++++++++ tests/rte_sw_solver_unit_tests.F90 | 346 ++++++++++++++++ ...test_zenith_angle_spherical_correction.F90 | 4 +- 14 files changed, 1625 insertions(+), 151 deletions(-) create mode 100644 tests/mo_gas_optics_defs_rrtmgp.F90 create mode 100644 tests/mo_testing_utils.F90 create mode 100644 tests/rte_lw_solver_unit_tests.F90 create mode 100644 tests/rte_optic_prop_unit_tests.F90 create mode 100644 tests/rte_sw_solver_unit_tests.F90 diff --git a/CITATION.cff b/CITATION.cff index fddcb6833..b574ec618 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -25,12 +25,14 @@ authors: given-names: Igor N. - family-names: Romero given-names: Nicols A. + - family-names: Kosukhin + given-names: Sergey S. - family-names: Wehe given-names: Andre type: software repository-code: "https://github.com/earth-system-radiaton/rte-rrtmgp" license: BSD-3-Clause -date-released: "2023-05-30" +date-released: "2023-11-27" version: 1.7 abstract: "RTE+RRTMGP is a set of codes for computing radiative fluxes in planetary atmospheres. diff --git a/Makefile b/Makefile index c36bd1179..668b88f93 100644 --- a/Makefile +++ b/Makefile @@ -2,7 +2,7 @@ # Top-level Makefile # .PHONY: libs tests check docs -all: libs tests check docs +all: libs tests check docs libs: make -C build -j @@ -11,21 +11,21 @@ libs: make -C examples/rfmip-clear-sky -j tests: + make -C tests tests make -C examples/rfmip-clear-sky tests make -C examples/all-sky tests - make -C tests tests check: + make -C tests check make -C examples/rfmip-clear-sky check make -C examples/all-sky check - make -C tests check docs: @cd doc; ./build_documentation.sh clean: make -C build clean + make -C tests clean make -C examples/rfmip-clear-sky clean make -C examples/all-sky clean - make -C tests clean rm -rf public diff --git a/rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 b/rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 index 9aeb476f5..c83601586 100644 --- a/rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 +++ b/rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 @@ -1728,10 +1728,10 @@ end function is_loaded ! subroutine finalize(this) class(ty_gas_optics_rrtmgp), intent(inout) :: this - real(wp), dimension(:), allocatable :: press_ref, press_ref_log, temp_ref if(this%is_loaded()) then !$acc exit data delete(this%gas_names, this%vmr_ref, this%flavor) & + !$acc delete(this%press_ref, this%press_ref_log, this%temp_ref) & !$acc delete(this%gpoint_flavor, this%kmajor) & !$acc delete(this%minor_limits_gpt_lower) & !$acc delete(this%minor_scales_with_density_lower, this%scale_by_complement_lower) & @@ -1742,6 +1742,7 @@ subroutine finalize(this) !$acc delete(this%idx_minor_upper, this%idx_minor_scaling_upper) & !$acc delete(this%kminor_start_upper, this%kminor_upper) !$omp target exit data map(release:this%gas_names, this%vmr_ref, this%flavor) & + !$omp map(release:this%press_ref, this%press_ref_log, this%temp_ref) !$omp map(release:this%gpoint_flavor, this%kmajor) & !$omp map(release:this%minor_limits_gpt_lower) & !$omp map(release:this%minor_scales_with_density_lower, this%scale_by_complement_lower) & @@ -1752,6 +1753,7 @@ subroutine finalize(this) !$omp map(release:this%idx_minor_upper, this%idx_minor_scaling_upper) & !$omp map(release:this%kminor_start_upper, this%kminor_upper) deallocate(this%gas_names, this%vmr_ref, this%flavor, this%gpoint_flavor, this%kmajor) + deallocate(this%press_ref, this%press_ref_log, this%temp_ref) deallocate(this%minor_limits_gpt_lower, & this%minor_scales_with_density_lower, this%scale_by_complement_lower, & this%idx_minor_lower, this%idx_minor_scaling_lower, this%kminor_start_lower, this%kminor_lower) diff --git a/rte-frontend/mo_rte_sw.F90 b/rte-frontend/mo_rte_sw.F90 index 961209543..22cdf3aeb 100644 --- a/rte-frontend/mo_rte_sw.F90 +++ b/rte-frontend/mo_rte_sw.F90 @@ -336,7 +336,7 @@ function rte_sw_mu0_full(atmos, top_at_1, & ! type is (ty_fluxes_broadband) if(associated(fluxes%flux_net)) then - !$acc parallel loop collapse(2) copyout(fluxes%flux_net) + !$acc parallel loop collapse(2) copyin(fluxes) copyout(fluxes%flux_net) !$omp target teams distribute parallel do simd collapse(2) do ilev = 1, nlay+1 do icol = 1, ncol diff --git a/tests/Makefile b/tests/Makefile index d433f4725..aaa159f47 100644 --- a/tests/Makefile +++ b/tests/Makefile @@ -32,9 +32,10 @@ VPATH += $(RRTMGP_ROOT)/rrtmgp-frontend:$(RRTMGP_ROOT)/extensions:$(RRTMGP_ROOT) # 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_load_coefficients.o mo_simple_netcdf.o mo_rfmip_io.o -ADDITIONS += mo_testing_io.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 @@ -45,7 +46,7 @@ ADDITIONS += mo_solar_variability.o # # Targets # -all: check_variants check_equivalence test_zenith_angle_spherical_correction +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 @@ -62,14 +63,28 @@ 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.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: 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/" diff --git a/tests/all_tests.sh b/tests/all_tests.sh index 9af6c7021..7752f4064 100644 --- a/tests/all_tests.sh +++ b/tests/all_tests.sh @@ -1,7 +1,10 @@ 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_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 +./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 40270178e..ed6022788 100644 --- a/tests/check_equivalence.F90 +++ b/tests/check_equivalence.F90 @@ -34,13 +34,13 @@ program rte_check_equivalence ty_optical_props_arry, & ty_optical_props_1scl, ty_optical_props_2str, ty_optical_props_nstr use mo_rte_util_array, only: zero_array - use mo_gas_optics_rrtmgp, only: ty_gas_optics_rrtmgp use mo_gas_concentrations, only: ty_gas_concs + use mo_gas_optics_defs, only: gas_optics, load_and_init use mo_source_functions, only: ty_source_func_lw use mo_fluxes, only: ty_fluxes_broadband use mo_rte_lw, only: rte_lw use mo_rte_sw, only: rte_sw - use mo_load_coefficients, only: load_and_init + use mo_testing_utils, only: increment_with_1scl, increment_with_2str, increment_with_nstr use mo_rfmip_io, only: read_size, read_and_block_pt, read_and_block_gases_ty, & read_and_block_lw_bc, read_and_block_sw_bc, determine_gas_names use mo_simple_netcdf, only: get_dim_size, read_field @@ -89,7 +89,6 @@ program rte_check_equivalence ! ! Derived types from the RTE and RRTMGP libraries ! - type(ty_gas_optics_rrtmgp) :: k_dist type(ty_gas_concs) :: gas_concs type(ty_gas_concs), dimension(:), allocatable & :: gas_conc_array @@ -111,7 +110,7 @@ program rte_check_equivalence character(len=32 ), & dimension(:), allocatable :: kdist_gas_names, rfmip_gas_games - character(len=256) :: input_file = "", k_dist_file = "" + character(len=256) :: input_file = "", gas_optics_file = "" ! ---------------------------------------------------------------------------------- ! Code ! ---------------------------------------------------------------------------------- @@ -120,10 +119,10 @@ program rte_check_equivalence ! failed = .false. nUserArgs = command_argument_count() - if (nUserArgs < 2) call stop_on_err("Need to supply input_file k_distribution_file ") + if (nUserArgs < 2) call stop_on_err("Need to supply input_file gas_optics_file ") if (nUserArgs > 3) print *, "Ignoring command line arguments beyond the first three..." call get_command_argument(1,input_file) - call get_command_argument(2,k_dist_file) + call get_command_argument(2,gas_optics_file) if(trim(input_file) == '-h' .or. trim(input_file) == "--help") then call stop_on_err("rte_check_equivalence input_file absorption_coefficients_file") end if @@ -132,7 +131,7 @@ program rte_check_equivalence ! Arrays are allocated as they are read ! call read_size (input_file, ncol, nlay, nexp) - call determine_gas_names(input_file, k_dist_file, 1, kdist_gas_names, rfmip_gas_games) + call determine_gas_names(input_file, gas_optics_file, 1, kdist_gas_names, rfmip_gas_games) call read_and_block_pt (input_file, ncol, p_lay_3d, p_lev_3d, t_lay_3d, t_lev_3d) ! ! Only do the first RFMIP experiment @@ -158,17 +157,17 @@ program rte_check_equivalence deallocate(gas_conc_array) ! ---------------------------------------------------------------------------- ! load data into classes - call load_and_init(k_dist, k_dist_file, gas_concs) - is_sw = k_dist%source_is_external() + call load_and_init(gas_optics, gas_optics_file, gas_concs) + is_sw = gas_optics%source_is_external() is_lw = .not. is_sw - print *, "k-distribution is for the " // merge("longwave ", "shortwave", is_lw) - print *, " pressure limits (Pa):", k_dist%get_press_min(), k_dist%get_press_max() - print *, " temperature limits (K):", k_dist%get_temp_min(), k_dist%get_temp_max() + print *, "gas optics is for the " // merge("longwave ", "shortwave", is_lw) + print *, " pressure limits (Pa):", gas_optics%get_press_min(), gas_optics%get_press_max() + print *, " temperature limits (K):", gas_optics%get_temp_min(), gas_optics%get_temp_max() ! ! Problem sizes ! - nbnd = k_dist%get_nband() - ngpt = k_dist%get_ngpt() + nbnd = gas_optics%get_nband() + ngpt = gas_optics%get_ngpt() top_at_1 = p_lay(1, 1) < p_lay(1, nlay) ! ---------------------------------------------------------------------------- ! @@ -191,7 +190,7 @@ program rte_check_equivalence mu0(:) = cos(abs(sza(:,1)) * acos(-1._wp)/180._wp) else allocate(sfc_t(ncol), sfc_emis(nbnd, ncol)) - call stop_on_err(lw_sources%alloc(ncol, nlay, k_dist)) + call stop_on_err(lw_sources%alloc(ncol, nlay, gas_optics)) call read_and_block_lw_bc(input_file, ncol, bc_3d, sfc_t_3d) ! ! Surface emissivity is spectrally uniform @@ -202,7 +201,7 @@ program rte_check_equivalence end if ! ---------------------------------------------------------------------------- ! - ! Fluxes, heat rates, Jacobians + ! Fluxes, heating rates, Jacobians ! allocate(ref_flux_up(ncol,nlay+1), ref_flux_dn(ncol,nlay+1), & tst_flux_up(ncol,nlay+1), tst_flux_dn(ncol,nlay+1), & @@ -221,9 +220,9 @@ program rte_check_equivalence ! ! initialization, finalization of optical properties ! - call make_optical_props_1scl(k_dist) + call make_optical_props_1scl(gas_optics) call atmos%finalize() - call make_optical_props_1scl(k_dist) + call make_optical_props_1scl(gas_optics) call atmos%set_name("gas only atmosphere") print *, " Intialized atmosphere twice" ! @@ -231,7 +230,7 @@ program rte_check_equivalence ! fluxes%flux_up => ref_flux_up(:,:) fluxes%flux_dn => ref_flux_dn(:,:) - call stop_on_err(k_dist%gas_optics(p_lay, p_lev, & + call stop_on_err(gas_optics%gas_optics(p_lay, p_lev, & t_lay, sfc_t, & gas_concs, & atmos, & @@ -302,7 +301,7 @@ program rte_check_equivalence .not. allclose(tst_flux_dn, ref_flux_dn) ) & call report_err(" halving/doubling fails") - call increment_with_1scl + call increment_with_1scl(atmos) call stop_on_err(rte_lw(atmos, top_at_1, & lw_sources, & sfc_emis, & @@ -311,7 +310,7 @@ program rte_check_equivalence .not. allclose(tst_flux_dn, ref_flux_dn) ) & call report_err(" Incrementing with 1scl fails") - call increment_with_2str + call increment_with_2str(atmos) call stop_on_err(rte_lw(atmos, top_at_1, & lw_sources, & sfc_emis, & @@ -320,7 +319,7 @@ program rte_check_equivalence .not. allclose(tst_flux_dn, ref_flux_dn) ) & call report_err(" Incrementing with 2str fails") - call increment_with_nstr + call increment_with_nstr(atmos) call stop_on_err(rte_lw(atmos, top_at_1, & lw_sources, & sfc_emis, & @@ -344,7 +343,7 @@ program rte_check_equivalence ! ! Increase surface temperature by 1K and recompute fluxes ! - call stop_on_err(k_dist%gas_optics(p_lay, p_lev, & + call stop_on_err(gas_optics%gas_optics(p_lay, p_lev, & t_lay, sfc_t + 1._wp, & gas_concs, & atmos, & @@ -370,9 +369,9 @@ program rte_check_equivalence ! ! initialization, finalization of optical properties ! - call make_optical_props_2str(k_dist) + call make_optical_props_2str(gas_optics) call atmos%finalize() - call make_optical_props_2str(k_dist) + call make_optical_props_2str(gas_optics) print *, " Intialized atmosphere twice" ! @@ -381,7 +380,7 @@ program rte_check_equivalence fluxes%flux_up => ref_flux_up (:,:) fluxes%flux_dn => ref_flux_dn (:,:) fluxes%flux_dn_dir => ref_flux_dir(:,:) - call stop_on_err(k_dist%gas_optics(p_lay, p_lev, & + call stop_on_err(gas_optics%gas_optics(p_lay, p_lev, & t_lay, & gas_concs, & atmos, & @@ -421,7 +420,7 @@ program rte_check_equivalence ! Incrementing ! Threshold of 4x spacing() works in double precision ! - call stop_on_err(k_dist%gas_optics(p_lay, p_lev, & + call stop_on_err(gas_optics%gas_optics(p_lay, p_lev, & t_lay, & gas_concs, & atmos, & @@ -437,7 +436,7 @@ program rte_check_equivalence .not. allclose(tst_flux_dir,ref_flux_dir,tol = 8._wp)) & call report_err(" halving/doubling fails") - call increment_with_1scl + call increment_with_1scl(atmos) call stop_on_err(rte_sw(atmos, top_at_1, & mu0, toa_flux, & sfc_alb_dir, sfc_alb_dif, & @@ -447,23 +446,23 @@ program rte_check_equivalence .not. allclose(tst_flux_dir,ref_flux_dir,tol = 6._wp)) & call report_err(" Incrementing with 1scl fails") - call stop_on_err(k_dist%gas_optics(p_lay, p_lev, & + call stop_on_err(gas_optics%gas_optics(p_lay, p_lev, & t_lay, & gas_concs, & atmos, & toa_flux)) - call increment_with_2str + call increment_with_2str(atmos) if(.not. allclose(tst_flux_up, ref_flux_up, tol = 8._wp) .or. & .not. allclose(tst_flux_dn, ref_flux_dn, tol = 6._wp) .or. & .not. allclose(tst_flux_dir,ref_flux_dir,tol = 6._wp)) & call report_err(" Incrementing with 2str fails") - call stop_on_err(k_dist%gas_optics(p_lay, p_lev, & + call stop_on_err(gas_optics%gas_optics(p_lay, p_lev, & t_lay, & gas_concs, & atmos, & toa_flux)) - call increment_with_nstr + call increment_with_nstr(atmos) call stop_on_err(rte_sw(atmos, top_at_1, & mu0, toa_flux, & sfc_alb_dir, sfc_alb_dif, & @@ -489,7 +488,6 @@ subroutine lw_clear_sky_vr character(32), & dimension(gas_concs%get_num_gases()) & :: gc_gas_names - ! ! Reverse the orientation of the problem ! @@ -510,7 +508,7 @@ subroutine lw_clear_sky_vr call stop_on_err(gas_concs_vr%set_vmr(gc_gas_names(i), vmr)) end do - call stop_on_err(k_dist%gas_optics(p_lay, p_lev, & + call stop_on_err(gas_optics%gas_optics(p_lay, p_lev, & t_lay, sfc_t, & gas_concs_vr, & atmos, & @@ -541,7 +539,7 @@ subroutine lw_clear_sky_subset :: up, dn integer :: i, colS, colE - call stop_on_err(k_dist%gas_optics(p_lay, p_lev, & + call stop_on_err(gas_optics%gas_optics(p_lay, p_lev, & t_lay, sfc_t, & gas_concs, & atmos, & @@ -597,7 +595,7 @@ subroutine sw_clear_sky_vr call stop_on_err(gas_concs_vr%set_vmr(gc_gas_names(i), vmr)) end do - call stop_on_err(k_dist%gas_optics(p_lay, p_lev, & + call stop_on_err(gas_optics%gas_optics(p_lay, p_lev, & t_lay, & gas_concs_vr, & atmos, & @@ -628,9 +626,9 @@ subroutine sw_clear_sky_tsi default_tsi = sum(toa_flux(1, :)) ! Set TSI to half the default - call stop_on_err(k_dist%set_tsi(tsi_scale*default_tsi)) + call stop_on_err(gas_optics%set_tsi(tsi_scale*default_tsi)) ! Redo gas optics - call stop_on_err(k_dist%gas_optics(p_lay, p_lev, & + call stop_on_err(gas_optics%gas_optics(p_lay, p_lev, & t_lay, & gas_concs, & atmos, & @@ -643,7 +641,7 @@ subroutine sw_clear_sky_tsi tst_flux_dn (:,:) = tst_flux_dn (:,:) / tsi_scale tst_flux_dir(:,:) = tst_flux_dir(:,:) / tsi_scale - call stop_on_err(k_dist%set_tsi(default_tsi)) + call stop_on_err(gas_optics%set_tsi(default_tsi)) toa_flux (:,:) = toa_flux(:,:) / tsi_scale end subroutine sw_clear_sky_tsi @@ -673,41 +671,10 @@ subroutine report_err(error_msg) failed = .true. end if end subroutine report_err - ! ---------------------------------------------------------------------------- - subroutine increment_with_1scl - type(ty_optical_props_1scl) :: transparent - - call stop_on_err(transparent%alloc_1scl(ncol, nlay, k_dist)) - call zero_array (ncol, nlay, ngpt, transparent%tau) - call stop_on_err(transparent%increment(atmos)) - call transparent%finalize() - end subroutine increment_with_1scl - ! ------- - subroutine increment_with_2str - type(ty_optical_props_2str) :: transparent - - call stop_on_err(transparent%alloc_2str(ncol, nlay, k_dist)) - call zero_array (ncol, nlay, ngpt, transparent%tau) - call zero_array (ncol, nlay, ngpt, transparent%ssa) - call zero_array (ncol, nlay, ngpt, transparent%g) - call stop_on_err(transparent%increment(atmos)) - call transparent%finalize() - end subroutine increment_with_2str - ! ------- - subroutine increment_with_nstr - type(ty_optical_props_nstr) :: transparent - integer, parameter :: nmom = 4 - call stop_on_err(transparent%alloc_nstr(nmom, ncol, nlay, k_dist)) - call zero_array ( ncol, nlay, ngpt, transparent%tau) - call zero_array ( ncol, nlay, ngpt, transparent%ssa) - call zero_array (nmom, ncol, nlay, ngpt, transparent%p) - call stop_on_err(transparent%increment(atmos)) - call transparent%finalize() - end subroutine increment_with_nstr ! ---------------------------------------------------------------------------- - subroutine make_optical_props_1scl(k_dist) - class (ty_optical_props), intent(in) :: k_dist + subroutine make_optical_props_1scl(gas_optics) + class (ty_optical_props), intent(in) :: gas_optics if(allocated(atmos)) then call atmos%finalize() @@ -719,14 +686,14 @@ subroutine make_optical_props_1scl(k_dist) ! select type(atmos) class is (ty_optical_props_1scl) - call stop_on_err(atmos%alloc_1scl(ncol, nlay, k_dist)) + call stop_on_err(atmos%alloc_1scl(ncol, nlay, gas_optics)) class default call stop_on_err("rte_check_equivalence: Don't recognize the kind of optical properties ") end select end subroutine make_optical_props_1scl ! ---------------------------------------------------------------------------- - subroutine make_optical_props_2str(k_dist) - class (ty_optical_props), intent(in) :: k_dist + subroutine make_optical_props_2str(gas_optics) + class (ty_optical_props), intent(in) :: gas_optics if(allocated(atmos)) then call atmos%finalize() deallocate(atmos) @@ -737,7 +704,7 @@ subroutine make_optical_props_2str(k_dist) ! select type(atmos) class is (ty_optical_props_2str) - call stop_on_err(atmos%alloc_2str(ncol, nlay, k_dist)) + call stop_on_err(atmos%alloc_2str(ncol, nlay, gas_optics)) class default call stop_on_err("rte_check_equivalence: Don't recognize the kind of optical properties ") end select diff --git a/tests/check_variants.F90 b/tests/check_variants.F90 index 72564c2d1..a57b9ec84 100644 --- a/tests/check_variants.F90 +++ b/tests/check_variants.F90 @@ -28,22 +28,21 @@ end subroutine stop_on_err ! Serves also to exercise various code paths ! Longwave: ! omiting level temperatures, use optimal angle, use three-angle integration, -! two-stream solution; reduced-resolution k-distribution +! two-stream solution; reduced-resolution gas optics ! Shortwave: -! reduced-resolution k-distribution +! reduced-resolution gas optics ! program rte_clear_sky_regression use mo_rte_kind, only: wp use mo_optical_props, only: ty_optical_props, & ty_optical_props_arry, & ty_optical_props_1scl, ty_optical_props_2str, ty_optical_props_nstr - use mo_gas_optics_rrtmgp, only: ty_gas_optics_rrtmgp use mo_gas_concentrations, only: ty_gas_concs + use mo_gas_optics_defs, only: gas_optics => gas_optics, load_and_init use mo_source_functions, only: ty_source_func_lw use mo_fluxes, only: ty_fluxes_broadband use mo_rte_lw, only: rte_lw use mo_rte_sw, only: rte_sw - use mo_load_coefficients, only: load_and_init use mo_rfmip_io, only: read_size, read_and_block_pt, read_and_block_gases_ty, & read_and_block_lw_bc, read_and_block_sw_bc, determine_gas_names use mo_simple_netcdf, only: get_dim_size, read_field @@ -91,7 +90,6 @@ program rte_clear_sky_regression ! ! Derived types from the RTE and RRTMGP libraries ! - type(ty_gas_optics_rrtmgp) :: k_dist, k_dist_2 type(ty_gas_concs) :: gas_concs type(ty_gas_concs), dimension(:), allocatable & :: gas_conc_array @@ -112,7 +110,7 @@ program rte_clear_sky_regression character(len=32 ), & dimension(:), allocatable :: kdist_gas_names, rfmip_gas_games - character(len=256) :: input_file = "", k_dist_file = "", k_dist_file_2 = "" + character(len=256) :: input_file = "", gas_optics_file = "", gas_optics_file_2 = "" ! ---------------------------------------------------------------------------------- ! Code ! ---------------------------------------------------------------------------------- @@ -120,10 +118,10 @@ 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 k_distribution_file [k_dist_file_2]") + if (nUserArgs < 2) call stop_on_err("Need to supply input_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,k_dist_file) - if (nUserArgs >= 3) call get_command_argument(3,k_dist_file_2) + 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(trim(input_file) == '-h' .or. trim(input_file) == "--help") then call stop_on_err("clear_sky_regression input_file absorption_coefficients_file") @@ -133,7 +131,7 @@ program rte_clear_sky_regression ! Arrays are allocated as they are read ! call read_size (input_file, ncol, nlay, nexp) - call determine_gas_names(input_file, k_dist_file, 1, kdist_gas_names, rfmip_gas_games) + call determine_gas_names(input_file, gas_optics_file, 1, kdist_gas_names, rfmip_gas_games) call read_and_block_pt (input_file, ncol, p_lay_3d, p_lev_3d, t_lay_3d, t_lev_3d) ! ! Only do the first RFMIP experiment @@ -159,15 +157,15 @@ program rte_clear_sky_regression deallocate(gas_conc_array) ! ---------------------------------------------------------------------------- ! load data into classes - call load_and_init(k_dist, k_dist_file, gas_concs) - is_sw = k_dist%source_is_external() + call load_and_init(gas_optics, gas_optics_file, gas_concs) + is_sw = gas_optics%source_is_external() is_lw = .not. is_sw - print *, "k-distribution is for the " // merge("longwave ", "shortwave", is_lw) + print *, "gas optics is for the " // merge("longwave ", "shortwave", is_lw) ! ! Problem sizes ! - nbnd = k_dist%get_nband() - ngpt = k_dist%get_ngpt() + nbnd = gas_optics%get_nband() + ngpt = gas_optics%get_ngpt() top_at_1 = p_lay(1, 1) < p_lay(1, nlay) ! ---------------------------------------------------------------------------- ! @@ -190,7 +188,7 @@ program rte_clear_sky_regression sun_up) else allocate(sfc_t(ncol), sfc_emis(nbnd, ncol)) - call stop_on_err(lw_sources%alloc(ncol, nlay, k_dist)) + call stop_on_err(lw_sources%alloc(ncol, nlay, gas_optics)) call read_and_block_lw_bc(input_file, ncol, bc_3d, sfc_t_3d) ! ! Surface emissivity is spectrally uniform @@ -211,36 +209,39 @@ program rte_clear_sky_regression fluxes%flux_up => flux_up(:,:) fluxes%flux_dn => flux_dn(:,:) if(is_lw) then - call make_optical_props_1scl(k_dist) + call make_optical_props_1scl(gas_optics) call atmos%set_name("gas only atmosphere") call lw_clear_sky_default call lw_clear_sky_notlev call lw_clear_sky_3ang call lw_clear_sky_optangle call lw_clear_sky_jaco - call make_optical_props_2str(k_dist) + call make_optical_props_2str(gas_optics) call lw_clear_sky_2str - if(len_trim(k_dist_file_2) > 0) then - call load_and_init(k_dist_2, k_dist_file_2, gas_concs) - print *, "Alternate k-distribution is for the " // merge("longwave ", "shortwave", .not. k_dist_2%source_is_external()) - print *, " Resolution :", k_dist_2%get_nband(), k_dist_2%get_ngpt() - ngpt = k_dist_2%get_ngpt() + ! + ! Replaces default gas optics with alternative + ! + if(len_trim(gas_optics_file_2) > 0) then + call load_and_init(gas_optics, gas_optics_file_2, gas_concs) + print *, "Alternate gas optics is for the " // merge("longwave ", "shortwave", gas_optics%source_is_internal()) + print *, " Resolution :", gas_optics%get_nband(), gas_optics%get_ngpt() + ngpt = gas_optics%get_ngpt() call atmos%finalize() - call make_optical_props_1scl(k_dist_2) - call stop_on_err(lw_sources%alloc(ncol, nlay, k_dist_2)) + call make_optical_props_1scl(gas_optics) + call stop_on_err(lw_sources%alloc(ncol, nlay, gas_optics)) call lw_clear_sky_alt end if else - call make_optical_props_2str(k_dist) + call make_optical_props_2str(gas_optics) call sw_clear_sky_default - if(len_trim(k_dist_file_2) > 0) then - call load_and_init(k_dist_2, k_dist_file_2, gas_concs) - print *, "Alternate k-distribution is for the " // merge("longwave ", "shortwave", .not. k_dist_2%source_is_external()) - print *, " Resolution :", k_dist_2%get_nband(), k_dist_2%get_ngpt() + if(len_trim(gas_optics_file_2) > 0) then + call load_and_init(gas_optics, gas_optics_file_2, gas_concs) + print *, "Alternate gas optics is for the " // merge("longwave ", "shortwave", gas_optics%source_is_internal()) + print *, " Resolution :", gas_optics%get_nband(), gas_optics%get_ngpt() call atmos%finalize() - call make_optical_props_2str(k_dist_2) + call make_optical_props_2str(gas_optics) deallocate(toa_flux) - allocate(toa_flux(ncol, k_dist_2%get_ngpt())) + allocate(toa_flux(ncol, gas_optics%get_ngpt())) call sw_clear_sky_alt end if end if @@ -254,7 +255,7 @@ subroutine lw_clear_sky_default real(wp), dimension(ncol, nlay) :: heating_rate fluxes%flux_net => flux_net - call stop_on_err(k_dist%gas_optics(p_lay, p_lev, & + call stop_on_err(gas_optics%gas_optics(p_lay, p_lev, & t_lay, sfc_t, & gas_concs, & atmos, & @@ -291,7 +292,7 @@ end subroutine lw_clear_sky_default ! Clear-sky longwave fluxes, level temperatures provided ! subroutine lw_clear_sky_notlev - call stop_on_err(k_dist%gas_optics(p_lay, p_lev, & + call stop_on_err(gas_optics%gas_optics(p_lay, p_lev, & t_lay, sfc_t, & gas_concs, & atmos, & @@ -308,7 +309,7 @@ end subroutine lw_clear_sky_notlev ! Clear-sky longwave fluxes, all info, three angles ! subroutine lw_clear_sky_3ang - call stop_on_err(k_dist%gas_optics(p_lay, p_lev, & + call stop_on_err(gas_optics%gas_optics(p_lay, p_lev, & t_lay, sfc_t, & gas_concs, & atmos, & @@ -327,13 +328,13 @@ end subroutine lw_clear_sky_3ang ! subroutine lw_clear_sky_optangle real(wp), dimension(ncol, ngpt) :: lw_Ds - call stop_on_err(k_dist%gas_optics(p_lay, p_lev, & + call stop_on_err(gas_optics%gas_optics(p_lay, p_lev, & t_lay, sfc_t, & gas_concs, & atmos, & lw_sources, & tlev = t_lev)) - call stop_on_err(k_dist%compute_optimal_angles(atmos, lw_Ds)) + call stop_on_err(gas_optics%compute_optimal_angles(atmos, lw_Ds)) call stop_on_err(rte_lw(atmos, top_at_1, & lw_sources, & sfc_emis, & @@ -348,7 +349,7 @@ end subroutine lw_clear_sky_optangle subroutine lw_clear_sky_jaco real(wp), dimension(ncol,nlay+1) :: jFluxUp - call stop_on_err(k_dist%gas_optics(p_lay, p_lev, & + call stop_on_err(gas_optics%gas_optics(p_lay, p_lev, & t_lay, sfc_t, & gas_concs, & atmos, & @@ -363,7 +364,7 @@ subroutine lw_clear_sky_jaco call write_broadband_field(input_file, flux_dn, "lw_flux_dn_jaco", "LW flux dn, computing Jaobians") call write_broadband_field(input_file, jFluxUp, "lw_jaco_up" , "Jacobian of LW flux up to surface temperature") - call stop_on_err(k_dist%gas_optics(p_lay, p_lev, & + call stop_on_err(gas_optics%gas_optics(p_lay, p_lev, & t_lay, sfc_t + 1._wp, & gas_concs, & atmos, & @@ -384,7 +385,7 @@ end subroutine lw_clear_sky_jaco ! The second uses the two-stream solver ! subroutine lw_clear_sky_2str - call stop_on_err(k_dist%gas_optics(p_lay, p_lev, & + call stop_on_err(gas_optics%gas_optics(p_lay, p_lev, & t_lay, sfc_t, & gas_concs, & atmos, & @@ -412,12 +413,12 @@ subroutine lw_clear_sky_alt real(wp), dimension(ncol, ngpt) :: lw_Ds fluxes%flux_net => flux_net - call stop_on_err(k_dist_2%gas_optics(p_lay, p_lev, & - t_lay, sfc_t, & - gas_concs, & - atmos, & - lw_sources, & - tlev = t_lev)) + call stop_on_err(gas_optics%gas_optics(p_lay, p_lev, & + t_lay, sfc_t, & + gas_concs, & + atmos, & + lw_sources, & + tlev = t_lev)) call stop_on_err(rte_lw(atmos, top_at_1, & lw_sources, & sfc_emis, & @@ -430,7 +431,7 @@ subroutine lw_clear_sky_alt "lw_flux_hr_alt", "LW heating rate, fewer g-points", & vert_dim_name = "layer") - call stop_on_err(k_dist_2%compute_optimal_angles(atmos, lw_Ds)) + 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, & @@ -442,7 +443,7 @@ subroutine lw_clear_sky_alt call write_broadband_field(input_file, heating_rate, & "lw_flux_hr_alt_oa", "LW heating rate, fewer g-points, opt. angle", & vert_dim_name = "layer") - call k_dist_2%finalize() + call gas_optics%finalize() end subroutine lw_clear_sky_alt ! ---------------------------------------------------------------------------- ! @@ -457,7 +458,7 @@ subroutine sw_clear_sky_default real(wp) :: rrtmgp_tsi type(ty_optical_props_2str) :: atmos2 - call stop_on_err(k_dist%gas_optics(p_lay, p_lev, & + call stop_on_err(gas_optics%gas_optics(p_lay, p_lev, & t_lay, & gas_concs, & atmos, & @@ -485,10 +486,10 @@ subroutine sw_clear_sky_default end subroutine sw_clear_sky_default ! ---------------------------------------------------------------------------- subroutine sw_clear_sky_alt - real(wp), dimension(ncol, k_dist_2%get_ngpt()) & + real(wp), dimension(ncol, gas_optics%get_ngpt()) & :: rfmip_tsi_scale real(wp) :: rrtmgp_tsi - call stop_on_err(k_dist_2%gas_optics(p_lay, p_lev, & + call stop_on_err(gas_optics%gas_optics(p_lay, p_lev, & t_lay, & gas_concs, & atmos, & @@ -497,7 +498,7 @@ subroutine sw_clear_sky_alt ! Scaling factor for column dependence of TSI in RFMIP ! rrtmgp_tsi = sum(toa_flux(1,:)) - rfmip_tsi_scale(:,:) = spread(tsi_3d(:,1)/rrtmgp_tsi, dim=2, ncopies=k_dist_2%get_ngpt()) + 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, & @@ -515,8 +516,8 @@ subroutine sw_clear_sky_alt call write_broadband_field(input_file, flux_dn, "sw_flux_dn_alt", "SW flux dn, fewer g-points") end subroutine sw_clear_sky_alt ! ---------------------------------------------------------------------------- - subroutine make_optical_props_1scl(k_dist) - class (ty_optical_props), intent(in) :: k_dist + subroutine make_optical_props_1scl(gas_optics) + class (ty_optical_props), intent(in) :: gas_optics if(allocated(atmos)) then call atmos%finalize() @@ -529,14 +530,14 @@ subroutine make_optical_props_1scl(k_dist) ! select type(atmos) class is (ty_optical_props_1scl) - call stop_on_err(atmos%alloc_1scl(ncol, nlay, k_dist)) + call stop_on_err(atmos%alloc_1scl(ncol, nlay, gas_optics)) class default call stop_on_err("rte_rrtmgp_atmos: Don't recognize the kind of optical properties ") end select end subroutine make_optical_props_1scl ! ---------------------------------------------------------------------------- - subroutine make_optical_props_2str(k_dist) - class (ty_optical_props), intent(in) :: k_dist + subroutine make_optical_props_2str(gas_optics) + class (ty_optical_props), intent(in) :: gas_optics if(allocated(atmos)) then call atmos%finalize() deallocate(atmos) @@ -548,7 +549,7 @@ subroutine make_optical_props_2str(k_dist) ! select type(atmos) class is (ty_optical_props_2str) - call stop_on_err(atmos%alloc_2str(ncol, nlay, k_dist)) + call stop_on_err(atmos%alloc_2str(ncol, nlay, gas_optics)) class default call stop_on_err("rte_rrtmgp_atmos: Don't recognize the kind of optical properties ") end select diff --git a/tests/mo_gas_optics_defs_rrtmgp.F90 b/tests/mo_gas_optics_defs_rrtmgp.F90 new file mode 100644 index 000000000..b99422f94 --- /dev/null +++ b/tests/mo_gas_optics_defs_rrtmgp.F90 @@ -0,0 +1,250 @@ +! This code is part of Radiative Transfer for Energetics (RTE) +! +! Contacts: Robert Pincus and Eli Mlawer +! email: rrtmgp@aer.com +! +! Copyright 2024- Atmospheric and Environmental Research, +! Regents of the University of Colorado, +! Trustees of Columbia University in the City of New York +! 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 +! ---------------------------------------------------------------------------- +! +! This module is intended to support generic codes using RTE-compatibale gas optics types. +! It defines variable `gas_optics` of class `ty_gas_optics` and subroutine +! `load_and_init()` that loads data into the gas_optics variable from a single netCDF file +! This module might be replaced by others that implement the same variable and procedure +! +! Gas optics classes need to be initialized with data; for RRTMGP data comes from a netCDF file. +! The gas optics classes themselves don't include methods for reading the data so we don't conflict with users' +! local environment. This module provides a straight-forward serial implementation of reading the data +! and calling gas_optics%load(). +! +! +module mo_gas_optics_defs + use mo_rte_kind, only: wp, wl + use mo_gas_concentrations, only: ty_gas_concs + use mo_gas_optics_rrtmgp, only: ty_gas_optics_rrtmgp + use mo_testing_utils, only: stop_on_err + ! -------------------------------------------------- + use mo_simple_netcdf, only: read_field, read_char_vec, read_logical_vec, var_exists, get_dim_size + use netcdf + implicit none + + type(ty_gas_optics_rrtmgp) :: gas_optics + private + public :: gas_optics, load_and_init + +contains + !-------------------------------------------------------------------------------------------------------------------- + ! read optical coefficients from NetCDF file + subroutine load_and_init(kdist, filename, available_gases) + class(ty_gas_optics_rrtmgp), intent(inout) :: kdist + character(len=*), intent(in ) :: filename + class(ty_gas_concs), intent(in ) :: available_gases ! Which gases does the host model have available? + ! -------------------------------------------------- + ! + ! Variables that will be passed to gas_optics%load() + ! + character(len=32), dimension(:), allocatable :: gas_names + integer, dimension(:,:,:), allocatable :: key_species + integer, dimension(:,: ), allocatable :: band2gpt + real(wp), dimension(:,: ), allocatable :: band_lims + real(wp) :: press_ref_trop, temp_ref_p, temp_ref_t + real(wp), dimension(: ), allocatable :: press_ref + real(wp), dimension(: ), allocatable :: temp_ref + real(wp), dimension(:,:,: ), allocatable :: vmr_ref + real(wp), dimension(:,:,:,:), allocatable :: kmajor + + character(len=32), dimension(:), allocatable :: gas_minor, identifier_minor + character(len=32), dimension(:), allocatable :: minor_gases_lower, minor_gases_upper + integer, dimension(:,:), allocatable :: minor_limits_gpt_lower, minor_limits_gpt_upper + logical(wl), dimension(:), allocatable :: minor_scales_with_density_lower, minor_scales_with_density_upper + character(len=32), dimension(:), allocatable :: scaling_gas_lower, scaling_gas_upper + logical(wl), dimension(:), allocatable :: scale_by_complement_lower, scale_by_complement_upper + integer, dimension(:), allocatable :: kminor_start_lower, kminor_start_upper + real(wp), dimension(:,:,:), allocatable :: kminor_lower, kminor_upper + + real(wp), dimension(:,:,: ), allocatable :: rayl_lower, rayl_upper + real(wp), dimension(: ), allocatable :: solar_quiet, solar_facular, solar_sunspot + real(wp) :: tsi_default, mg_default, sb_default + real(wp), dimension(:,: ), allocatable :: totplnk + real(wp), dimension(:,:,:,:), allocatable :: planck_frac + real(wp), dimension(:,:) , allocatable :: optimal_angle_fit + + ! ----------------- + ! + ! Book-keeping variables + ! + integer :: ncid + integer :: ntemps, & + npress, & + nabsorbers, & + nextabsorbers, & + nminorabsorbers, & + nmixingfracs, & + nlayers, & + nbnds, & + ngpts, & + npairs, & + nminor_absorber_intervals_lower, & + nminor_absorber_intervals_upper, & + ncontributors_lower, & + ncontributors_upper, & + ninternalSourcetemps, & + nfit_coeffs + ! -------------------------------------------------- + ! + ! How big are the various arrays? + ! + if(nf90_open(trim(fileName), NF90_NOWRITE, ncid) /= NF90_NOERR) & + call stop_on_err("load_and_init(): can't open file " // trim(fileName)) + ntemps = get_dim_size(ncid,'temperature') + npress = get_dim_size(ncid,'pressure') + nabsorbers = get_dim_size(ncid,'absorber') + nminorabsorbers = get_dim_size(ncid,'minor_absorber') + nextabsorbers = get_dim_size(ncid,'absorber_ext') + nmixingfracs = get_dim_size(ncid,'mixing_fraction') + nlayers = get_dim_size(ncid,'atmos_layer') + nbnds = get_dim_size(ncid,'bnd') + ngpts = get_dim_size(ncid,'gpt') + npairs = get_dim_size(ncid,'pair') + nminor_absorber_intervals_lower & + = get_dim_size(ncid,'minor_absorber_intervals_lower') + nminor_absorber_intervals_upper & + = get_dim_size(ncid,'minor_absorber_intervals_upper') + ninternalSourcetemps & + = get_dim_size(ncid,'temperature_Planck') + ncontributors_lower = get_dim_size(ncid,'contributors_lower') + ncontributors_upper = get_dim_size(ncid,'contributors_upper') + nfit_coeffs = get_dim_size(ncid,'fit_coeffs') ! Will be 0 for SW + + ! ----------------- + ! + ! Read the many arrays + ! + gas_names = read_char_vec(ncid, 'gas_names', nabsorbers) + key_species = int(read_field(ncid, 'key_species', 2, nlayers, nbnds)) + band_lims = read_field(ncid, 'bnd_limits_wavenumber', 2, nbnds) + band2gpt = int(read_field(ncid, 'bnd_limits_gpt', 2, nbnds)) + press_ref = read_field(ncid, 'press_ref', npress) + temp_ref = read_field(ncid, 'temp_ref', ntemps) + temp_ref_p = read_field(ncid, 'absorption_coefficient_ref_P') + temp_ref_t = read_field(ncid, 'absorption_coefficient_ref_T') + press_ref_trop = read_field(ncid, 'press_ref_trop') + kminor_lower = read_field(ncid, 'kminor_lower', & + ncontributors_lower, nmixingfracs, ntemps) + kminor_upper = read_field(ncid, 'kminor_upper', & + ncontributors_upper, nmixingfracs, ntemps) + gas_minor = read_char_vec(ncid, 'gas_minor', nminorabsorbers) + identifier_minor = read_char_vec(ncid, 'identifier_minor', nminorabsorbers) + minor_gases_lower = read_char_vec(ncid, 'minor_gases_lower', nminor_absorber_intervals_lower) + minor_gases_upper = read_char_vec(ncid, 'minor_gases_upper', nminor_absorber_intervals_upper) + minor_limits_gpt_lower & + = int(read_field(ncid, 'minor_limits_gpt_lower', npairs,nminor_absorber_intervals_lower)) + minor_limits_gpt_upper & + = int(read_field(ncid, 'minor_limits_gpt_upper', npairs,nminor_absorber_intervals_upper)) + minor_scales_with_density_lower & + = read_logical_vec(ncid, 'minor_scales_with_density_lower', nminor_absorber_intervals_lower) + minor_scales_with_density_upper & + = read_logical_vec(ncid, 'minor_scales_with_density_upper', nminor_absorber_intervals_upper) + scale_by_complement_lower & + = read_logical_vec(ncid, 'scale_by_complement_lower', nminor_absorber_intervals_lower) + scale_by_complement_upper & + = read_logical_vec(ncid, 'scale_by_complement_upper', nminor_absorber_intervals_upper) + scaling_gas_lower & + = read_char_vec(ncid, 'scaling_gas_lower', nminor_absorber_intervals_lower) + scaling_gas_upper & + = read_char_vec(ncid, 'scaling_gas_upper', nminor_absorber_intervals_upper) + kminor_start_lower & + = int(read_field(ncid, 'kminor_start_lower', nminor_absorber_intervals_lower)) + kminor_start_upper & + = int(read_field(ncid, 'kminor_start_upper', nminor_absorber_intervals_upper)) + vmr_ref = read_field(ncid, 'vmr_ref', nlayers, nextabsorbers, ntemps) + + kmajor = read_field(ncid, 'kmajor', ngpts, nmixingfracs, npress+1, ntemps) + if(var_exists(ncid, 'rayl_lower')) then + rayl_lower = read_field(ncid, 'rayl_lower', ngpts, nmixingfracs, ntemps) + rayl_upper = read_field(ncid, 'rayl_upper', ngpts, nmixingfracs, ntemps) + end if + ! -------------------------------------------------- + ! + ! Initialize the gas optics class with data. The calls look slightly different depending + ! on whether the radiation sources are internal to the atmosphere (longwave) or external (shortwave) + ! gas_optics%load() returns a string; a non-empty string indicates an error. + ! + if(var_exists(ncid, 'totplnk')) then + ! + ! If there's a totplnk variable in the file it's a longwave (internal sources) type + ! + totplnk = read_field(ncid, 'totplnk', ninternalSourcetemps, nbnds) + planck_frac = read_field(ncid, 'plank_fraction', ngpts, nmixingfracs, npress+1, ntemps) + optimal_angle_fit = read_field(ncid, 'optimal_angle_fit', nfit_coeffs, nbnds) + call stop_on_err(kdist%load(available_gases, & + gas_names, & + key_species, & + band2gpt, & + band_lims, & + press_ref, & + press_ref_trop, & + temp_ref, & + temp_ref_p, temp_ref_t, & + vmr_ref, kmajor, & + kminor_lower, kminor_upper, & + gas_minor,identifier_minor, & + minor_gases_lower, minor_gases_upper, & + minor_limits_gpt_lower, & + minor_limits_gpt_upper, & + minor_scales_with_density_lower, & + minor_scales_with_density_upper, & + scaling_gas_lower, scaling_gas_upper, & + scale_by_complement_lower, & + scale_by_complement_upper, & + kminor_start_lower, & + kminor_start_upper, & + totplnk, planck_frac, & + rayl_lower, rayl_upper, & + optimal_angle_fit)) + else + ! + ! Solar source doesn't have an dependencies yet + ! + solar_quiet = read_field(ncid, 'solar_source_quiet', ngpts) + solar_facular = read_field(ncid, 'solar_source_facular', ngpts) + solar_sunspot = read_field(ncid, 'solar_source_sunspot', ngpts) + tsi_default = read_field(ncid, 'tsi_default') + mg_default = read_field(ncid, 'mg_default') + sb_default = read_field(ncid, 'sb_default') + call stop_on_err(kdist%load(available_gases, & + gas_names, & + key_species, & + band2gpt, & + band_lims, & + press_ref, & + press_ref_trop, & + temp_ref, & + temp_ref_p, temp_ref_t, & + vmr_ref, kmajor, & + kminor_lower, kminor_upper, & + gas_minor,identifier_minor,& + minor_gases_lower, minor_gases_upper, & + minor_limits_gpt_lower, & + minor_limits_gpt_upper, & + minor_scales_with_density_lower, & + minor_scales_with_density_upper, & + scaling_gas_lower, scaling_gas_upper, & + scale_by_complement_lower, & + scale_by_complement_upper, & + kminor_start_lower, & + kminor_start_upper, & + solar_quiet, solar_facular, solar_sunspot, & + tsi_default, mg_default, sb_default, & + rayl_lower, rayl_upper)) + end if + ! -------------------------------------------------- + ncid = nf90_close(ncid) + end subroutine load_and_init + !-------------------------------------------------------------------------------------------------------------------- +end module mo_gas_optics_defs \ No newline at end of file diff --git a/tests/mo_testing_utils.F90 b/tests/mo_testing_utils.F90 new file mode 100644 index 000000000..07adffa39 --- /dev/null +++ b/tests/mo_testing_utils.F90 @@ -0,0 +1,276 @@ +! This code is part of Radiative Transfer for Energetics (RTE) +! +! Contacts: Robert Pincus and Eli Mlawer +! email: rrtmgp@aer.com +! +! Copyright 2023- Atmospheric and Environmental Research, +! Regents of the University of Colorado, +! Trustees of Columbia University in the City of New York +! 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_testing_utils + use iso_fortran_env, only: error_unit + use mo_rte_kind, only: wp + use mo_rte_util_array, only: zero_array + use mo_optical_props, only: ty_optical_props_arry, ty_optical_props_1scl, & + ty_optical_props_2str, ty_optical_props_nstr + use mo_source_functions, only: ty_source_func_lw + implicit none + private + public :: allclose, ops_match, check_fluxes + public :: stop_on_err, report_err + public :: increment_with_1scl, increment_with_2str, increment_with_nstr, vr + + interface allclose + module procedure allclose_1, allclose_2, allclose_3, allclose_4 + end interface allclose + + interface ops_match + module procedure ops_match_1scl, ops_match_2str, ops_match_nstr + end interface ops_match + + interface check_fluxes + module procedure check_fluxes_1pair, check_fluxes_2pair + end interface check_fluxes +contains + ! ---------------------------------------------------------------------------- + ! + ! Compare two arrays; return false if abs(x-y) > tol*spacing(x) for any element + ! + ! ---------------------------------------------------------------------------- + logical function allclose_1(tst_array, ref_array, tol) + real(wp), dimension(:), intent(in) :: tst_array, ref_array + real(wp), optional, intent(in) :: tol + + real(wp) :: tolerance + if (present(tol)) then + tolerance = tol + else + tolerance = 2._wp + end if + + allclose_1 = all(abs(tst_array-ref_array) <= tolerance * spacing(ref_array)) + end function allclose_1 + ! ---------------------------------------------------------------------------- + logical function allclose_2(tst_array, ref_array, tol) + real(wp), dimension(:,:), intent(in) :: tst_array, ref_array + real(wp), optional, intent(in) :: tol + + real(wp) :: tolerance + if (present(tol)) then + tolerance = tol + else + tolerance = 2._wp + end if + + allclose_2= all(abs(tst_array-ref_array) <= tolerance * spacing(ref_array)) + end function allclose_2 + ! ---------------------------------------------------------------------------- + logical function allclose_3(tst_array, ref_array, tol) + real(wp), dimension(:,:,:), intent(in) :: tst_array, ref_array + real(wp), optional, intent(in) :: tol + + real(wp) :: tolerance + if (present(tol)) then + tolerance = tol + else + tolerance = 2._wp + end if + + allclose_3= all(abs(tst_array-ref_array) <= tolerance * spacing(ref_array)) + end function allclose_3 + ! ---------------------------------------------------------------------------- + logical function allclose_4(tst_array, ref_array, tol) + real(wp), dimension(:,:,:,:), intent(in) :: tst_array, ref_array + real(wp), optional, intent(in) :: tol + + real(wp) :: tolerance + if (present(tol)) then + tolerance = tol + else + tolerance = 2._wp + end if + + allclose_4= all(abs(tst_array-ref_array) <= tolerance * spacing(ref_array)) + end function allclose_4 + ! ---------------------------------------------------------------------------- + ! + ! Compare two sets of optical properties; return false if abs(x-y) > tol*spacing(x) for any element + ! + ! ---------------------------------------------------------------------------- + logical function ops_match_1scl(tst_values, ref_values, tol) + class(ty_optical_props_1scl), intent(in) :: tst_values, ref_values + real(wp), optional, intent(in) :: tol + + ops_match_1scl = allclose(tst_values%tau, ref_values%tau, tol) + end function ops_match_1scl + ! ---------------------------------------------------------------------------- + logical function ops_match_2str(tst_values, ref_values, tol) + class(ty_optical_props_2str), intent(in) :: tst_values, ref_values + real(wp), optional, intent(in) :: tol + + ops_match_2str = allclose(tst_values%tau, ref_values%tau, tol) .and. & + allclose(tst_values%ssa, ref_values%ssa, tol) .and. & + allclose(tst_values%g , ref_values%g , tol) + end function ops_match_2str + ! ---------------------------------------------------------------------------- + logical function ops_match_nstr(tst_values, ref_values, tol) + class(ty_optical_props_nstr), intent(in) :: tst_values, ref_values + real(wp), optional, intent(in) :: tol + + ops_match_nstr = allclose(tst_values%tau, ref_values%tau, tol) .and. & + allclose(tst_values%ssa, ref_values%ssa, tol) .and. & + allclose(tst_values%p , ref_values%p , tol) + end function ops_match_nstr + ! ---------------------------------------------------------------------------- + + ! ---------------------------------------------------------------------------- + ! + ! Error report - print to screen with or without exit + ! + ! ---------------------------------------------------------------------------- + subroutine report_err(error_msg) + use iso_fortran_env, only : error_unit + character(len=*), intent(in) :: error_msg + + if(error_msg /= "") then + write (error_unit,*) trim(error_msg) + end if + end subroutine report_err + ! ---------------------------------------------------------------------------- + subroutine stop_on_err(error_msg) + use iso_fortran_env, only : error_unit + character(len=*), intent(in) :: error_msg + + if(error_msg /= "") then + write (error_unit,*) trim(error_msg) + write (error_unit,*) "unit tests stopping" + error stop 1 + end if + end subroutine stop_on_err + ! ---------------------------------------------------------------------------- + subroutine check_fluxes_1pair(flux_1, flux_2, status, message) + real(wp), dimension(:,:), intent(in) :: flux_1, flux_2 + logical :: status + character(len=*), intent(in) :: message + + if(.not. allclose(flux_1, flux_2)) then + status = .false. + print *, "check_fluxes: max diffs rel. to scaling: ", & + maxval(abs(flux_1 - flux_2)/spacing(flux_1)) + call report_err(" " // trim(message)) + end if + end subroutine check_fluxes_1pair + ! ---------------------------------------------------------------------------- + subroutine check_fluxes_2pair(flux_1, flux_2, flux_3, flux_4, status, message) + real(wp), dimension(:,:), intent(in) :: flux_1, flux_2, flux_3, flux_4 + logical :: status + character(len=*), intent(in) :: message + + if(.not. (allclose(flux_1, flux_2) .and. & + allclose(flux_3, flux_4))) then + status = .false. + print *, "check_fluxes: max diffs rel. to scaling: ", & + maxval(abs(flux_1 - flux_2)/spacing(flux_1)), & + maxval(abs(flux_3 - flux_4)/spacing(flux_3)) + call report_err(" " // trim(message)) + end if + end subroutine check_fluxes_2pair + ! ---------------------------------------------------------------------------- + ! + ! Adding transparent (tau = 0) optical properties + ! These routines test allocation, validation, incrementing, and + ! finalization for optical properties + ! Fluxes should not change + ! Should these be extended to test end-to-end with GPUs? + ! + ! ---------------------------------------------------------------------------- + subroutine increment_with_1scl(atmos) + class(ty_optical_props_arry), intent(inout) :: atmos + + ! Local variable + type(ty_optical_props_1scl) :: transparent + integer :: ncol, nlay, ngpt + ncol = atmos%get_ncol() + nlay = atmos%get_nlay() + ngpt = atmos%get_ngpt() + + call stop_on_err(transparent%alloc_1scl(ncol, nlay, atmos)) + call zero_array (ncol, nlay, ngpt, transparent%tau) + call stop_on_err(transparent%increment(atmos)) + call stop_on_err(atmos%validate()) + call transparent%finalize() + end subroutine increment_with_1scl + ! ------- + subroutine increment_with_2str(atmos) + class(ty_optical_props_arry), intent(inout) :: atmos + + ! Local variable + type(ty_optical_props_2str) :: transparent + integer :: ncol, nlay, ngpt + ncol = atmos%get_ncol() + nlay = atmos%get_nlay() + ngpt = atmos%get_ngpt() + + call stop_on_err(transparent%alloc_2str(ncol, nlay, atmos)) + call zero_array (ncol, nlay, ngpt, transparent%tau) + call zero_array (ncol, nlay, ngpt, transparent%ssa) + call zero_array (ncol, nlay, ngpt, transparent%g) + call stop_on_err(transparent%increment(atmos)) + call stop_on_err(atmos%validate()) + call transparent%finalize() + end subroutine increment_with_2str + ! ------- + subroutine increment_with_nstr(atmos) + class(ty_optical_props_arry), intent(inout) :: atmos + + ! Local variable + type(ty_optical_props_nstr) :: transparent + integer, parameter :: nmom = 4 + integer :: ncol, nlay, ngpt + ncol = atmos%get_ncol() + nlay = atmos%get_nlay() + ngpt = atmos%get_ngpt() + + call stop_on_err(transparent%alloc_nstr(nmom, ncol, nlay, atmos)) + call zero_array ( ncol, nlay, ngpt, transparent%tau) + call zero_array ( ncol, nlay, ngpt, transparent%ssa) + call zero_array (nmom, ncol, nlay, ngpt, transparent%p) + call stop_on_err(transparent%increment(atmos)) + call stop_on_err(atmos%validate()) + call transparent%finalize() + end subroutine increment_with_nstr + ! ---------------------------------------------------------------------------- + ! + ! Vertically reverse optical properties + ! + subroutine vr(atmos, sources) + class(ty_optical_props_arry), intent(inout) :: atmos + type(ty_source_func_lw), optional, & + intent(inout) :: sources + + integer :: nlay + ! ----------------------- + nlay = atmos%get_nlay() + + atmos%tau(:,:,:) = atmos%tau(:,nlay:1:-1,:) + + select type (atmos) + type is (ty_optical_props_2str) + atmos%ssa(:,:,:) = atmos%ssa(:,nlay:1:-1,:) + atmos%g (:,:,:)= atmos%g (:,nlay:1:-1,:) + type is (ty_optical_props_nstr) + atmos%ssa(:,:,:) = atmos%ssa(:,nlay:1:-1,:) + atmos%p(:,:,:,:) = atmos%p(:,:,nlay:1:-1,:) + end select + + if(present(sources)) then + sources%lev_source(:,:,:) = sources%lev_source(:,nlay+1:1:-1,:) + sources%lay_source(:,:,:) = sources%lay_source(:,nlay :1:-1,:) + end if + end subroutine vr + ! ---------------------------------------------------------------------------- +end module mo_testing_utils diff --git a/tests/rte_lw_solver_unit_tests.F90 b/tests/rte_lw_solver_unit_tests.F90 new file mode 100644 index 000000000..fc18afc44 --- /dev/null +++ b/tests/rte_lw_solver_unit_tests.F90 @@ -0,0 +1,383 @@ +! This code is part of Radiative Transfer for Energetics (RTE) +! +! Contacts: Robert Pincus and Eli Mlawer +! email: rrtmgp@aer.com +! +! Copyright 2023- Atmospheric and Environmental Research, +! Regents of the University of Colorado, +! Trustees of Columbia University in the City of New York +! 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 +! ---------------------------------------------------------------------------- +program rte_lw_solver_unit_tests + ! + ! Exercise various paths through RTE LW solvers + ! Tests are run on an idealized problem (radiative equilibirum) + ! and checked for correctness against known analytic solution + ! Beyond correctness tests check invariance, e.g. with respect to vertical ordering + ! + use mo_rte_kind, only: wp + use mo_optical_props, only: ty_optical_props, & + ty_optical_props_arry, & + ty_optical_props_1scl, ty_optical_props_2str, ty_optical_props_nstr + use mo_rte_util_array, only: zero_array + use mo_source_functions, only: ty_source_func_lw + use mo_fluxes, only: ty_fluxes_broadband + use mo_rte_lw, only: rte_lw + use mo_testing_utils, only: allclose, stop_on_err, report_err, check_fluxes, & + vr, & + increment_with_1scl, increment_with_2str, increment_with_nstr + implicit none + ! ---------------------------------------------------------------------------------- + ! + ! Longwave tests use gray radiative equilibrium from + ! e.g. Weaver and Rmanathan 1995 https://doi.org/10.1029/95JD00770 + ! Net flux is constant with height, OLR is known from surface temperature + ! Tests include + ! Solutions match analytic results + ! Net fluxes = down-up when computed in various combos (net only, up/down only, all three) + ! using ty_broadband_flues + ! Answers are invariant to + ! Extracting subsets + ! Vertical orientation + ! Adding transparent optical properties + ! Longwave specific tests: + ! Computing the Jacibian doesn't change fluxes + ! Fluxes inferred from Jacobian are close to fluxes with perturbed surface T. (TODO) + ! + ! Other possibilites: + ! Vertical discretization? Maybe just check boundary fluxes + ! Test the application of the boundary condition? + + real(wp), parameter :: pi = acos(-1._wp) + integer, parameter :: ncol = 8, nlay = 16 + integer :: icol, ilay + ! + ! Longwave tests - gray radiative equilibrium + ! + real(wp), parameter :: sigma = 5.670374419e-8_wp, & ! Stefan-Boltzmann constant + D = 1.66_wp ! Diffusivity angle, from single-angle RRTMGP solver + real(wp), dimension( ncol), parameter :: sfc_t = [(285._wp, icol = 1,ncol/2), & + (310._wp, icol = 1,ncol/2)] + real(wp), dimension( ncol), parameter :: total_tau = [0.1_wp, 1._wp, 10._wp, 50._wp, & + 0.1_wp, 1._wp, 10._wp, 50._wp] ! Would be nice to parameterize + real(wp), dimension(1,ncol), parameter :: sfc_emis = 1._wp + real(wp), dimension(ncol,1), parameter :: lw_Ds = D ! Diffusivity angle - use default value for all columns + + type(ty_optical_props_1scl) :: lw_atmos + type(ty_source_func_lw) :: lw_sources + type(ty_optical_props_2str) :: sw_atmos + type(ty_fluxes_broadband) :: fluxes + logical :: top_at_1 + real(wp), dimension(ncol,nlay+1), target :: & + ref_flux_up, ref_flux_dn, ref_flux_net, & + tst_flux_up, tst_flux_dn, tst_flux_net, & + jFluxUp + + logical :: passed + + ! ------------------------------------------------------------------------------------------------------ + top_at_1 = .true. + ! ------------------------------------------------------------------------------------------------------ + ! + ! Longwave tests + ! + ! ------------------------------------------------------------------------------------------------------ + print *, "RTE LW solver unit tests" + ! + ! Gray radiative equillibrium + ! + print *, "Using gray radiative equilibrium" + call gray_rad_equil(sfc_t(1:ncol), total_tau(1:ncol), nlay, top_at_1, lw_atmos, lw_sources) + + 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, & + lw_sources, & + sfc_emis, & + fluxes)) + ! + ! Is the solution correct (does it satisfy the profile for radiative equilibrium?) + ! Error reporting happens inside check_gray_rad_equil() + ! + passed = check_gray_rad_equil(sfc_t, total_tau, top_at_1, ref_flux_up, ref_flux_net) + ! ------------------------------------------------------------------------------------ + ! + ! Net fluxes on- vs off-line + ! Are the net fluxes correct? + ! + print *, " Net flux variants" + call check_fluxes(ref_flux_net, ref_flux_dn-ref_flux_up, passed, "net fluxes don't match down-up") + ! + ! Compute only net fluxes + ! + nullify(fluxes%flux_up) + nullify(fluxes%flux_dn) + call stop_on_err(rte_lw(lw_atmos, top_at_1, & + lw_sources, sfc_emis,& + fluxes)) + call check_fluxes(ref_flux_net, ref_flux_dn-ref_flux_up, & + passed, "Net fluxes computed alone doesn'tt match down-up computed separately") + ! + ! Compute only up and down fluxes + ! + fluxes%flux_up => tst_flux_up (:,:) + fluxes%flux_dn => tst_flux_dn (:,:) + call stop_on_err(rte_lw(lw_atmos, top_at_1, & + lw_sources, sfc_emis, & + fluxes)) + call check_fluxes(ref_flux_net, tst_flux_dn-tst_flux_up, & + passed, "LW net fluxes don't match down-up computed together") + ! ------------------------------------------------------- + ! + ! Subsets of atmospheric columns + ! + print *, " Subsetting invariance" + call gray_rad_equil(sfc_t, total_tau, nlay, top_at_1, lw_atmos, lw_sources) + call clear_sky_subset(lw_atmos, lw_sources, sfc_emis, tst_flux_up, tst_flux_dn) + call check_fluxes(tst_flux_up, ref_flux_up, & + tst_flux_dn, ref_flux_dn, & + passed, "Doing problem in subsets fails") + + ! ------------------------------------------------------- + ! + ! Vertically-reverse + ! + 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, & + lw_sources, sfc_emis, & + fluxes)) + ! + ! Seems like these fluxes should be bitwise identical regardless of orientation + ! but nvfortran 22.5 on Levante has differences of up to 3*spacing() + ! + if (.not. allclose(tst_flux_up(:,nlay+1:1:-1), ref_flux_up) .and. & + allclose(tst_flux_dn(:,nlay+1:1:-1), ref_flux_dn, tol=3._wp)) then + passed = .false. + call report_err(" " // "Doing problem upside down fails") + end if + call vr(lw_atmos, lw_sources) + ! ------------------------------------------------------- + ! + ! Computing Jacobian shouldn't change net fluxes + ! + 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, & + lw_sources, & + sfc_emis, & + fluxes, & + flux_up_Jac = jFluxUp)) + call check_fluxes(tst_flux_up, ref_flux_up, tst_flux_dn, ref_flux_dn, & + passed, "Computing Jacobian changes fluxes fails") + ! + ! Increase surface temperature in source function by 1K and recompute fluxes + ! + 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, & + lw_sources, & + sfc_emis, & + fluxes)) + ! + ! Comparision of fluxes with increased surface T aren't expected to match + ! fluxes + their Jacobian w.r.t. surface T exactly + ! + print '(" Jacobian accurate to within ", f6.2, "%")', & + maxval((tst_flux_up - ref_flux_up + jFluxUp)/tst_flux_up * 100._wp) + ! ------------------------------------------------------------------------------------ + ! + ! Using Tang approach for purely absorbing problem should be the same + ! + print *, " Two-stream optical properties" + call gray_rad_equil(sfc_t, total_tau, nlay, top_at_1, lw_atmos, lw_sources) + call stop_on_err(sw_atmos%alloc_2str(ncol, nlay, lw_atmos)) + sw_atmos%tau = lw_atmos%tau + sw_atmos%ssa = 0._wp + sw_atmos%g = 0._wp + + call stop_on_err(rte_lw(sw_atmos, top_at_1, & + lw_sources, & + sfc_emis, & + fluxes, & + flux_up_Jac = jFluxUp)) + call check_fluxes(tst_flux_up, ref_flux_up, tst_flux_dn, ref_flux_dn, & + passed, "Using two-stream properties fails") + call sw_atmos%finalize() + ! ------------------------------------------------------------------------------------ + ! + ! Specifying diffusivity angle + ! + print *, " Specified transport angle" + call stop_on_err(rte_lw(lw_atmos, top_at_1, & + lw_sources, & + sfc_emis, & + fluxes, & + lw_Ds = lw_Ds)) + call check_fluxes(tst_flux_up, ref_flux_up, tst_flux_dn, ref_flux_dn, & + passed, "Specifying diffusivity angle D fails") + + ! ------------------------------------------------------------------------------------ + ! Done + ! + print *, "RTE LW solver unit tests done" + print * + if(.not. passed) error stop 1 + ! ------------------------------------------------------------------------------------ +contains + ! ------------------------------------------------------------------------------------ + ! + ! Define an atmosphere in gray radiative equillibrium + ! See, for example, section 2 of Weaver and Rmanathan 1995 https://doi.org/10.1029/95JD00770 + ! + subroutine gray_rad_equil(sfc_t, total_tau, nlay, top_at_1, atmos, sources) + real(wp), dimension(:), intent(in) :: sfc_t, total_tau + integer, intent(in) :: nlay + logical, intent(in) :: top_at_1 + type(ty_optical_props_1scl), intent(inout) :: atmos + type(ty_source_func_lw), intent(inout) :: sources + + integer :: ncol + real(wp), dimension(size(sfc_t)) :: olr + + ncol = size(sfc_t) + ! + ! Set up a gray spectral distribution - one band, one g-point + ! + call stop_on_err(atmos%init(band_lims_wvn = reshape([0._wp, 3250._wp], shape = [2, 1]), & + band_lims_gpt = reshape([1, 1], shape = [2, 1]), & + name = "Gray atmosphere")) + call stop_on_err(atmos%alloc_1scl(ncol, nlay)) + + ! + ! 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) + + ! + ! Longwave sources - for broadband these are sigma/pi T^4 + ! (isotropic radiation) + ! + olr(:) = gray_rad_equil_olr(sfc_t, total_tau) + + call stop_on_err(sources%alloc(ncol, nlay, atmos)) + sources%sfc_source (:,1) = sigma/pi * sfc_t**4 + sources%sfc_source_Jac(:,1) = 4._wp * sigma/pi * sfc_t**3 + ! + ! 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 + if (.not. top_at_1) then + ! + ! Reverse vertical ordering of source functions + ! + sources%lev_source(:,1:nlay+1,1) = sources%lev_source(:,nlay+1:1:-1,1) + sources%lay_source(:,1:nlay, 1) = sources%lay_source(:,nlay :1:-1,1) + end if + end subroutine gray_rad_equil + ! ------------------------------------------------------------------------------------ + ! + ! Check that solutions are in gray radiative equilibrium + ! We could use this to check heating rates but we'd have to make up pressure levels... + ! + function check_gray_rad_equil(sfc_T, lw_tau, top_at_1, up_flux, net_flux) + real(wp), dimension(:), intent(in) :: sfc_T, lw_tau + real(wp), dimension(:,:), intent(in) :: up_flux, net_flux + logical, intent(in) :: top_at_1 + logical :: check_gray_rad_equil + + logical :: passed + integer :: toa + ! ------------------------------ + check_gray_rad_equil = .true. + toa = merge(1, size(up_flux, 2), top_at_1) + + ! + ! Check top-of-atmosphere energy balance + ! + if(.not. allclose(up_flux(:,toa), & + gray_rad_equil_olr(sfc_t, lw_tau), tol=4._wp)) then + call report_err("OLR is not consistent with gray radiative equilibrium") + check_gray_rad_equil = .false. + end if + ! + ! Check that net fluxes are constant with height + ! Fairly relaxed threshold w.r.t. spacing() because net flux is small relative to + ! large up and down fluxes that vary with tau + ! + if(.not. allclose(net_flux(:,:), & + spread(net_flux(:,1), dim=2, ncopies=size(net_flux,2)), & + tol = 70._wp)) then + call report_err("Net flux not constant with tau in gray radiative equilibrium") + check_gray_rad_equil = .false. + end if + end function check_gray_rad_equil + ! ------------------------------------------------------------------------------------ + ! + ! Incoming energy = OLR in gray radiative equilibirum + ! Equation 6b of Weaver and Rmanathan 1995 https://doi.org/10.1029/95JD00770 with with f0 = OLR + ! + function gray_rad_equil_olr(T, tau) + real(wp), dimension(:), intent(in) :: T, tau + real(wp), dimension(size(T)) :: gray_rad_equil_olr + + gray_rad_equil_olr(:) = (2._wp * sigma * T(:)**4)/(2 + D * tau(:)) + end function gray_rad_equil_olr + ! ------------------------------------------------------------------------------------ + ! + ! Invariance tests + ! + ! ------------------------------------------------------------------------------------ + ! + ! Clear-sky longwave fluxes, half the columns at a time + ! We're counting on ncol being even + ! + subroutine clear_sky_subset(atmos, sources, sfc_emis, flux_up, flux_dn) + type(ty_optical_props_1scl), intent(inout) :: atmos + type(ty_source_func_lw), intent(inout) :: sources + real(wp), dimension(:,:), intent(in ) :: sfc_emis + real(wp), dimension(:,:), intent( out) :: flux_up, flux_dn + + type(ty_optical_props_1scl) :: atmos_subset + type(ty_source_func_lw) :: sources_subset + type(ty_fluxes_broadband) :: fluxes ! Use local variable + real(wp), dimension(atmos%get_ncol()/2, & + atmos%get_nlay()+1), target & + :: up, dn + integer :: i, colS, colE + integer :: ncol + ! ------------------------------ + ncol = atmos%get_ncol() + call stop_on_err(atmos_subset%init(atmos)) + fluxes%flux_up => up + fluxes%flux_dn => dn + + do i = 1, 2 + colS = ((i-1) * ncol/2) + 1 + colE = i * ncol/2 + call stop_on_err(atmos%get_subset (colS, ncol/2, atmos_subset)) + call stop_on_err(sources%get_subset(colS, ncol/2, sources_subset)) + call stop_on_err(rte_lw(atmos_subset, top_at_1, & + sources_subset, & + sfc_emis(:,colS:colE), & + fluxes)) + flux_up(colS:colE,:) = up + flux_dn(colS:colE,:) = dn + end do + end subroutine clear_sky_subset +end program rte_lw_solver_unit_tests \ No newline at end of file diff --git a/tests/rte_optic_prop_unit_tests.F90 b/tests/rte_optic_prop_unit_tests.F90 new file mode 100644 index 000000000..4fea31ff5 --- /dev/null +++ b/tests/rte_optic_prop_unit_tests.F90 @@ -0,0 +1,231 @@ +! This code is part of Radiative Transfer for Energetics (RTE) +! +! Contacts: Robert Pincus and Eli Mlawer +! email: rrtmgp@aer.com +! +! Copyright 2023- Atmospheric and Environmental Research, +! Regents of the University of Colorado, +! Trustees of Columbia University in the City of New York +! 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 +! ---------------------------------------------------------------------------- +program optical_prop_unit_tests + ! + ! Unit tests for RTE optical properties + ! Incrementing with tranparent medium (tau=0) doesn't change optical props + ! + use mo_rte_kind, only: wp + use mo_optical_props, only: ty_optical_props_arry, & + ty_optical_props_1scl, ty_optical_props_2str, ty_optical_props_nstr + use mo_rte_util_array, only: zero_array + use mo_testing_utils, only: allclose, ops_match, stop_on_err, report_err, & + increment_with_1scl, increment_with_2str, increment_with_nstr + + type(ty_optical_props_1scl) :: ref_1scl, tst_1scl + type(ty_optical_props_2str) :: ref_2str, tst_2str + type(ty_optical_props_nstr) :: ref_nstr, tst_nstr + integer, parameter :: ncol = 4, nlay = 8, nmom = 4 + integer :: icol, ilay, imom + logical :: passed + real(wp), dimension( ncol), parameter :: total_tau = [0.1_wp, 1._wp, 10._wp, 50._wp] + real(wp), parameter :: g = 0.85_wp, ssa = 1._wp - 1.e-4_wp + ! ---------------------------------------------------------------------------- + print *, "Optical properties unit testing" + ! + ! Set up a gray spectral distribution - one band, one g-point + ! + call stop_on_err(ref_1scl%init(band_lims_wvn = reshape([0._wp, 3250._wp], shape = [2, 1]), & + band_lims_gpt = reshape([1, 1], shape = [2, 1]), & + name = "Gray atmosphere")) + call stop_on_err(ref_1scl%alloc_1scl(ncol, nlay)) + print '(" Problem size: (ncol, nlay, nband, ngpt): ", 4(i2, 2x))', & + ref_1scl%get_ncol(), ref_1scl%get_nlay(), ref_1scl%get_nband(), ref_1scl%get_ngpt() + ! + ! Divide optical depth evenly among layers + ! + ref_1scl%tau(1:ncol,1:nlay,1) = spread(total_tau(1:ncol)/real(nlay, wp), dim=2, ncopies=nlay) + ! + ! 2- and n-stream optical properties + ! + call stop_on_err(ref_2str%alloc_2str(ncol, nlay, ref_1scl)) + ref_2str%tau = ref_1scl%tau + ref_2str%ssa = ssa + ref_2str%g = g + + call stop_on_err(ref_nstr%alloc_nstr(nmom, ncol, nlay, ref_1scl)) + ref_nstr%tau = ref_1scl%tau + ref_nstr%ssa = ssa + ! Henyey-Greenstein phase function + do imom = 1, nmom + ref_nstr%p(imom,:,:,:) = g**imom + end do + + passed = .true. + ! ---------------------------------------------------------------------------- + ! + ! Incrementing with transparent (tau=0) sets of optical properties + ! + ! ---------------------------------------------------------------------------- + print *, " Incrementing 1scl" + ! + ! Increment 1scl + ! + call make_copy_1scl + call increment_with_1scl(tst_1scl) + if(.not. ops_match(tst_1scl, ref_1scl)) then + call report_err("1scl+1scl fails") + passed = .false. + end if + + call make_copy_1scl + call increment_with_2str(tst_1scl) + if(.not. ops_match(tst_1scl, ref_1scl)) then + call report_err("1scl+2str fails") + passed = .false. + end if + + call make_copy_1scl + call increment_with_nstr(tst_1scl) + if(.not. ops_match(tst_1scl, ref_1scl)) then + call report_err("1scl+nstr fails") + passed = .false. + end if + + call tst_1scl%finalize() + ! ---------------------------------------------------------------------------- + print *, " Incrementing 2str" + ! + ! Increment 2str + ! + call make_copy_2str + call increment_with_1scl(tst_2str) + if(.not. ops_match(tst_2str, ref_2str)) then + call report_err("2str+1scl fails") + passed = .false. + end if + + call make_copy_2str + call increment_with_2str(tst_2str) + if(.not. ops_match(tst_2str, ref_2str)) then + call report_err("2str+2str fails") + passed = .false. + end if + + call make_copy_2str + call increment_with_nstr(tst_2str) + if(.not. ops_match(tst_2str, ref_2str)) then + call report_err("2str+nstr fails") + passed = .false. + end if + + call tst_2str%finalize() + ! ---------------------------------------------------------------------------- + print *, " Incrementing nstr" + ! + ! Increment nstr + ! + call make_copy_nstr + call increment_with_1scl(tst_nstr) + if(.not. ops_match(tst_nstr, ref_nstr)) then + call report_err("nstr+1scl fails") + passed = .false. + end if + + call make_copy_nstr + call increment_with_2str(tst_nstr) + if(.not. ops_match(tst_nstr, ref_nstr)) then + call report_err("nstr+2str fails") + passed = .false. + end if + + call make_copy_nstr + call increment_with_nstr(tst_nstr) + if(.not. ops_match(tst_nstr, ref_nstr)) then + call report_err("nstr+nstr fails") + passed = .false. + end if + + call tst_2str%finalize() + ! ---------------------------------------------------------------------------- + print *, " Halving/doubling optical thickness" + ! + ! Adding two media of half optical thickness to recover original values + ! + call make_copy_1scl + tst_1scl%tau = 0.5_wp * tst_1scl%tau + call stop_on_err(tst_1scl%increment(tst_1scl)) + if(.not. ops_match(tst_1scl, ref_1scl)) then + call report_err("1scl half/double fails") + passed = .false. + end if + + call make_copy_2str + tst_2str%tau = 0.5_wp * tst_2str%tau + call stop_on_err(tst_2str%increment(tst_2str)) + if(.not. ops_match(tst_2str, ref_2str)) then + call report_err("2str half/double fails") + passed = .false. + end if + + call make_copy_nstr + tst_nstr%tau = 0.5_wp * tst_nstr%tau + call stop_on_err(tst_nstr%increment(tst_nstr)) + if(.not. ops_match(tst_nstr, ref_nstr)) then + call report_err("nstr half/double fails") + passed = .false. + end if + ! ---------------------------------------------------------------------------- + print *, " Delta scaling" + ! + ! Delta-scale with forward-fraction f=0 (i.e. Rayleigh scattering) + ! + call make_copy_2str + call stop_on_err(tst_2str%delta_scale(spread(spread(spread(0._wp, 1, ncol), 2, nlay), 3, 1))) + if(.not. ops_match(tst_2str, ref_2str)) then + call report_err("2str delta-scaling with f=0 fails") + passed = .false. + end if + ! ---------------------------------------------------------------------------- + if (.not. passed) call stop_on_err("Optical props unit tests fail") + print *, "Optical properties unit testing finished" + print * + ! ---------------------------------------------------------------------------- +contains + ! ---------------------------------------------------------------------------- + ! + ! Make copies of the existing optical depth + ! + subroutine make_copy_1scl + call tst_1scl%finalize() + call stop_on_err(tst_1scl%alloc_1scl(ref_1scl%get_ncol(), ref_1scl%get_nlay(), & + ref_1scl)) + tst_1scl%tau = ref_1scl%tau + end subroutine make_copy_1scl + ! ---------------------------------------------------------------------------- + ! + ! Make copies of the existing optical depth + ! + subroutine make_copy_2str + call tst_2str%finalize() + call stop_on_err(tst_2str%alloc_2str(ref_2str%get_ncol(), ref_2str%get_nlay(), & + ref_2str)) + tst_2str%tau = ref_2str%tau + tst_2str%ssa = ref_2str%ssa + tst_2str%g = ref_2str%g + end subroutine make_copy_2str + ! ---------------------------------------------------------------------------- + ! + ! Make copies of the existing optical depth + ! + subroutine make_copy_nstr + call tst_nstr%finalize() + call stop_on_err(tst_nstr%alloc_nstr(ref_nstr%get_nmom(), ref_nstr%get_ncol(), ref_nstr%get_nlay(), & + ref_2str)) + tst_nstr%tau = ref_nstr%tau + tst_nstr%ssa = ref_nstr%ssa + tst_nstr%p = ref_nstr%p + end subroutine make_copy_nstr + ! ---------------------------------------------------------------------------- +end program optical_prop_unit_tests \ No newline at end of file diff --git a/tests/rte_sw_solver_unit_tests.F90 b/tests/rte_sw_solver_unit_tests.F90 new file mode 100644 index 000000000..196fd1f4a --- /dev/null +++ b/tests/rte_sw_solver_unit_tests.F90 @@ -0,0 +1,346 @@ +! This code is part of Radiative Transfer for Energetics (RTE) +! +! Contacts: Robert Pincus and Eli Mlawer +! email: rrtmgp@aer.com +! +! Copyright 2023- Atmospheric and Environmental Research, +! Regents of the University of Colorado, +! Trustees of Columbia University in the City of New York +! 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 +! ---------------------------------------------------------------------------- +program rte_sw_solver_unit_tests + ! + ! Exercise various paths through RTE code including solvers, optical properties, fluxes + ! Tests are run on idealized problems with analytic solutions (e.g. radiative equilibrium) + ! Solutions are checked for correctness where possible + ! Some tests check invariance, e.g. with respect to vertical ordering + ! + use mo_rte_kind, only: wp + use mo_optical_props, only: ty_optical_props, & + ty_optical_props_arry, & + ty_optical_props_1scl, ty_optical_props_2str, ty_optical_props_nstr + use mo_rte_util_array, only: zero_array + use mo_fluxes, only: ty_fluxes_broadband + use mo_rte_sw, only: rte_sw + use mo_testing_utils, only: allclose, stop_on_err, report_err, check_fluxes, & + vr, & + increment_with_1scl, increment_with_2str, increment_with_nstr + implicit none + ! ---------------------------------------------------------------------------------- + ! + ! Exercise various paths through RTE SW solvers + ! For the moment tests are run using thin, scattering, gray atmospheres + ! and checked for correctness against known analytic solution (not a great approximation) + ! Beyond correctness tests check invariance, e.g. with respect to vertical ordering + ! Tests include + ! Net fluxes = down-up when computed in various combos (net only, up/down only, all three) + ! using ty_broadband_fluxes + ! Answers are invariant to + ! Extracting subsets + ! Vertical orientation + ! Adding transparent optical properties + ! Shortwave specific tests: + ! Solutions are linear in TOA flux + ! + ! Other possibilites: + ! Vertical discretization? Maybe just check boundary fluxes + ! Test the application of the boundary condition? + + real(wp), parameter :: pi = acos(-1._wp) + integer, parameter :: ncol = 8, nlay = 16 + integer, parameter :: nmu0 = 2 + integer :: icol, ilay, imu0 + + ! + ! Shorteave tests - thin atmosphere + ! + real(wp), dimension(2), parameter :: g = [0.85_wp, 0.65_wp], & + tau = [1.e-4_wp, 1.e-2_wp], & + ssa = 1._wp - & + [1.e-4_wp, 1.e-2_wp] + real(wp), dimension(nmu0), parameter :: mu0 = [1._wp, 0.5_wp] + real(wp), dimension(1,ncol), parameter :: sfc_albedo = 0._wp + real(wp), dimension(ncol,1), parameter :: toa_flux = 1._wp + real(wp), parameter :: factor = 2._wp ! for checking linearity + real(wp), dimension(ncol) :: mu0_arr + + type(ty_optical_props_2str) :: atmos + type(ty_fluxes_broadband) :: fluxes + logical :: top_at_1 + 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 +#ifdef __NVCOMPILER + ! We need the following temporary variables to circumvent the -Mbounds issue: + real(wp), dimension(ncol, nlay+1) :: tst_flux_up_reversed, ref_flux_dn_reversed +#endif + real(wp), dimension(:), pointer :: sfc + + logical :: passed + + ! ------------------------------------------------------------------------------------------------------ + top_at_1 = .true. + ! ------------------------------------------------------------------------------------ + ! + ! Shortwave tests - thin atmospheres + ! + ! ------------------------------------------------------------------------------------ + print *, "RTE SW solver unit tests" + + print *, "Thin, scattering atmospheres" + call stop_on_err(atmos%init(band_lims_wvn = reshape([3250._wp, 10000._wp], shape = [2, 1]), & + band_lims_gpt = reshape([1, 1 ], shape = [2, 1]), & + name = "Gray atmosphere")) + call stop_on_err(atmos%alloc_2str(ncol, nlay)) + call thin_scattering(tau, ssa, g, nlay, atmos) + + do imu0 = 1, nmu0 + print '(" mu0 = ", f4.2)', mu0(imu0) + mu0_arr = mu0(imu0) + fluxes%flux_up => ref_flux_up (:,:) + 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, & + mu0_arr, & + toa_flux, & + sfc_albedo, sfc_albedo, & + fluxes)) + ! + ! Is the solution correct? + ! WIP - differences are up to 25%, so skip correctness test for the moment + ! + ! passed = check_thin_scattering(atmos, spread(mu0(1), 1, ncol), top_at_1, & + ! ref_flux_up, ref_flux_dn, ref_flux_dir) + if(imu0 == 1) passed = .true. + ! ------------------------------------------------------------------------------------ + ! + ! Check direct beam for correctness with Beer-Lambert-Bouguier + ! + if(top_at_1) then + sfc => ref_flux_dir(:,nlay+1) + else + sfc => ref_flux_dir(:, 1) + end if + if(.not. allclose(sfc, & + toa_flux(:,1)*mu0_arr*exp(-sum(atmos%tau(:,:,1),dim=2)/mu0_arr), & + tol=20._wp)) then ! Tolerances as big as 20 needed for GPU implementations + passed = .false. + call report_err("Direct flux doesn't match") + end if + ! ------------------------------------------------------------------------------------ + ! + ! Net fluxes on- vs off-line + ! Are the net fluxes correct? + ! + print *, " Net flux variants" + call check_fluxes(ref_flux_net, ref_flux_dn-ref_flux_up, passed, "net fluxes don't match down-up") + ! + ! Compute only net fluxes + ! + nullify(fluxes%flux_up) + nullify(fluxes%flux_dn) + call stop_on_err(rte_sw(atmos, top_at_1, & + mu0_arr, & + toa_flux, & + sfc_albedo, sfc_albedo, & + fluxes)) + call check_fluxes(ref_flux_net, ref_flux_dn-ref_flux_up, & + passed, "Net fluxes computed alone doesn't match down-up computed separately") + ! + ! Compute only up and down fluxes + ! + fluxes%flux_up => tst_flux_up (:,:) + fluxes%flux_dn => tst_flux_dn (:,:) + call stop_on_err(rte_sw(atmos, top_at_1, & + mu0_arr, & + toa_flux, & + sfc_albedo, sfc_albedo, & + fluxes)) + call check_fluxes(ref_flux_net, tst_flux_dn-tst_flux_up, & + passed, "Net fluxes don't match down-up computed together") + ! ------------------------------------------------------- + ! + ! Subsets of atmospheric columns + ! + print *, " Subsetting invariance" + call clear_sky_subset(atmos, mu0_arr, toa_flux, sfc_albedo, tst_flux_up, tst_flux_dn) + call check_fluxes(tst_flux_up, ref_flux_up, & + tst_flux_dn, ref_flux_dn, & + passed, "Doing problem in subsets fails") + ! ------------------------------------------------------- + ! + ! Vertically-reverse + ! + 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, & + mu0_arr, & + toa_flux, & + sfc_albedo, sfc_albedo, & + fluxes)) +#ifdef __NVCOMPILER + ! The -Mbounds flag breaks the actual arrays that are passed to check_fluxes + tst_flux_up_reversed = tst_flux_up(:,nlay+1:1:-1) + ref_flux_dn_reversed = tst_flux_dn(:,nlay+1:1:-1) + call check_fluxes(tst_flux_up_reversed, ref_flux_up, & + ref_flux_dn_reversed, ref_flux_dn, & +#else + call check_fluxes(tst_flux_up(:,nlay+1:1:-1), ref_flux_up, & + tst_flux_dn(:,nlay+1:1:-1), ref_flux_dn, & +#endif + passed, "Doing problem upside down fails") + call vr(atmos) + ! ------------------------------------------------------- + ! + ! Linear in TOA flux + ! + print *, " Linear in TOA flux" + call thin_scattering(tau, ssa, g, nlay, atmos) + call stop_on_err(rte_sw(atmos, top_at_1, & + mu0_arr, & + toa_flux * factor, & + sfc_albedo, sfc_albedo, & + fluxes)) + call check_fluxes(tst_flux_up/factor, ref_flux_up, & + tst_flux_dn/factor, ref_flux_dn, & + passed, "Linearity in TOA flux fails") + ! ------------------------------------------------------------------------------------ + end do + ! Done + ! + print *, "RTE SW solver unit tests done" + print * + if(.not. passed) error stop 1 + ! ------------------------------------------------------------------------------------ +contains + ! ------------------------------------------------------------------------------------ + ! + ! Define an atmosphere in gray radiative equillibrium + ! See, for example, section 2 of Weaver and Rmanathan 1995 https://doi.org/10.1029/95JD00770 + ! + subroutine thin_scattering(tau, ssa, g, nlay, atmos) + real(wp), dimension(:), intent(in) :: tau, ssa, g + integer, intent(in) :: nlay + type(ty_optical_props_2str), intent(inout) :: atmos + + integer :: ntau, nssa, ng, ncol + integer :: i, j, k + real(wp), dimension(size(tau)*size(ssa)*size(g)) & + :: temp + + ntau = size(tau); nssa = size(ssa); ng = size(g) + ncol = ntau*nssa*ng + if(ncol /= atmos%get_ncol()) call stop_on_err("Number of SW columns incompatible") + ! + ! Set up a gray spectral distribution - one band, one g-point + ! + call stop_on_err(atmos%init(band_lims_wvn = reshape([3250._wp, 1.e5_wp], shape = [2, 1]), & + band_lims_gpt = reshape([1, 1], shape = [2, 1]), & + name = "Gray SW atmosphere")) + call stop_on_err(atmos%alloc_2str(ncol, nlay)) + + temp = [(((tau(i), k = 1, 1 ), i = 1, ntau), j = 1, nssa*ng)] + ! + ! Divide optical depth evenly among layers + ! + atmos%tau(1:ncol,1:nlay,1) = spread(temp(1:ncol)/real(nlay, wp), dim=2, ncopies=nlay) + ! + ! ... and make the medium uniform + ! + temp = [(((ssa(i), k = 1, ntau ), i = 1, nssa), j = 1, ng)] + atmos%ssa(1:ncol,1:nlay,1) = spread(temp(1:ncol), dim=2, ncopies=nlay) + temp = [(((g (i), k = 1, ntau*ng), i = 1, ng ), j = 1, 1 )] + atmos%g (1:ncol,1:nlay,1) = spread(temp(1:ncol), dim=2, ncopies=nlay) + + if(.false.) then + print *, "Original values" + print '("tau: ", 8(e9.3,2x))', sum(atmos%tau(:,:,1),dim=2) + print '("ssa: ", 8(e9.3,2x))', atmos%ssa(:,1,1) + print '("g : ", 8(e9.3,2x))', atmos%g (:,1,1) + print * + end if + ! + ! Delta-scale + ! + call stop_on_err(atmos%delta_scale()) + + end subroutine thin_scattering + ! ------------------------------------------------------------------------------------ + function check_thin_scattering(atmos, mu0, top_at_1, ref_flux_up, ref_flux_dn, ref_flux_dir) + type(ty_optical_props_2str), intent(in) :: atmos + real(wp), dimension(:), intent(in) :: mu0 + logical, intent(in) :: top_at_1 + real(wp), dimension(:,:), intent(in) :: ref_flux_up, ref_flux_dn, ref_flux_dir + logical :: check_thin_scattering + + real(wp), dimension(atmos%get_ncol()) :: gamma3, R, T ! Reflectance, transmittance + + check_thin_scattering = .true. + ! + ! Solutions for small tau + ! Meador and Weaver 1980, 10.1175/1520-0469(1980)037<0630:TSATRT>2.0.CO;2 + ! Equation 19 using the same gamma3 as in RTE (and gamma3+gamma4=1) + ! ssa and g are assumed to be vertically uniform + ! + gamma3(:) = (2._wp - 3._wp * mu0(:) * atmos%g(:,1,1)) * .25_wp + R(:) = (atmos%ssa(:,1,1)*sum(atmos%tau(:,:,1),dim=2))/mu0(:) * gamma3(:) + if(.false.) then + print '("tau: ", 8(e9.3,2x))', sum(atmos%tau(:,:,1),dim=2) + print '("ssa: ", 8(e9.3,2x))', atmos%ssa(:,1,1) + print '("g : ", 8(e9.3,2x))', atmos%g (:,1,1) + print * + print '("R : ", 8(e9.3,2x))', R + print '("RTE: ", 8(e9.3,2x))', ref_flux_up(:,1) + print '("Dif: ", 8(e9.3,2x))', R(:) - ref_flux_up(:,1) + print '("Rel: ", 8(f9.2,2x))', (R(:) - ref_flux_up(:,1))/ref_flux_up(:,1) * 100._wp + end if + end function check_thin_scattering + ! ------------------------------------------------------------------------------------ + ! + ! Invariance tests + ! + ! ------------------------------------------------------------------------------------ + ! + ! Clear-sky fluxes, half the columns at a time + ! We're counting on ncol being even + ! + subroutine clear_sky_subset(atmos, mu0, toa_flux, sfc_albedo, flux_up, flux_dn) + type(ty_optical_props_2str), intent(inout) :: atmos + real(wp), dimension(:), intent(in ) :: mu0 + real(wp), dimension(:,:), intent(in ) :: toa_flux, sfc_albedo + real(wp), dimension(:,:), intent( out) :: flux_up, flux_dn + + type(ty_optical_props_2str) :: atmos_subset + type(ty_fluxes_broadband) :: fluxes ! Use local variable + real(wp), dimension(atmos%get_ncol()/2, & + atmos%get_nlay()+1), target & + :: up, dn + integer :: i, colS, colE + integer :: ncol + ! ------------------------------ + ncol = atmos%get_ncol() + call stop_on_err(atmos_subset%init(atmos)) + fluxes%flux_up => up + fluxes%flux_dn => dn + + do i = 1, 2 + colS = ((i-1) * ncol/2) + 1 + colE = i * ncol/2 + call stop_on_err(atmos%get_subset(colS, ncol/2, atmos_subset)) + call stop_on_err(rte_sw(atmos_subset, top_at_1, & + mu0(colS:colE), & + toa_flux(colS:colE,:), & + sfc_albedo(:,colS:colE), sfc_albedo(:,colS:colE), & + fluxes)) + flux_up(colS:colE,:) = up + flux_dn(colS:colE,:) = dn + end do + end subroutine clear_sky_subset + ! ------------------------------------------------------------------------------------ + +end program rte_sw_solver_unit_tests \ No newline at end of file diff --git a/tests/test_zenith_angle_spherical_correction.F90 b/tests/test_zenith_angle_spherical_correction.F90 index 851281521..26ad615e0 100644 --- a/tests/test_zenith_angle_spherical_correction.F90 +++ b/tests/test_zenith_angle_spherical_correction.F90 @@ -16,11 +16,10 @@ program test_solar_zenith_angle use mo_rcemip_profiles, only: make_rcemip_profiles use mo_gas_concentrations, only: ty_gas_concs - use mo_gas_optics_rrtmgp, only: ty_gas_optics_rrtmgp + use mo_gas_optics_defs, only: gas_optics, load_and_init use mo_optical_props, only: ty_optical_props_2str use mo_fluxes, only: ty_fluxes_broadband use mo_rte_sw, only: rte_sw - use mo_load_coefficients, only: load_and_init use mo_heating_rates, only: compute_heating_rate implicit none @@ -31,7 +30,6 @@ program test_solar_zenith_angle real(wp), dimension(ncol, 0:nlay) :: alt, mu type(ty_gas_concs) :: gas_concs - type(ty_gas_optics_rrtmgp) :: gas_optics type(ty_optical_props_2str) :: atmos type(ty_fluxes_broadband) :: fluxes real(wp) :: p_lay(ncol, nlay), t_lay(ncol, nlay), & From 3274155d41bfc51dc71ea91672e84f6219080038 Mon Sep 17 00:00:00 2001 From: Sergey Kosukhin Date: Mon, 19 Feb 2024 18:00:55 +0100 Subject: [PATCH 12/20] Fix cce-gpu-openacc-SP --- .gitlab/lumi.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.gitlab/lumi.yml b/.gitlab/lumi.yml index 6899ac95d..2b105d19c 100644 --- a/.gitlab/lumi.yml +++ b/.gitlab/lumi.yml @@ -126,5 +126,5 @@ cce-gpu-openacc-DP: cce-gpu-openacc-SP: extends: - - .dp + - .sp - .cce-gpu-openacc From e24246637482d2352114537c1c7c1405296e6dfd Mon Sep 17 00:00:00 2001 From: Sergey Kosukhin Date: Mon, 19 Feb 2024 18:38:33 +0100 Subject: [PATCH 13/20] Small refactoring for the build system (#266) Streamline and add flexibility to the build process. --- .github/workflows/containerized-ci.yml | 10 +++---- .github/workflows/continuous-integration.yml | 8 +++--- .github/workflows/self-hosted-ci.yml | 10 +++---- .gitlab/levante.yml | 8 +++--- .gitlab/lumi.yml | 8 +++--- Makefile | 28 ++++++++++---------- examples/all-sky/Makefile | 4 +-- examples/rfmip-clear-sky/Makefile | 2 +- tests/Makefile | 2 +- 9 files changed, 40 insertions(+), 40 deletions(-) diff --git a/.github/workflows/containerized-ci.yml b/.github/workflows/containerized-ci.yml index 01343d395..a43718cf1 100644 --- a/.github/workflows/containerized-ci.yml +++ b/.github/workflows/containerized-ci.yml @@ -80,8 +80,8 @@ jobs: if: matrix.fortran-compiler != 'ifx' || matrix.rte-kernels != 'accel' run: | $FC --version - make libs - make -C build separate-libs + make -j4 libs + make -j4 -C build separate-libs # # Build libraries, examples and tests (expect failure) # @@ -90,7 +90,7 @@ jobs: shell: bash run: | $FC --version - make libs 2> >(tee make.err >&2) && { + make -j4 libs 2> >(tee make.err >&2) && { echo "Unexpected success" exit 1 } || { @@ -106,7 +106,7 @@ jobs: # - name: Run examples and tests if: steps.build-success.outcome != 'skipped' - run: make tests + run: make -j4 tests # # Relax failure thresholds for single precision # @@ -118,7 +118,7 @@ jobs: # - name: Compare the results if: steps.build-success.outcome != 'skipped' - run: make check + run: make -j4 check # # Generate validation plots # diff --git a/.github/workflows/continuous-integration.yml b/.github/workflows/continuous-integration.yml index 2c5a931bb..ed8ccd964 100644 --- a/.github/workflows/continuous-integration.yml +++ b/.github/workflows/continuous-integration.yml @@ -89,15 +89,15 @@ jobs: - name: Build libraries, examples and tests run: | $FC --version - make libs - make -C build separate-libs + make -j4 libs + make -j4 -C build separate-libs # # Run examples and tests # - name: Run examples and tests - run: make tests + run: make -j4 tests # # Compare the results # - name: Compare the results - run: make check + run: make -j4 check diff --git a/.github/workflows/self-hosted-ci.yml b/.github/workflows/self-hosted-ci.yml index 87dd7e32d..3e4810eff 100644 --- a/.github/workflows/self-hosted-ci.yml +++ b/.github/workflows/self-hosted-ci.yml @@ -96,22 +96,22 @@ jobs: - name: Build libraries, examples and tests run: | $FC --version - make libs - make -C build separate-libs + make -j8 libs + make -j8 -C build separate-libs # # Run examples and tests (expect success) # - name: Run examples and tests (expect success) id: run-success if: matrix.config-name != 'cce-gpu-openmp' - run: make tests + run: make -j8 tests # # Run examples and tests (expect failure) # - name: Run examples and tests (expect failure) if: steps.run-success.outcome == 'skipped' run: | - make tests && { + make -j8 tests && { echo "Unexpected success" exit 1 } || echo "Expected failure" @@ -126,4 +126,4 @@ jobs: # - name: Compare the results if: steps.run-success.outcome != 'skipped' - run: make check + run: make -j8 check diff --git a/.gitlab/levante.yml b/.gitlab/levante.yml index 8104acd52..78de62697 100644 --- a/.gitlab/levante.yml +++ b/.gitlab/levante.yml @@ -77,8 +77,8 @@ variables: # Build libraries, examples and tests # - ${FC} ${VERSION_FCFLAGS} - - make libs - - make -C build separate-libs + - make -j8 libs + - make -j8 -C build separate-libs # # Check out data # @@ -86,11 +86,11 @@ variables: # # Run examples and tests # - - make tests + - make -j8 tests # # Compare the results # - - make check + - make -j8 check .nvhpc-gpu-openacc: extends: diff --git a/.gitlab/lumi.yml b/.gitlab/lumi.yml index 6899ac95d..d8894c1be 100644 --- a/.gitlab/lumi.yml +++ b/.gitlab/lumi.yml @@ -92,8 +92,8 @@ setup-python: # Build libraries, examples and tests # - ${FC} ${VERSION_FCFLAGS} - - make libs - - make -C build separate-libs + - make -j8 libs + - make -j8 -C build separate-libs # # Check out data # @@ -101,11 +101,11 @@ setup-python: # # Run examples and tests # - - make tests + - make -j8 tests # # Compare the results # - - make check + - make -j8 check .cce-gpu-openacc: extends: diff --git a/Makefile b/Makefile index 668b88f93..ec6eb75a3 100644 --- a/Makefile +++ b/Makefile @@ -5,27 +5,27 @@ all: libs tests check docs libs: - make -C build -j - make -C tests -j 1 - make -C examples/all-sky -j - make -C examples/rfmip-clear-sky -j + $(MAKE) -C build + $(MAKE) -C tests + $(MAKE) -C examples/all-sky + $(MAKE) -C examples/rfmip-clear-sky tests: - make -C tests tests - make -C examples/rfmip-clear-sky tests - make -C examples/all-sky tests + $(MAKE) -C examples/rfmip-clear-sky $@ + $(MAKE) -C examples/all-sky $@ + $(MAKE) -C tests $@ check: - make -C tests check - make -C examples/rfmip-clear-sky check - make -C examples/all-sky check + $(MAKE) -C examples/rfmip-clear-sky $@ + $(MAKE) -C examples/all-sky $@ + $(MAKE) -C tests $@ docs: @cd doc; ./build_documentation.sh clean: - make -C build clean - make -C tests clean - make -C examples/rfmip-clear-sky clean - make -C examples/all-sky clean + $(MAKE) -C build $@ + $(MAKE) -C examples/rfmip-clear-sky $@ + $(MAKE) -C examples/all-sky $@ + $(MAKE) -C tests $@ rm -rf public diff --git a/examples/all-sky/Makefile b/examples/all-sky/Makefile index b221dd0a2..a3da6beb1 100644 --- a/examples/all-sky/Makefile +++ b/examples/all-sky/Makefile @@ -50,10 +50,10 @@ tests: $(RUN_CMD) bash all_tests.sh check: - python ${RRTMGP_ROOT}/examples/compare-to-reference.py --ref_dir ${RRTMGP_DATA}/examples/all-sky/reference --tst_dir ${RRTMGP_ROOT}/examples/all-sky \ + $${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 ${RRTMGP_ROOT}/examples/compare-to-reference.py --ref_dir ${RRTMGP_DATA}/examples/all-sky/reference --tst_dir ${RRTMGP_ROOT}/examples/all-sky \ + $${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 diff --git a/examples/rfmip-clear-sky/Makefile b/examples/rfmip-clear-sky/Makefile index 6f2bc3621..0d90cea7a 100644 --- a/examples/rfmip-clear-sky/Makefile +++ b/examples/rfmip-clear-sky/Makefile @@ -57,7 +57,7 @@ tests: multiple_input4MIPs_radiation_RFMIP_UColorado-RFMIP-1-2_none.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 ${RRTMGP_ROOT}/examples/compare-to-reference.py \ + $${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 diff --git a/tests/Makefile b/tests/Makefile index aaa159f47..1b90ed26c 100644 --- a/tests/Makefile +++ b/tests/Makefile @@ -66,7 +66,7 @@ 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.o: $(LIB_DEPS) mo_simple_netcdf.o mo_rfmip_io.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 From dc88744316e9eb0620f5d9c4ccab7d7e4f959f8f Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Wed, 28 Feb 2024 10:21:33 -0500 Subject: [PATCH 14/20] Make front-end/kernel libraries by default --- .github/workflows/containerized-ci.yml | 1 - .github/workflows/continuous-integration.yml | 1 - .github/workflows/self-hosted-ci.yml | 1 - .gitlab/levante.yml | 1 - .gitlab/lumi.yml | 1 - build/Makefile | 3 ++- 6 files changed, 2 insertions(+), 6 deletions(-) diff --git a/.github/workflows/containerized-ci.yml b/.github/workflows/containerized-ci.yml index a43718cf1..0f247ba05 100644 --- a/.github/workflows/containerized-ci.yml +++ b/.github/workflows/containerized-ci.yml @@ -81,7 +81,6 @@ jobs: run: | $FC --version make -j4 libs - make -j4 -C build separate-libs # # Build libraries, examples and tests (expect failure) # diff --git a/.github/workflows/continuous-integration.yml b/.github/workflows/continuous-integration.yml index ed8ccd964..3cb19a757 100644 --- a/.github/workflows/continuous-integration.yml +++ b/.github/workflows/continuous-integration.yml @@ -90,7 +90,6 @@ jobs: run: | $FC --version make -j4 libs - make -j4 -C build separate-libs # # Run examples and tests # diff --git a/.github/workflows/self-hosted-ci.yml b/.github/workflows/self-hosted-ci.yml index 3e4810eff..21fe2116e 100644 --- a/.github/workflows/self-hosted-ci.yml +++ b/.github/workflows/self-hosted-ci.yml @@ -97,7 +97,6 @@ jobs: run: | $FC --version make -j8 libs - make -j8 -C build separate-libs # # Run examples and tests (expect success) # diff --git a/.gitlab/levante.yml b/.gitlab/levante.yml index 78de62697..4835a4383 100644 --- a/.gitlab/levante.yml +++ b/.gitlab/levante.yml @@ -78,7 +78,6 @@ variables: # - ${FC} ${VERSION_FCFLAGS} - make -j8 libs - - make -j8 -C build separate-libs # # Check out data # diff --git a/.gitlab/lumi.yml b/.gitlab/lumi.yml index 997ccf4b9..245129d5c 100644 --- a/.gitlab/lumi.yml +++ b/.gitlab/lumi.yml @@ -93,7 +93,6 @@ setup-python: # - ${FC} ${VERSION_FCFLAGS} - make -j8 libs - - make -j8 -C build separate-libs # # Check out data # diff --git a/build/Makefile b/build/Makefile index d00adcb3a..800bbed3e 100644 --- a/build/Makefile +++ b/build/Makefile @@ -9,7 +9,8 @@ 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 +all: librte.a librrtmgp.a \ + librtekernels.a librtef.a librrtmgpkernels.a librrtmgpf.a separate-libs: librtekernels.a librtef.a librrtmgpkernels.a librrtmgpf.a COMPILE = $(FC) $(FCFLAGS) $(FCINCLUDE) -c From f6d4a63118f535480def9319661afc3c163ae27a Mon Sep 17 00:00:00 2001 From: Jian Sun Date: Fri, 22 Mar 2024 11:41:37 -0600 Subject: [PATCH 15/20] fix the gpu runtime error on nvidia a100 with nvhpc/24.3 --- rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 | 22 ++++++++++++++++------ 1 file changed, 16 insertions(+), 6 deletions(-) diff --git a/rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 b/rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 index c83601586..906c7cccb 100644 --- a/rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 +++ b/rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 @@ -801,6 +801,7 @@ function set_tsi(this, tsi) result(error_msg) character(len=128) :: error_msg !! Empty if successful real(wp) :: norm + integer :: igpt, length ! ---------------------------------------------------------- error_msg = "" if(tsi < 0._wp) then @@ -809,12 +810,21 @@ function set_tsi(this, tsi) result(error_msg) ! ! Scale the solar source function to the input tsi ! - !$acc kernels - !$omp target - norm = 1._wp/sum(this%solar_source(:)) - this%solar_source(:) = this%solar_source(:) * tsi * norm - !$acc end kernels - !$omp end target + norm = 0._wp + length = size(this%solar_source) + !$acc parallel loop gang vector reduction(+:norm) + !$omp target teams distribute parallel do simd reduction(+:norm) + do igpt = 1, length + norm = norm + this%solar_source(igpt) + end do + + norm = 1._wp/norm + + !$acc parallel loop gang vector + !$omp target teams distribute parallel do simd + do igpt = 1, length + this%solar_source(igpt) = this%solar_source(igpt) * tsi * norm + end do end if end function set_tsi From 36f34c6d5bb4c6462eb31473ca151b0520a1c252 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Tue, 2 Apr 2024 10:25:06 -0400 Subject: [PATCH 16/20] Add kernel API (#272) Adds C and Fortran headers for the kernels in *-kernels/api subdirectories. Kernel headers are tested for consistency against Fortran front-end. Changes the makefiles a) to not build the documentation by defaults, and b) to group the building and running of tests, so that make libs makes, um, only the libraries. --------- Co-authored-by: Alexander Soklev --- .github/workflows/check-api.yml | 34 ++ .github/workflows/continuous-integration.yml | 4 +- .github/workflows/self-hosted-ci.yml | 8 +- Makefile | 15 +- build/Makefile | 26 +- examples/all-sky/Makefile | 3 +- examples/rfmip-clear-sky/Makefile | 13 +- extensions/mo_fluxes_byband.F90 | 4 +- .../api/mo_gas_optics_rrtmgp_kernels.F90 | 245 +++++++++++ rrtmgp-kernels/api/rrtmgp_kernels.h | 123 ++++++ .../api/mo_fluxes_broadband_kernels.F90 | 74 ++++ rte-kernels/api/mo_optical_props_kernels.F90 | 377 +++++++++++++++++ rte-kernels/api/mo_rte_solver_kernels.F90 | 202 +++++++++ rte-kernels/api/mo_rte_util_array.F90 | 47 +++ rte-kernels/api/rte_kernels.h | 389 ++++++++++++++++++ rte-kernels/api/rte_types.h | 34 ++ tests/Makefile | 13 +- 17 files changed, 1575 insertions(+), 36 deletions(-) create mode 100644 .github/workflows/check-api.yml create mode 100644 rrtmgp-kernels/api/mo_gas_optics_rrtmgp_kernels.F90 create mode 100644 rrtmgp-kernels/api/rrtmgp_kernels.h create mode 100644 rte-kernels/api/mo_fluxes_broadband_kernels.F90 create mode 100644 rte-kernels/api/mo_optical_props_kernels.F90 create mode 100644 rte-kernels/api/mo_rte_solver_kernels.F90 create mode 100644 rte-kernels/api/mo_rte_util_array.F90 create mode 100644 rte-kernels/api/rte_kernels.h create mode 100644 rte-kernels/api/rte_types.h diff --git a/.github/workflows/check-api.yml b/.github/workflows/check-api.yml new file mode 100644 index 000000000..b3681e666 --- /dev/null +++ b/.github/workflows/check-api.yml @@ -0,0 +1,34 @@ +name: Check Fortran API +on: + push: + branches: + - main + - develop + pull_request: + branches-ignore: + - documentation + workflow_dispatch: + + +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 }} + RTE_KERNELS: extern + steps: + # + # Check out repository under $GITHUB_WORKSPACE + # + - name: Check out code + uses: actions/checkout@v4 + # + # Build libraries + # + - name: Build libraries + run: | + $FC --version + make -j4 libs \ No newline at end of file diff --git a/.github/workflows/continuous-integration.yml b/.github/workflows/continuous-integration.yml index 3cb19a757..138994fd4 100644 --- a/.github/workflows/continuous-integration.yml +++ b/.github/workflows/continuous-integration.yml @@ -86,14 +86,14 @@ jobs: # # Build libraries, examples and tests # - - name: Build libraries, examples and tests + - name: Build libraries run: | $FC --version make -j4 libs # # Run examples and tests # - - name: Run examples and tests + - name: Build and run examples and tests run: make -j4 tests # # Compare the results diff --git a/.github/workflows/self-hosted-ci.yml b/.github/workflows/self-hosted-ci.yml index 21fe2116e..5ef3cd98b 100644 --- a/.github/workflows/self-hosted-ci.yml +++ b/.github/workflows/self-hosted-ci.yml @@ -42,7 +42,7 @@ jobs: env: # Core variables: FC: ftn - FCFLAGS: ${{ matrix.fcflags }} -DRTE_USE_${{ matrix.fpmodel}} + FCFLAGS: ${{ matrix.fcflags }} -DRTE_USE_${{ matrix.fpmodel }} # Make variables: RRTMGP_ROOT: ${{ github.workspace }} RRTMGP_DATA: ${{ github.workspace }}/rrtmgp-data @@ -93,21 +93,21 @@ jobs: # # Build libraries, examples and tests # - - name: Build libraries, examples and tests + - name: Build libraries run: | $FC --version make -j8 libs # # Run examples and tests (expect success) # - - name: 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: 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 && { diff --git a/Makefile b/Makefile index ec6eb75a3..96c1e2edd 100644 --- a/Makefile +++ b/Makefile @@ -1,31 +1,28 @@ # # Top-level Makefile # -.PHONY: libs tests check docs -all: libs tests check docs +.PHONY: libs tests check +all: libs tests check libs: - $(MAKE) -C build - $(MAKE) -C tests - $(MAKE) -C examples/all-sky - $(MAKE) -C examples/rfmip-clear-sky + $(MAKE) -C build $@ tests: + $(MAKE) -C tests $@ $(MAKE) -C examples/rfmip-clear-sky $@ $(MAKE) -C examples/all-sky $@ - $(MAKE) -C tests $@ check: + $(MAKE) -C tests $@ $(MAKE) -C examples/rfmip-clear-sky $@ $(MAKE) -C examples/all-sky $@ - $(MAKE) -C tests $@ docs: @cd doc; ./build_documentation.sh clean: $(MAKE) -C build $@ + $(MAKE) -C tests $@ $(MAKE) -C examples/rfmip-clear-sky $@ $(MAKE) -C examples/all-sky $@ - $(MAKE) -C tests $@ rm -rf public diff --git a/build/Makefile b/build/Makefile index 800bbed3e..efd0341d2 100644 --- a/build/Makefile +++ b/build/Makefile @@ -12,6 +12,7 @@ RRTMGP_KERNEL_DIR = ../rrtmgp-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 @@ -22,23 +23,38 @@ include $(RRTMGP_DIR)/Make.depends include $(RTE_KERNEL_DIR)/Make.depends include $(RRTMGP_KERNEL_DIR)/Make.depends -VPATH = $(RTE_DIR):$(RTE_KERNEL_DIR):$(RRTMGP_DIR):$(RRTMGP_KERNEL_DIR):$(GAS_OPTICS_DIR) # # If using OpenACC/OpenMP files in *-kernels/accel take precendence # ifeq ($(RTE_KERNELS), accel) - VPATH = $(RTE_DIR):$(RTE_KERNEL_DIR)/accel:$(RTE_KERNEL_DIR):$(RRTMGP_DIR):$(RRTMGP_KERNEL_DIR)/accel:$(RRTMGP_KERNEL_DIR):$(GAS_OPTICS_DIR) + 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) diff --git a/examples/all-sky/Makefile b/examples/all-sky/Makefile index a3da6beb1..d404181c1 100644 --- a/examples/all-sky/Makefile +++ b/examples/all-sky/Makefile @@ -1,3 +1,4 @@ +#!/usr/bin/env make # # Location of RTE+RRTMGP libraries, module files. # @@ -46,7 +47,7 @@ mo_load_coefficients.o: mo_simple_netcdf.o 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: +tests: rrtmgp_allsky $(RUN_CMD) bash all_tests.sh check: diff --git a/examples/rfmip-clear-sky/Makefile b/examples/rfmip-clear-sky/Makefile index 0d90cea7a..bad1bc80f 100644 --- a/examples/rfmip-clear-sky/Makefile +++ b/examples/rfmip-clear-sky/Makefile @@ -1,3 +1,4 @@ +#!/usr/bin/env make # # Location of RTE+RRTMGP libraries, module files. # @@ -5,11 +6,10 @@ RRTMGP_BUILD = $(RRTMGP_ROOT)/build # # RRTMGP library, module files # -# LDFLAGS += -L$(RRTMGP_BUILD) -# LIBS += -lrrtmgp -lrte +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 @@ -38,11 +38,11 @@ ADDITIONS = mo_simple_netcdf.o mo_rfmip_io.o mo_load_coefficients.o all: rrtmgp_rfmip_lw rrtmgp_rfmip_sw -rrtmgp_rfmip_lw: rrtmgp_rfmip_lw.o $(ADDITIONS) $(RRTMGP_BUILD)/librrtmgp.a $(RRTMGP_BUILD)/librte.a +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_BUILD)/librrtmgp.a $(RRTMGP_BUILD)/librte.a +rrtmgp_rfmip_sw: rrtmgp_rfmip_sw.o $(ADDITIONS) rrtmgp_rfmip_sw.o: rrtmgp_rfmip_sw.F90 $(ADDITIONS) @@ -50,7 +50,8 @@ 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: multiple_input4MIPs_radiation_RFMIP_UColorado-RFMIP-1-2_none.nc \ +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 diff --git a/extensions/mo_fluxes_byband.F90 b/extensions/mo_fluxes_byband.F90 index 587ab42bf..7bf334b66 100644 --- a/extensions/mo_fluxes_byband.F90 +++ b/extensions/mo_fluxes_byband.F90 @@ -156,7 +156,7 @@ end function are_desired_byband ! ! Spectral reduction over all points ! - subroutine sum_byband(ncol, nlev, ngpt, nbnd, band_lims, spectral_flux, byband_flux) bind (C) + subroutine sum_byband(ncol, nlev, ngpt, nbnd, band_lims, spectral_flux, byband_flux) bind(C, name="rte_sum_byband") integer, intent(in ) :: ncol, nlev, ngpt, nbnd integer, dimension(2, nbnd), intent(in ) :: band_lims real(wp), dimension(ncol, nlev, ngpt), intent(in ) :: spectral_flux @@ -181,7 +181,7 @@ end subroutine sum_byband ! ! Net flux: Spectral reduction over all points ! - subroutine net_byband_full(ncol, nlev, ngpt, nbnd, band_lims, spectral_flux_dn, spectral_flux_up, byband_flux_net) bind (C) + subroutine net_byband_full(ncol, nlev, ngpt, nbnd, band_lims, spectral_flux_dn, spectral_flux_up, byband_flux_net) bind(C, name="rte_net_byband_full") integer, intent(in ) :: ncol, nlev, ngpt, nbnd integer, dimension(2, nbnd), intent(in ) :: band_lims real(wp), dimension(ncol, nlev, ngpt), intent(in ) :: spectral_flux_dn, spectral_flux_up diff --git a/rrtmgp-kernels/api/mo_gas_optics_rrtmgp_kernels.F90 b/rrtmgp-kernels/api/mo_gas_optics_rrtmgp_kernels.F90 new file mode 100644 index 000000000..9272bb16e --- /dev/null +++ b/rrtmgp-kernels/api/mo_gas_optics_rrtmgp_kernels.F90 @@ -0,0 +1,245 @@ +module mo_gas_optics_rrtmgp_kernels + use mo_rte_kind, only : wp, wl + use mo_rte_util_array,only : zero_array + implicit none + private + public :: interpolation, compute_tau_absorption, compute_tau_rayleigh, compute_Planck_source + ! ------------------------------------------------------------------------------------------------------------------ + interface + subroutine interpolation( & + ncol,nlay,ngas,nflav,neta, npres, ntemp, & + flavor, & + press_ref_log, temp_ref,press_ref_log_delta, & + temp_ref_min,temp_ref_delta,press_ref_trop_log, & + vmr_ref, & + play,tlay,col_gas, & + jtemp,fmajor,fminor,col_mix,tropo,jeta,jpress) bind(C, name="rrtmgp_interpolation") + use mo_rte_kind, only : wp, wl + ! input dimensions + integer, intent(in) :: ncol,nlay + !! physical domain size + integer, intent(in) :: ngas,nflav,neta,npres,ntemp + !! k-distribution table dimensions + integer, dimension(2,nflav), intent(in) :: flavor + !! index into vmr_ref of major gases for each flavor + real(wp), dimension(npres), intent(in) :: press_ref_log + !! log of pressure dimension in RRTMGP tables + real(wp), dimension(ntemp), intent(in) :: temp_ref + !! temperature dimension in RRTMGP tables + real(wp), intent(in) :: press_ref_log_delta, & + temp_ref_min, temp_ref_delta, & + press_ref_trop_log + !! constants related to RRTMGP tables + real(wp), dimension(2,0:ngas,ntemp), intent(in) :: vmr_ref + !! reference volume mixing ratios used in compute "binary species parameter" eta + + ! inputs from profile or parent function + real(wp), dimension(ncol,nlay), intent(in) :: play, tlay + !! input pressure (Pa?) and temperature (K) + real(wp), dimension(ncol,nlay,0:ngas), intent(in) :: col_gas + !! input column gas amount - molecules/cm^2 + ! outputs + integer, dimension(ncol,nlay), intent(out) :: jtemp, jpress + !! temperature and pressure interpolation indexes + 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 >= 2021 + ! 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 + 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 + !! Interpolation fraction in eta, strat/trop + end subroutine interpolation + end interface + ! ------------------------------------------------------------------------------------------------------------------ + interface + subroutine compute_tau_absorption( & + ncol,nlay,nbnd,ngpt, & ! dimensions + ngas,nflav,neta,npres,ntemp, & + nminorlower, nminorklower, & ! number of minor contributors, total num absorption coeffs + nminorupper, nminorkupper, & + idx_h2o, & + gpoint_flavor, & + band_lims_gpt, & + kmajor, & + kminor_lower, & + kminor_upper, & + minor_limits_gpt_lower, & + minor_limits_gpt_upper, & + minor_scales_with_density_lower, & + minor_scales_with_density_upper, & + scale_by_complement_lower, & + scale_by_complement_upper, & + idx_minor_lower, & + idx_minor_upper, & + idx_minor_scaling_lower, & + idx_minor_scaling_upper, & + kminor_start_lower, & + kminor_start_upper, & + tropo, & + col_mix,fmajor,fminor, & + play,tlay,col_gas, & + jeta,jtemp,jpress, & + tau) bind(C, name="rrtmgp_compute_tau_absorption") + ! --------------------- + use mo_rte_kind, only : wp, wl + ! input dimensions + integer, intent(in) :: ncol,nlay,nbnd,ngpt !! array sizes + integer, intent(in) :: ngas,nflav,neta,npres,ntemp !! tables sizes + integer, intent(in) :: nminorlower, nminorklower,nminorupper, nminorkupper + !! table sizes + integer, intent(in) :: idx_h2o !! index of water vapor in col_gas + ! --------------------- + ! inputs from object + integer, dimension(2,ngpt), intent(in) :: gpoint_flavor + !! major gas flavor (pair) by upper/lower, g-point + integer, dimension(2,nbnd), intent(in) :: band_lims_gpt + !! beginning and ending g-point for each band + real(wp), dimension(ntemp,neta,npres+1,ngpt), intent(in) :: kmajor + !! absorption coefficient table - major gases + real(wp), dimension(ntemp,neta,nminorklower), intent(in) :: kminor_lower + !! absorption coefficient table - minor gases, lower atmosphere + real(wp), dimension(ntemp,neta,nminorkupper), intent(in) :: kminor_upper + !! absorption coefficient table - minor gases, upper atmosphere + integer, dimension(2,nminorlower), intent(in) :: minor_limits_gpt_lower + !! beginning and ending g-point for each minor gas + integer, dimension(2,nminorupper), intent(in) :: minor_limits_gpt_upper + logical(wl), dimension( nminorlower), intent(in) :: minor_scales_with_density_lower + !! generic treatment of minor gases - scales with density (e.g. continuum, collision-induced absorption)? + logical(wl), dimension( nminorupper), intent(in) :: minor_scales_with_density_upper + logical(wl), dimension( nminorlower), intent(in) :: scale_by_complement_lower + !! generic treatment of minor gases - scale by density (e.g. self-continuum) or complement? + logical(wl), dimension( nminorupper), intent(in) :: scale_by_complement_upper + integer, dimension( nminorlower), intent(in) :: idx_minor_lower + !! index of each minor gas in col_gas + integer, dimension( nminorupper), intent(in) :: idx_minor_upper + integer, dimension( nminorlower), intent(in) :: idx_minor_scaling_lower + !! for this minor gas, index of the "scaling gas" in col_gas + integer, dimension( nminorupper), intent(in) :: idx_minor_scaling_upper + integer, dimension( nminorlower), intent(in) :: kminor_start_lower + !! starting g-point index in minor gas absorption table + integer, dimension( nminorupper), intent(in) :: kminor_start_upper + logical(wl), dimension(ncol,nlay), intent(in) :: tropo + !! use upper- or lower-atmospheric tables? + ! --------------------- + ! inputs from profile or parent function + real(wp), dimension(2, ncol,nlay,nflav ), intent(in) :: col_mix + !! combination of major species's column amounts - computed in interpolation() + real(wp), dimension(2,2,2,ncol,nlay,nflav ), intent(in) :: fmajor + !! interpolation weights for major gases - computed in interpolation() + real(wp), dimension(2,2, ncol,nlay,nflav ), intent(in) :: fminor + !! interpolation weights for minor gases - computed in interpolation() + real(wp), dimension( ncol,nlay ), intent(in) :: play, tlay + !! input temperature and pressure + real(wp), dimension( ncol,nlay,0:ngas), intent(in) :: col_gas + !! input column gas amount (molecules/cm^2) + integer, dimension(2, ncol,nlay,nflav ), intent(in) :: jeta + !! interpolation indexes in eta - computed in interpolation() + integer, dimension( ncol,nlay ), intent(in) :: jtemp + !! interpolation indexes in temperature - computed in interpolation() + integer, dimension( ncol,nlay ), intent(in) :: jpress + !! interpolation indexes in pressure - computed in interpolation() + ! --------------------- + ! output - optical depth + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau !! aborption optional depth + end subroutine compute_tau_absorption + end interface + ! ------------------------------------------------------------------------------------------------------------------ + interface + subroutine compute_tau_rayleigh(ncol,nlay,nbnd,ngpt, & + ngas,nflav,neta,npres,ntemp, & + gpoint_flavor,band_lims_gpt, & + krayl, & + idx_h2o, col_dry,col_gas, & + fminor,jeta,tropo,jtemp, & + tau_rayleigh) bind(C, name="rrtmgp_compute_tau_rayleigh") + use mo_rte_kind, only : wp, wl + integer, intent(in ) :: ncol,nlay,nbnd,ngpt + !! input dimensions + integer, intent(in ) :: ngas,nflav,neta,npres,ntemp + !! table dimensions + integer, dimension(2,ngpt), intent(in ) :: gpoint_flavor + !! major gas flavor (pair) by upper/lower, g-point + integer, dimension(2,nbnd), intent(in ) :: band_lims_gpt + !! start and end g-point for each band + real(wp), dimension(ntemp,neta,ngpt,2), intent(in ) :: krayl + !! Rayleigh scattering coefficients + integer, intent(in ) :: idx_h2o + !! index of water vapor in col_gas + real(wp), dimension(ncol,nlay), intent(in ) :: col_dry + !! column amount of dry air + real(wp), dimension(ncol,nlay,0:ngas), intent(in ) :: col_gas + !! input column gas amount (molecules/cm^2) + real(wp), dimension(2,2,ncol,nlay,nflav), intent(in ) :: fminor + !! interpolation weights for major gases - computed in interpolation() + integer, dimension(2, ncol,nlay,nflav), intent(in ) :: jeta + !! interpolation indexes in eta - computed in interpolation() + logical(wl), dimension(ncol,nlay), intent(in ) :: tropo + !! use upper- or lower-atmospheric tables? + integer, dimension(ncol,nlay), intent(in ) :: jtemp + !! interpolation indexes in temperature - computed in interpolation() + ! outputs + real(wp), dimension(ncol,nlay,ngpt), intent(out) :: tau_rayleigh + !! Rayleigh optical depth + end subroutine compute_tau_rayleigh + end interface + ! ------------------------------------------------------------------------------------------------------------------ + interface + subroutine compute_Planck_source( & + ncol, nlay, nbnd, ngpt, & + nflav, neta, npres, ntemp, nPlanckTemp,& + tlay, tlev, tsfc, sfc_lay, & + fmajor, jeta, tropo, jtemp, jpress, & + gpoint_bands, band_lims_gpt, & + pfracin, temp_ref_min, totplnk_delta, totplnk, gpoint_flavor, & + sfc_src, lay_src, lev_src, sfc_source_Jac) bind(C, name="rrtmgp_compute_Planck_source") + use mo_rte_kind, only : wp, wl + integer, intent(in) :: ncol, nlay, nbnd, ngpt + !! input dimensions + integer, intent(in) :: nflav, neta, npres, ntemp, nPlanckTemp + !! table dimensions + real(wp), dimension(ncol,nlay ), intent(in) :: tlay !! temperature at layer centers (K) + real(wp), dimension(ncol,nlay+1), intent(in) :: tlev !! temperature at interfaces (K) + real(wp), dimension(ncol ), intent(in) :: tsfc !! surface temperture + integer, intent(in) :: sfc_lay !! index into surface layer + ! Interpolation variables + real(wp), dimension(2,2,2,ncol,nlay,nflav), intent(in) :: fmajor + !! interpolation weights for major gases - computed in interpolation() + integer, dimension(2, ncol,nlay,nflav), intent(in) :: jeta + !! interpolation indexes in eta - computed in interpolation() + logical(wl), dimension( ncol,nlay), intent(in) :: tropo + !! use upper- or lower-atmospheric tables? + integer, dimension( ncol,nlay), intent(in) :: jtemp, jpress + !! interpolation indexes in temperature and pressure - computed in interpolation() + ! Table-specific + integer, dimension(ngpt), intent(in) :: gpoint_bands !! band to which each g-point belongs + integer, dimension(2, nbnd), intent(in) :: band_lims_gpt !! start and end g-point for each band + real(wp), dimension(ntemp,neta,npres+1,ngpt), intent(in) :: pfracin !! Fraction of the Planck function in each g-point + real(wp), intent(in) :: temp_ref_min, totplnk_delta !! interpolation constants + real(wp), dimension(nPlanckTemp,nbnd), intent(in) :: totplnk !! Total Planck function by band at each temperature + integer, dimension(2,ngpt), intent(in) :: gpoint_flavor !! major gas flavor (pair) by upper/lower, g-point + + real(wp), dimension(ncol, ngpt), intent(out) :: sfc_src !! Planck emssion from the surface + real(wp), dimension(ncol,nlay, ngpt), intent(out) :: lay_src !! Planck emssion from layer centers + real(wp), dimension(ncol,nlay+1,ngpt), intent(out) :: lev_src !! Planck emission at layer boundaries + real(wp), dimension(ncol, ngpt), intent(out) :: sfc_source_Jac + !! Jacobian (derivative) of the surface Planck source with respect to surface temperature + end subroutine compute_Planck_source + end interface + ! ------------------------------------------------------------------------------------------------------------------ +end module mo_gas_optics_rrtmgp_kernels diff --git a/rrtmgp-kernels/api/rrtmgp_kernels.h b/rrtmgp-kernels/api/rrtmgp_kernels.h new file mode 100644 index 000000000..92c36957c --- /dev/null +++ b/rrtmgp-kernels/api/rrtmgp_kernels.h @@ -0,0 +1,123 @@ +/* This code is part of RRTMGP + +Contacts: Robert Pincus and Eli Mlawer +email: rrtmgp@aer.com + +Copyright 2024- + Trustees of Columbia University in the City of New York + 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 + +This header files defines the C bindings for the kernels used in RRTMGP + Adapted from code written by Chiel van Heerwaarden at Wageningen University and Research + +*/ +#pragma once +#include "rte_types.h" + +extern "C" +{ + void rrtmgp_interpolation( + const int& ncol, const int& nlay, + const int& ngas, const int& nflav, const int& neta, + const int& npres, const int& ntemp, + const int* flavor, // (2,nflav) + const Float* press_ref_log, // (npres) + const Float* temp_ref, // (ntemp) + const Float& press_ref_log_delta, + const Float& temp_ref_min, + const Float& temp_ref_delta, + const Float& press_ref_trop_log, + const Float* vmr_ref, //(2,ngas+1,ntemp) + const Float* play, // (ncol,nlay) + const Float* tlay, // (ncol,nlay) + const Float* col_gas, // (ncol,nlay,ngas+1) + int* jtemp, // [out] (ncol*nlay) + Float* fmajor, // [out] (2,2,2,ncol,nlay,nflav) + Float* fminor, // [out[ (2,2, ncol,nlay,nflav)) + Float* col_mix, // [out] (2, ncol,nlay,nflav) + Bool* tropo, // [out] size (ncol*nlay) + int* jeta,// [out] size (2*ncol*nlay*nflav) + int* jpress // [out] size (ncol*nlay) + ); + + void rrtmgp_compute_tau_absorption( + 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& nminorlower, const int& nminorklower, + const int& nminorupper, const int& nminorkupper, + const int& idx_h2o, + const int* gpoint_flavor, // (2,ngpt) + const int* band_lims_gpt, // (2,nbnd) + const Float* kmajor, // (ntemp,neta,npres+1,ngpt) + const Float* kminor_lower, // (ntemp,neta,nminorklower) + const Float* kminor_upper, // (ntemp,neta,nminorkupper) + const int* minor_limits_gpt_lower, // (2,nminorlower) + const int* minor_limits_gpt_upper, // (2,nminorupper) + const Bool* minor_scales_with_density_lower, // ( nminorlower) + const Bool* minor_scales_with_density_upper,// ( nminorupper) + const Bool* scale_by_complement_lower,// ( nminorlower) + const Bool* scale_by_complement_upper,// ( nminorupper) + const int* idx_minor_lower, // ( nminorlower) + const int* idx_minor_upper, // ( nminorupper) + const int* idx_minor_scaling_lower,// ( nminorlower) + const int* idx_minor_scaling_upper,// ( nminorupper) + const int* kminor_start_lower, // ( nminorlower) + const int* kminor_start_upper,// ( nminorupper) + const Bool* tropo, // (ncol,nlay) + const Float* col_mix, // (2, ncol,nlay,nflav ) + const Float* fmajor, // (2,2,2,ncol,nlay,nflav ) + const Float* fminor, // (2,2, ncol,nlay,nflav ) + const Float* play, // (ncol,nlay) + const Float* tlay, // (ncol,nlay) + const Float* col_gas, // (ncol,nlay,ngas+1) + const int* jeta, // (2, ncol,nlay,nflav ) + const int* jtemp, // (ncol,nlay) + const int* jpress, // (ncol,nlay) + Float* tau // [inout] (ncol,nlay.ngpt) + ); + + 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_Planck_source( + const int& ncol, const int& nlay, const int& nbnd, const int& ngpt, + const int& nflav, const int& neta, const int& npres, const int& ntemp, + const int& nPlanckTemp, + const Float* tlay, // (ncol,nlay ) + const Float* tlev, // (ncol,nlay+1) + const Float* tsfc, //(ncol ) + const int& sfc_lay, + const Float* fmajor, // (2,2,2,ncol,nlay,nflav) + const int* jeta, // (2, ncol,nlay,nflav) + const Bool* tropo, // ( ncol,nlay) + const int* jtemp, // ( ncol,nlay) + const int* jpress, // ( ncol,nlay) + const int* gpoint_bands, // (ngpt) + const int* band_lims_gpt, // (2, nbnd) + const Float* pfracin, // (ntemp,neta,npres+1,ngpt) + const Float& temp_ref_min, const Float& totplnk_delta, + const Float* totplnk, // (nPlanckTemp,nbnd) + const int* gpoint_flavor, // (2,ngpt) + Float* sfc_src, // [out] (ncol, ngpt) + Float* lay_src, // [out] (ncol,nlay, ngpt) + Float* lev_src, // [out] (ncol,nlay+1,ngpt) + Float* sfc_src_jac // [out] (ncol, ngpt) + ); +} diff --git a/rte-kernels/api/mo_fluxes_broadband_kernels.F90 b/rte-kernels/api/mo_fluxes_broadband_kernels.F90 new file mode 100644 index 000000000..d64c1fa29 --- /dev/null +++ b/rte-kernels/api/mo_fluxes_broadband_kernels.F90 @@ -0,0 +1,74 @@ +! This code is part of Radiative Transfer for Energetics (RTE) +! +! Contacts: Robert Pincus and Eli Mlawer +! email: rrtmgp@aer.com +! +! Copyright 2015-, Atmospheric and Environmental Research, +! Regents of the University of Colorado, 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 +! ------------------------------------------------------------------------------------------------- +!> +!> ## Kernels for computing broadband fluxes +!> +! ------------------------------------------------------------------------------------------------- +module mo_fluxes_broadband_kernels + use, intrinsic :: iso_c_binding + use mo_rte_kind, only: wp + implicit none + private + public :: sum_broadband, net_broadband + + ! ---------------------------------------------------------------------------- + !> + !> Spectral reduction over all points + !> + interface + subroutine sum_broadband(ncol, nlev, ngpt, spectral_flux, broadband_flux) bind(C, name="rte_sum_broadband") + use mo_rte_kind, only: wp + integer, intent(in ) :: ncol, nlev, ngpt + !! Array sizes + real(wp), dimension(ncol, nlev, ngpt), intent(in ) :: spectral_flux + !! Spectrally-resolved flux + real(wp), dimension(ncol, nlev), intent(out) :: broadband_flux + !! Sum of spectrally-resolved flux over `ngpt` + end subroutine sum_broadband + end interface + ! ---------------------------------------------------------------------------- + !> + !> Spectral reduction over all points for net flux + !> Overloaded - which routine is called depends on arguments + !> + interface net_broadband + ! ---------------------------------------------------------------------------- + !> + !> Net flux from g-point fluxes up and down + !> + subroutine net_broadband_full(ncol, nlev, ngpt, spectral_flux_dn, spectral_flux_up, broadband_flux_net) & + bind(C, name="rte_net_broadband_full") + use mo_rte_kind, only: wp + integer, intent(in ) :: ncol, nlev, ngpt + !! Array sizes + real(wp), dimension(ncol, nlev, ngpt), intent(in ) :: spectral_flux_dn, spectral_flux_up + !! Spectrally-resolved flux up and down + real(wp), dimension(ncol, nlev), intent(out) :: broadband_flux_net + !! Net (down minus up) summed over `ngpt` + end subroutine net_broadband_full + ! ---------------------------------------------------------------------------- + !> + !> Net flux when bradband flux up and down are already available + !> + subroutine net_broadband_precalc(ncol, nlev, flux_dn, flux_up, broadband_flux_net) & + bind(C, name="rte_net_broadband_precalc") + use mo_rte_kind, only: wp + integer, intent(in ) :: ncol, nlev + !! Array sizes + real(wp), dimension(ncol, nlev), intent(in ) :: flux_dn, flux_up + !! Broadband downward and upward fluxes + real(wp), dimension(ncol, nlev), intent(out) :: broadband_flux_net + !! Net (down minus up) + end subroutine net_broadband_precalc + end interface net_broadband + ! ---------------------------------------------------------------------------- +end module mo_fluxes_broadband_kernels diff --git a/rte-kernels/api/mo_optical_props_kernels.F90 b/rte-kernels/api/mo_optical_props_kernels.F90 new file mode 100644 index 000000000..bb7e5116f --- /dev/null +++ b/rte-kernels/api/mo_optical_props_kernels.F90 @@ -0,0 +1,377 @@ +! This code is part of Radiative Transfer for Energetics (RTE) +! +! Contacts: Robert Pincus and Eli Mlawer +! email: rrtmgp@aer.com +! +! Copyright 2015-, Atmospheric and Environmental Research, +! Regents of the University of Colorado, 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 +! ------------------------------------------------------------------------------------------------- +! +!> ## Kernels for arrays of optical properties: +!> - delta-scaling +!> - adding two sets of properties +!> - extracting subsets along the column dimension +! +! ------------------------------------------------------------------------------------------------- + +module mo_optical_props_kernels + use, intrinsic :: iso_c_binding + use mo_rte_kind, only: wp, wl + implicit none + + public + + ! ------------------------------------------------------------------------------------------------- + ! + ! Delta-scaling is provided only for two-stream properties at present + ! + interface delta_scale_2str_kernel + ! ------------------------------------------------------------------------------------------------- + !> Delta-scale two-stream optical properties given user-provided value of \(f\) (forward scattering) + ! + pure subroutine delta_scale_2str_f_k(ncol, nlay, ngpt, tau, ssa, g, f) & + bind(C, name="rte_delta_scale_2str_f_k") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt + !! Array sizes + real(wp), dimension(ncol, nlay, ngpt), intent(inout) :: tau, ssa, g + !! Optical depth, single-scattering albedo, asymmetry parameter + real(wp), dimension(ncol, nlay, ngpt), intent(in ) :: f + !! User-provided forward-scattering fraction + end subroutine delta_scale_2str_f_k + ! --------------------------------- + !> Delta-scale assuming forward-scatternig fraction is the square of the asymmetry parameter + !> i.e. \(f = g^2\) + ! + pure subroutine delta_scale_2str_k(ncol, nlay, ngpt, tau, ssa, g) & + bind(C, name="rte_delta_scale_2str_k") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt + !! Array sizes + real(wp), dimension(ncol, nlay, ngpt), intent(inout) :: tau, ssa, g + !! Optical depth, single-scattering albedo, asymmetry parameter + end subroutine delta_scale_2str_k + end interface delta_scale_2str_kernel + ! ------------------------------------------------------------------------------------------------- + ! + ! Addition of optical properties: the first set are incremented by the second set. + ! + ! There are three possible representations of optical properties (scalar = optical depth only; + ! two-stream = tau, single-scattering albedo, and asymmetry factor g, and + ! n-stream = tau, ssa, and phase function moments p.) Thus we need nine routines, three for + ! each choice of representation on the left hand side times three representations of the + ! optical properties to be added. + ! + ! There are two sets of these nine routines. In the first the two sets of optical + ! properties are defined at the same spectral resolution. There is also a set of routines + ! to add properties defined at lower spectral resolution to a set defined at higher spectral + ! resolution (adding properties defined by band to those defined by g-point) + ! + ! ------------------------------------------------------------------------------------------------- + !> increase one absorption optical depth by a second value + interface + pure subroutine increment_1scalar_by_1scalar(ncol, nlay, ngpt, & + tau1, & + tau2) bind(C, name="rte_increment_1scalar_by_1scalar") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt !! array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1 !! optical properties to be modified + real(wp), dimension(ncol,nlay,ngpt), intent(in ) :: tau2 !! optical properties to be added to original + end subroutine increment_1scalar_by_1scalar + end interface + ! --------------------------------- + !> increase absorption optical depth with extinction optical depth (2-stream form) + interface + pure subroutine increment_1scalar_by_2stream(ncol, nlay, ngpt, & + tau1, & + tau2, ssa2) bind(C, name="rte_increment_1scalar_by_2stream") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt !! array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1 !! optical properties to be modified + real(wp), dimension(ncol,nlay,ngpt), intent(in ) :: tau2, ssa2 !! optical properties to be added to original + end subroutine increment_1scalar_by_2stream + end interface + ! --------------------------------- + !> increase absorption optical depth with extinction optical depth (n-stream form) + interface + pure subroutine increment_1scalar_by_nstream(ncol, nlay, ngpt, & + tau1, & + tau2, ssa2) bind(C, name="rte_increment_1scalar_by_nstream") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt !! array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1 !! optical properties to be modified + real(wp), dimension(ncol,nlay,ngpt), intent(in ) :: tau2, ssa2 !! optical properties to be added to original + end subroutine increment_1scalar_by_nstream + end interface + ! --------------------------------- + ! --------------------------------- + !> increment two-stream optical properties \(\tau, \omega_0, g\) with absorption optical depth + interface + pure subroutine increment_2stream_by_1scalar(ncol, nlay, ngpt, & + tau1, ssa1, & + tau2) bind(C, name="rte_increment_2stream_by_1scalar") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt !! array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1 !! optical properties to be modified + real(wp), dimension(ncol,nlay,ngpt), intent(in ) :: tau2 !! optical properties to be added to original + end subroutine increment_2stream_by_1scalar + end interface + ! --------------------------------- + !> increment two-stream optical properties \(\tau, \omega_0, g\) with a second set + interface + pure subroutine increment_2stream_by_2stream(ncol, nlay, ngpt, & + tau1, ssa1, g1, & + tau2, ssa2, g2) bind(C, name="rte_increment_2stream_by_2stream") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt !! array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1, g1 !! optical properties to be modified + real(wp), dimension(ncol,nlay,ngpt), intent(in ) :: tau2, ssa2, g2 !! optical properties to be added to original + end subroutine increment_2stream_by_2stream + end interface + ! --------------------------------- + !> increment two-stream optical properties \(\tau, \omega_0, g\) with _n_-stream + interface + pure subroutine increment_2stream_by_nstream(ncol, nlay, ngpt, nmom2, & + tau1, ssa1, g1, & + tau2, ssa2, p2) bind(C, name="rte_increment_2stream_by_nstream") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt, nmom2 !! array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1, g1 !! optical properties to be modified + real(wp), dimension(ncol,nlay,ngpt), intent(in ) :: tau2, ssa2 !! optical properties to be added to original + real(wp), dimension(nmom2, & + ncol,nlay,ngpt), intent(in ) :: p2 !! moments of the phase function to be added + end subroutine increment_2stream_by_nstream + end interface + ! --------------------------------- + ! --------------------------------- + !> increment _n_-stream optical properties \(\tau, \omega_0, p\) with absorption optical depth + interface + pure subroutine increment_nstream_by_1scalar(ncol, nlay, ngpt, & + tau1, ssa1, & + tau2) bind(C, name="rte_increment_nstream_by_1scalar") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt !! array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1 !! optical properties to be modified + real(wp), dimension(ncol,nlay,ngpt), intent(in ) :: tau2 !! optical properties to be added to original + end subroutine increment_nstream_by_1scalar + end interface + ! --------------------------------- + !> increment _n_-stream optical properties \(\tau, \omega_0, p\) with two-stream values + interface + pure subroutine increment_nstream_by_2stream(ncol, nlay, ngpt, nmom1, & + tau1, ssa1, p1, & + tau2, ssa2, g2) bind(C, name="rte_increment_nstream_by_2stream") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt, nmom1 !! array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1 !! optical properties to be modified + real(wp), dimension(nmom1, & + ncol,nlay,ngpt), intent(inout) :: p1 !! moments of the phase function be modified + real(wp), dimension(ncol,nlay,ngpt), intent(in ) :: tau2, ssa2, g2 !! optical properties to be added to original + end subroutine increment_nstream_by_2stream + end interface + ! --------------------------------- + !> increment one set of _n_-stream optical properties with another set + interface + pure subroutine increment_nstream_by_nstream(ncol, nlay, ngpt, nmom1, nmom2, & + tau1, ssa1, p1, & + tau2, ssa2, p2) bind(C, name="rte_increment_nstream_by_nstream") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt, nmom1, nmom2 !! array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1 !! optical properties to be modified + real(wp), dimension(nmom1, & + ncol,nlay,ngpt), intent(inout) :: p1 !! moments of the phase function be modified + real(wp), dimension(ncol,nlay,ngpt), intent(in ) :: tau2, ssa2 !! optical properties to be added to original + real(wp), dimension(nmom2, & + ncol,nlay,ngpt), intent(in ) :: p2 !! moments of the phase function to be added + end subroutine increment_nstream_by_nstream + end interface + ! ------------------------------------------------------------------------------------------------- + ! + ! Incrementing when the second set of optical properties is defined at lower spectral resolution + ! (e.g. by band instead of by gpoint) + ! + ! ------------------------------------------------------------------------------------------------- + !> increase one absorption optical depth defined on g-points by a second value defined on bands + interface + pure subroutine inc_1scalar_by_1scalar_bybnd(ncol, nlay, ngpt, & + tau1, & + tau2, & + nbnd, gpt_lims) bind(C, name="rte_inc_1scalar_by_1scalar_bybnd") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt, nbnd !! array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1 !! optical properties to be modified (defined on g-points) + real(wp), dimension(ncol,nlay,nbnd), intent(in ) :: tau2 !! optical properties to be added to original (defined on bands) + integer, dimension(2,nbnd), intent(in ) :: gpt_lims !! Starting and ending gpoint for each band + end subroutine inc_1scalar_by_1scalar_bybnd + end interface + ! --------------------------------- + !> increase absorption optical depth defined on g-points with extinction optical depth (2-stream form) defined on bands + interface + pure subroutine inc_1scalar_by_2stream_bybnd(ncol, nlay, ngpt, & + tau1, & + tau2, ssa2, & + nbnd, gpt_lims) bind(C, name="rte_inc_1scalar_by_2stream_bybnd") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt, nbnd !! array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1 !! optical properties to be modified (defined on g-points) + real(wp), dimension(ncol,nlay,nbnd), intent(in ) :: tau2, ssa2 !! optical properties to be added to original (defined on bands) + integer, dimension(2,nbnd), intent(in ) :: gpt_lims !! Starting and ending gpoint for each band + end subroutine inc_1scalar_by_2stream_bybnd + end interface + ! --------------------------------- + !> increase absorption optical depth defined on g-points with extinction optical depth (n-stream form) defined on bands + interface + pure subroutine inc_1scalar_by_nstream_bybnd(ncol, nlay, ngpt, & + tau1, & + tau2, ssa2, & + nbnd, gpt_lims) bind(C, name="rte_inc_1scalar_by_nstream_bybnd") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt, nbnd !! array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1 !! optical properties to be modified (defined on g-points) + real(wp), dimension(ncol,nlay,nbnd), intent(in ) :: tau2, ssa2 !! optical properties to be added to original (defined on bands) + integer, dimension(2,nbnd), intent(in ) :: gpt_lims !! Starting and ending gpoint for each band + end subroutine inc_1scalar_by_nstream_bybnd + end interface + ! --------------------------------- + !> increment two-stream optical properties \(\tau, \omega_0, g\) defined on g-points with absorption optical depth defined on bands + interface + pure subroutine inc_2stream_by_1scalar_bybnd(ncol, nlay, ngpt, & + tau1, ssa1, & + tau2, & + nbnd, gpt_lims) bind(C, name="rte_inc_2stream_by_1scalar_bybnd") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt, nbnd !! array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1 !! optical properties to be modified (defined on g-points) + real(wp), dimension(ncol,nlay,nbnd), intent(in ) :: tau2 !! optical properties to be added to original (defined on bands) + integer, dimension(2,nbnd), intent(in ) :: gpt_lims !! Starting and ending gpoint for each band + end subroutine inc_2stream_by_1scalar_bybnd + end interface + ! --------------------------------- + !> increment 2-stream optical properties defined on g-points with another set defined on bands + interface + pure subroutine inc_2stream_by_2stream_bybnd(ncol, nlay, ngpt, & + tau1, ssa1, g1, & + tau2, ssa2, g2, & + nbnd, gpt_lims) bind(C, name="rte_inc_2stream_by_2stream_bybnd") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt, nbnd !! array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1, g1 !! optical properties to be modified (defined on g-points) + real(wp), dimension(ncol,nlay,nbnd), intent(in ) :: tau2, ssa2, g2 !! optical properties to be added to original (defined on bands) + integer, dimension(2,nbnd), intent(in ) :: gpt_lims !! Starting and ending gpoint for each band + end subroutine inc_2stream_by_2stream_bybnd + end interface + ! --------------------------------- + !> increment 2-stream optical properties defined on g-points with _n_-stream properties set defined on bands + interface + pure subroutine inc_2stream_by_nstream_bybnd(ncol, nlay, ngpt, nmom2, & + tau1, ssa1, g1, & + tau2, ssa2, p2, & + nbnd, gpt_lims) bind(C, name="rte_inc_2stream_by_nstream_bybnd") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt, nmom2, nbnd !! array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1, g1 !! optical properties to be modified (defined on g-points) + real(wp), dimension(ncol,nlay,nbnd), intent(in ) :: tau2, ssa2 !! optical properties to be added to original (defined on bands) + real(wp), dimension(nmom2, & + ncol,nlay,nbnd), intent(in ) :: p2 !! moments of the phase function to be added + integer, dimension(2,nbnd), intent(in ) :: gpt_lims !! Starting and ending gpoint for each band + end subroutine inc_2stream_by_nstream_bybnd + end interface + ! --------------------------------- + ! --------------------------------- + !> increment _n_-stream optical properties defined on g-points with absorption optical depth defined on bands + interface + pure subroutine inc_nstream_by_1scalar_bybnd(ncol, nlay, ngpt, & + tau1, ssa1, & + tau2, & + nbnd, gpt_lims) bind(C, name="rte_inc_nstream_by_1scalar_bybnd") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt, nbnd !! array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1 !! optical properties to be modified (defined on g-points) + real(wp), dimension(ncol,nlay,nbnd), intent(in ) :: tau2 !! optical properties to be added to original (defined on bands) + integer, dimension(2,nbnd), intent(in ) :: gpt_lims !! Starting and ending gpoint for each band + end subroutine inc_nstream_by_1scalar_bybnd + end interface + ! --------------------------------- + !> increment n-stream optical properties defined on g-points with 2-stream properties set defined on bands + interface + pure subroutine inc_nstream_by_2stream_bybnd(ncol, nlay, ngpt, nmom1, & + tau1, ssa1, p1, & + tau2, ssa2, g2, & + nbnd, gpt_lims) bind(C, name="rte_inc_nstream_by_2stream_bybnd") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt, nmom1, nbnd !! array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1 !! optical properties to be modified (defined on g-points) + real(wp), dimension(nmom1, & + ncol,nlay,ngpt), intent(inout) :: p1 !! moments of the phase function be modified + real(wp), dimension(ncol,nlay,nbnd), intent(in ) :: tau2, ssa2, g2 !! optical properties to be added to original (defined on bands) + integer, dimension(2,nbnd), intent(in ) :: gpt_lims !! Starting and ending gpoint for each band + end subroutine inc_nstream_by_2stream_bybnd + end interface + ! --------------------------------- + !> increment _n_-stream optical properties defined on g-points with a second set defined on bands + interface + pure subroutine inc_nstream_by_nstream_bybnd(ncol, nlay, ngpt, nmom1, nmom2, & + tau1, ssa1, p1, & + tau2, ssa2, p2, & + nbnd, gpt_lims) bind(C, name="rte_inc_nstream_by_nstream_bybnd") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt, nmom1, nmom2, nbnd !! array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1 !! optical properties to be modified (defined on g-points) + real(wp), dimension(nmom1, & + ncol,nlay,ngpt), intent(inout) :: p1 !! moments of the phase function be modified + real(wp), dimension(ncol,nlay,nbnd), intent(in ) :: tau2, ssa2 !! optical properties to be added to original (defined on bands) + real(wp), dimension(nmom2, & + ncol,nlay,nbnd), intent(in ) :: p2 !! moments of the phase function to be added + integer, dimension(2,nbnd), intent(in ) :: gpt_lims !! Starting and ending gpoint for each band + end subroutine inc_nstream_by_nstream_bybnd + end interface + ! ------------------------------------------------------------------------------------------------- + ! + ! Subsetting, meaning extracting some portion of the 3D domain + ! + ! ------------------------------------------------------------------------------------------------- + !> + !> Extract a subset from the first dimension (normally columns) of a 3D field. + !> Applicable to most variables e.g. tau, ssa, g + !> + interface extract_subset + pure subroutine extract_subset_dim1_3d(ncol, nlay, ngpt, array_in, colS, colE, array_out) & + bind (C, name="rte_extract_subset_dim1_3d") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt !! Array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(in ) :: array_in !! Array to subset + integer, intent(in ) :: colS, colE !! Starting and ending index + real(wp), dimension(colE-colS+1,& + nlay,ngpt), intent(out) :: array_out !! subset of the input array + end subroutine extract_subset_dim1_3d + ! --------------------------------- + !> Extract a subset from the second dimension (normally columns) of a 4D field. + !> Applicable to phase function moments, where the first dimension is the moment + pure subroutine extract_subset_dim2_4d(nmom, ncol, nlay, ngpt, array_in, colS, colE, array_out) & + bind (C, name="rte_extract_subset_dim2_4d") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: nmom, ncol, nlay, ngpt !! Array sizes + real(wp), dimension(nmom,ncol,nlay,ngpt), intent(in ) :: array_in !! Array to subset + integer, intent(in ) :: colS, colE !! Starting and ending index + real(wp), dimension(nmom,colE-colS+1,& + nlay,ngpt), intent(out) :: array_out !! subset of the input array + end subroutine extract_subset_dim2_4d + ! --------------------------------- + ! + !> Extract the absorption optical thickness \(\tau_{abs} = 1 - \omega_0 \tau_{ext}\) + ! + pure subroutine extract_subset_absorption_tau(ncol, nlay, ngpt, tau_in, ssa_in, & + colS, colE, tau_out) & + bind (C, name="rte_extract_subset_absorption_tau") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt !! Array sizes + real(wp), dimension(ncol,nlay,ngpt), intent(in ) :: tau_in, ssa_in !! Optical thickness, single scattering albedo + integer, intent(in ) :: colS, colE !! Starting and ending index + real(wp), dimension(colE-colS+1,& + nlay,ngpt), intent(out) :: tau_out !! absorption optical thickness subset + end subroutine extract_subset_absorption_tau + end interface extract_subset +end module mo_optical_props_kernels diff --git a/rte-kernels/api/mo_rte_solver_kernels.F90 b/rte-kernels/api/mo_rte_solver_kernels.F90 new file mode 100644 index 000000000..da8d0f706 --- /dev/null +++ b/rte-kernels/api/mo_rte_solver_kernels.F90 @@ -0,0 +1,202 @@ +! This code is part of Radiative Transfer for Energetics (RTE) +! +! Contacts: Robert Pincus and Eli Mlawer +! email: rrtmgp@aer.com +! +! Copyright 2015-, Atmospheric and Environmental Research, +! Regents of the University of Colorado, 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 +! ------------------------------------------------------------------------------------------------- +! +!>## Numeric calculations for radiative transfer solvers +!> - Emission/absorption (no-scattering) calculations +!> - solver for multi-angle Gaussian quadrature +!> - Extinction-only calculation (direct solar beam) +!> - Two-stream calculations: +!> solvers for LW and SW with different boundary conditions and source functions +! +! ------------------------------------------------------------------------------------------------- +module mo_rte_solver_kernels + use, intrinsic :: iso_c_binding + use mo_rte_kind, only: wp, wl + implicit none + private + + public :: lw_solver_noscat, lw_solver_2stream, & + sw_solver_noscat, sw_solver_2stream + ! ------------------------------------------------------------------------------------------------- + ! + ! Top-level longwave kernels + ! + ! ------------------------------------------------------------------------------------------------- + ! + !> LW transport, no scattering, multi-angle quadrature + !> Users provide a set of weights and quadrature angles + !> Routine sums over single-angle solutions for each sets of angles/weights + ! + ! --------------------------------------------------------------- + interface + subroutine lw_solver_noscat(ncol, nlay, ngpt, top_at_1, & + nmus, Ds, weights, & + tau, & + lay_source, lev_source, & + sfc_emis, sfc_src, & + inc_flux, & + flux_up, flux_dn, & + do_broadband, broadband_up, broadband_dn, & + do_Jacobians, sfc_srcJac, flux_upJac, & + do_rescaling, ssa, g) bind(C, name="rte_lw_solver_noscat") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt + !! Number of columns, layers, g-points + logical(wl), intent(in ) :: top_at_1 + !! ilay = 1 is the top of the atmosphere? + integer, intent(in ) :: nmus + !! number of quadrature angles + real(wp), dimension (ncol, ngpt, & + nmus), intent(in ) :: Ds + !! quadrature secants + real(wp), dimension(nmus), intent(in ) :: weights + !! quadrature weights + real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: tau + !! Absorption optical thickness [] + real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: lay_source + !! Planck source at layer average temperature [W/m2] + real(wp), dimension(ncol,nlay+1,ngpt), intent(in ) :: lev_source + !! Planck source at layer edge for radiation [W/m2] + real(wp), dimension(ncol, ngpt), intent(in ) :: sfc_emis + !! Surface emissivity [] + real(wp), dimension(ncol, ngpt), intent(in ) :: sfc_src + !! Surface source function [W/m2] + real(wp), dimension(ncol, ngpt), intent(in ) :: inc_flux + !! Incident diffuse flux, probably 0 [W/m2] + real(wp), dimension(ncol,nlay+1,ngpt), target, & + intent( out) :: flux_up, flux_dn + !! Fluxes [W/m2] + ! + ! Optional variables - arrays aren't referenced if corresponding logical == False + ! + logical(wl), intent(in ) :: do_broadband + real(wp), dimension(ncol,nlay+1 ), target, & + intent( out) :: broadband_up, broadband_dn + !! Spectrally-integrated fluxes [W/m2] + logical(wl), intent(in ) :: do_Jacobians + !! compute Jacobian with respect to surface temeprature? + real(wp), dimension(ncol ,ngpt), intent(in ) :: sfc_srcJac + !! surface temperature Jacobian of surface source function [W/m2/K] + real(wp), dimension(ncol,nlay+1 ), target, & + intent( out) :: flux_upJac + !! surface temperature Jacobian of Radiances [W/m2-str / K] + logical(wl), intent(in ) :: do_rescaling + !! Approximate treatment of scattering (10.1175/JAS-D-18-0014.1) + real(wp), dimension(ncol,nlay ,ngpt), intent(in ) :: ssa, g + !! single-scattering albedo, asymmetry parameter + end subroutine lw_solver_noscat + end interface + ! ------------------------------------------------------------------------------------------------- + ! + !> Longwave two-stream calculation: + !> - combine RRTMGP-specific sources at levels + !> - compute layer reflectance, transmittance + !> - compute total source function at levels using linear-in-tau + !> - transport + ! + ! ------------------------------------------------------------------------------------------------- + interface + subroutine lw_solver_2stream (ncol, nlay, ngpt, top_at_1, & + tau, ssa, g, & + lay_source, lev_source, sfc_emis, sfc_src, & + inc_flux, & + flux_up, flux_dn) bind(C, name="rte_lw_solver_2stream") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt + !! Number of columns, layers, g-points + logical(wl), intent(in ) :: top_at_1 + !! ilay = 1 is the top of the atmosphere? + real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: tau, ssa, g + !! Optical thickness, single-scattering albedo, asymmetry parameter [] + real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: lay_source + !! Planck source at layer average temperature [W/m2] + real(wp), dimension(ncol,nlay+1,ngpt), intent(in ) :: lev_source + !! Planck source at layer edge for radiation [W/m2] + real(wp), dimension(ncol, ngpt), intent(in ) :: sfc_emis + !! Surface emissivity [] + real(wp), dimension(ncol, ngpt), intent(in ) :: sfc_src + !! Surface source function [W/m2] + real(wp), dimension(ncol, ngpt), intent(in ) :: inc_flux + !! Incident diffuse flux, probably 0 [W/m2] + real(wp), dimension(ncol,nlay+1,ngpt), intent( out) :: flux_up, flux_dn + !! Fluxes [W/m2] + end subroutine lw_solver_2stream + end interface + ! ------------------------------------------------------------------------------------------------- + ! + ! Top-level shortwave kernels + ! + ! ------------------------------------------------------------------------------------------------- + ! + ! !> Extinction-only shortwave solver i.e. solar direct beam + ! + ! ------------------------------------------------------------------------------------------------- + interface + pure subroutine sw_solver_noscat(ncol, nlay, ngpt, top_at_1, & + tau, mu0, inc_flux_dir, flux_dir) bind(C, name="rte_sw_solver_noscat") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt ! Number of columns, layers, g-points + !! Number of columns, layers, g-points + logical(wl), intent(in ) :: top_at_1 + !! ilay = 1 is the top of the atmosphere? + real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: tau + !! Absorption optical thickness [] + real(wp), dimension(ncol,nlay ), intent(in ) :: mu0 + !! cosine of solar zenith angle + real(wp), dimension(ncol, ngpt), intent(in ) :: inc_flux_dir + !! Direct beam incident flux [W/m2] + real(wp), dimension(ncol,nlay+1,ngpt), intent(out) :: flux_dir + end subroutine sw_solver_noscat + end interface + ! ------------------------------------------------------------------------------------------------- + ! + !> Shortwave two-stream calculation: + !> compute layer reflectance, transmittance + !> compute solar source function for diffuse radiation + !> transport + ! + ! ------------------------------------------------------------------------------------------------- + interface + subroutine sw_solver_2stream (ncol, nlay, ngpt, top_at_1, & + tau, ssa, g, mu0, & + sfc_alb_dir, sfc_alb_dif, & + inc_flux_dir, & + flux_up, flux_dn, flux_dir, & + has_dif_bc, inc_flux_dif, & + do_broadband, broadband_up, & + broadband_dn, broadband_dir) bind(C, name="rte_sw_solver_2stream") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ncol, nlay, ngpt + !! Number of columns, layers, g-points + logical(wl), intent(in ) :: top_at_1 + !! ilay = 1 is the top of the atmosphere? + real(wp), dimension(ncol,nlay, ngpt), intent(in ) :: tau, ssa, g + !! Optical thickness, single-scattering albedo, asymmetry parameter [] + real(wp), dimension(ncol,nlay ), intent(in ) :: mu0 + !! cosine of solar zenith angle + real(wp), dimension(ncol, ngpt), intent(in ) :: sfc_alb_dir, sfc_alb_dif + !! Spectral surface albedo for direct and diffuse radiation + real(wp), dimension(ncol, ngpt), intent(in ) :: inc_flux_dir + !! Direct beam incident flux + real(wp), dimension(ncol,nlay+1,ngpt), target, & + intent( out) :: flux_up, flux_dn, flux_dir + !! Fluxes [W/m2] + logical(wl), intent(in ) :: has_dif_bc + !! Is a boundary condition for diffuse flux supplied? + real(wp), dimension(ncol, ngpt), intent(in ) :: inc_flux_dif + !! Boundary condition for diffuse flux [W/m2] + logical(wl), intent(in ) :: do_broadband + !! Provide broadband-integrated, not spectrally-resolved, fluxes? + real(wp), dimension(ncol,nlay+1 ), intent( out) :: broadband_up, broadband_dn, broadband_dir + end subroutine sw_solver_2stream + end interface +end module mo_rte_solver_kernels diff --git a/rte-kernels/api/mo_rte_util_array.F90 b/rte-kernels/api/mo_rte_util_array.F90 new file mode 100644 index 000000000..cdae473c4 --- /dev/null +++ b/rte-kernels/api/mo_rte_util_array.F90 @@ -0,0 +1,47 @@ +! This code is part of Radiative Transfer for Energetics (RTE) +! +! Contacts: Robert Pincus and Eli Mlawer +! email: rrtmgp@aer.com +! +! Copyright 2015- Atmospheric and Environmental Research, +! Regents of the University of Colorado, +! Trustees of Columbia University in the City of New York +! 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_rte_util_array + use mo_rte_kind, only: wp, wl + implicit none + public :: zero_array + + !------------------------------------------------------------------------------------------------- + ! Initializing arrays to 0 + !------------------------------------------------------------------------------------------------- + interface zero_array + subroutine zero_array_1D(ni, array) bind(C, name="zero_array_1D") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ni + real(wp), dimension(ni), intent(out) :: array + end subroutine zero_array_1D + ! ---------------------------------------------------------- + subroutine zero_array_2D(ni, nj, array) bind(C, name="zero_array_2D") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ni, nj + real(wp), dimension(ni, nj), intent(out) :: array + end subroutine zero_array_2D + ! ---------------------------------------------------------- + subroutine zero_array_3D(ni, nj, nk, array) bind(C, name="zero_array_3D") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ni, nj, nk + real(wp), dimension(ni, nj, nk), intent(out) :: array + end subroutine zero_array_3D + ! ---------------------------------------------------------- + subroutine zero_array_4D(ni, nj, nk, nl, array) bind(C, name="zero_array_4D") + use mo_rte_kind, only: wp, wl + integer, intent(in ) :: ni, nj, nk, nl + real(wp), dimension(ni, nj, nk, nl), intent(out) :: array + end subroutine zero_array_4D + end interface zero_array +end module mo_rte_util_array diff --git a/rte-kernels/api/rte_kernels.h b/rte-kernels/api/rte_kernels.h new file mode 100644 index 000000000..71f86b81e --- /dev/null +++ b/rte-kernels/api/rte_kernels.h @@ -0,0 +1,389 @@ +/* This code is part of Radiative Transfer for Energetics (RTE) + +Contacts: Robert Pincus and Eli Mlawer +email: rrtmgp@aer.com + +Copyright 2024- + Trustees of Columbia University in the City of New York + 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 + +This header files defines the C bindings for the kernels used in RTE + Adapted from code written by Chiel van Heerwaarden at Wageningen University and Research + +*/ + +#include "rte_types.h" + +extern "C" +{ + // + // Shortwave solvers + // + void rte_sw_solver_noscat( + const int& ncol, const int& nlay, const int& ngpt, const Bool& top_at_1, + const Float* tau, // (ncol,nlay, ngpt) + const Float* mu0, // (ncol,nlay) + const Float* inc_flux_dir, // (ncol, ngpt) + Float* flux_dir); // [out] (ncol,nlay+1,ngpt) + + void rte_sw_solver_2stream( + const int& ncol, + const int& nlay, + const int& ngpt, + const Bool& top_at_1, + const Float* tau, // (ncol,nlay, ngpt) + const Float* ssa, // (ncol,nlay, ngpt) + const Float* g, // (ncol,nlay, ngpt) + const Float* mu0, // (ncol,nlay) + const Float* sfc_alb_dir, // (ncol, ngpt) + const Float* sfc_alb_dif, // (ncol, ngpt) + const Float* inc_flux_dir, // (ncol, ngpt) + Float* flux_up, // [out] (ncol,nlay+1,ngpt) + Float* flux_dn, // [out] (ncol,nlay+1,ngpt) + Float* flux_dir, // [out] (ncol,nlay+1,ngpt) + const Bool& has_dif_bc, + const Float* inc_flux_dif, // (ncol, ngpt) + const Bool& do_broadband, + Float* broadband_up, // [out] (ncol,nlay+1) + Float* broadband_dn, // [out] (ncol,nlay+1) + Float* broadband_dir); // [out] (ncol,nlay+1) + + void rte_lw_solver_noscat( + const int& ncol, + const int& nlay, + const int& ngpt, + const Bool& top_at_1, + const int& nmus, + const Float* secants, // (nmus) + const Float* weights, // (nmus) + const Float* tau, // (ncol,nlay, ngpt) + const Float* lay_source, // (ncol,nlay, ngpt) + const Float* lev_source, // (ncol,nlay+1,ngpt) + const Float* sfc_emis, // (ncol, ngpt) + const Float* sfc_src, // (ncol, ngpt) + const Float* inc_flux, // (ncol, ngpt) + Float* flux_up, // [out] (ncol,nlay+1,ngpt) + Float* flux_dn, // [out] (ncol,nlay+1,ngpt) + const Bool& do_broadband, + Float* broadband_up, + // [out] (ncol,nlay+1) + Float* broadband_dn, + // [out] (ncol,nlay+1) + const Bool& do_jacobians, + const Float* sfc_src_jac, + // (ncol, ngpt) + Float* flux_up_jac, + // [out] (ncol,nlay+1,ngpt) + const Bool& do_rescaling, + const Float* ssa, // (ncol,nlay, ngpt) + const Float* g); // (ncol,nlay, ngpt) + + + void rte_lw_solver_2stream( + const int& ncol, + const int& nlay, + const int& ngpt, + const Bool& top_at_1, + const Float* tau, // (ncol,nlay, ngpt) + const Float* ssa, // (ncol,nlay, ngpt) + const Float* g, // (ncol,nlay, ngpt) + const Float* lay_source, // (ncol,nlay, ngpt) + const Float* lev_source, // (ncol,nlay+1,ngpt) + const Float* sfc_emis, // (ncol, ngpt) + const Float* sfc_src, // (ncol, ngpt) + const Float* inc_flux, // (ncol, ngpt) + Float* flux_up, // [out] (ncol,nlay+1,ngpt) + Float* flux_dn // [out] (ncol,nlay+1,ngpt) + ); + + // + // OPTICAL PROPS - INCREMENT + // + void rte_increment_1scalar_by_1scalar( + const int& ncol, + const int& nlay, + const int& ngpt, + Float* tau_inout, // [inout] (ncol,nlay,ngpt) + const Float* tau_in); // (ncol,nlay,ngpt) + + + void rte_increment_1scalar_by_2stream( + const int& ncol, + const int& nlay, + const int& ngpt, + Float* tau_inout, // [inout] (ncol,nlay,ngpt) + const Float* tau_in, // (ncol,nlay,ngpt) + const Float* ssa_in); // (ncol,nlay,ngpt) + + void rte_increment_1scalar_by_nstream( + const int& ncol, + const int& nlay, + const int& ngpt, + Float* tau_inout, // [inout] (ncol,nlay,ngpt) + const Float* tau_in, // (ncol,nlay,ngpt) + const Float* ssa_in); // (ncol,nlay,ngpt) + + void rte_increment_2stream_by_1scalar( + const int& ncol, + const int& nlay, + const int& ngpt, + Float* tau_inout, // [inout] (ncol,nlay,ngpt) + Float* ssa_inout, // [inout] (ncol,nlay,ngpt) + const Float* tau_in); // (ncol,nlay,ngpt) + + void rte_increment_2stream_by_2stream( + const int& ncol, + const int& nlay, + const int& ngpt, + Float* tau_inout, // [inout] (ncol,nlay,ngpt) + Float* ssa_inout, // [inout] (ncol,nlay,ngpt) + Float* g_inout, // [inout] (ncol,nlay,ngpt) + const Float* tau_in, // (ncol,nlay,ngpt) + const Float* ssa_in, // (ncol,nlay,ngpt) + const Float* g_in); // (ncol,nlay,ngpt) + + + void rte_increment_2stream_by_nstream( + const int& ncol, const int& nlay, const int& ngpt, const int& nmom, + Float* tau_inout, // [inout] (ncol,nlay,ngpt) + Float* ssa_inout, // [inout] (ncol,nlay,ngpt) + Float* g_inout, // [inout] (ncol,nlay,ngpt) + const Float* tau_in, // (ncol,nlay,ngpt) + const Float* ssa_in, // (ncol,nlay,ngpt) + const Float* p_in); //(nmom,ncol,nlay,ngpt) + + void rte_increment_nstream_by_1scalar( + const int& ncol, + const int& nlay, + const int& ngpt, + Float* tau_inout, // [inout] (ncol,nlay,ngpt) + const Float* tau_in, // (ncol,nlay,ngpt) + const Float* ssa_in); // (ncol,nlay,ngpt) + + void rte_increment_nstream_by_2stream( + const int& ncol, + const int& nlay, + const int& ngpt, + const int& nmom1, + Float* tau_inout, // [inout] (ncol,nlay,ngpt) + Float* ssa_inout, // [inout] (ncol,nlay,ngpt) + Float* p_inout, // [inout] (nmom,ncol,nlay,ngpt) + const Float* tau_in, // (ncol,nlay,ngpt) + const Float* ssa_in, // (ncol,nlay,ngpt) + const Float* g_in); // (ncol,nlay,ngpt) + + void rte_increment_nstream_by_nstream( + const int& ncol, + const int& nlay, + const int& ngpt, + const int& nmom1, + const int& nmom2, + Float* tau_inout, // [inout] (ncol,nlay,ngpt) + Float* ssa_inout, // [inout] (ncol,nlay,ngpt) + Float* p_inout, // [inout](nmom1,ncol,nlay,ngpt) + const Float* tau_in, // (ncol,nlay,ngpt) + const Float* ssa_in, // (ncol,nlay,ngpt) + const Float* p_in); // (nmom2,ncol,nlay,ngpt) + + // + // OPTICAL PROPS - INCREMENT BYBND + // + void rte_inc_1scalar_by_1scalar_bybnd( + const int& ncol, + const int& nlay, + const int& ngpt, + Float* tau_inout,// [inout] (ncol,nlay,ngpt) + const Float* tau_in, // (ncol,nlay,nbnd) + const int& nbnd, + const int* band_lims_gpoint); // (2,nbnd) + + void rte_inc_1scalar_by_2stream_bybnd( + const int& ncol, + const int& nlay, + const int& ngpt, + Float* tau_inout, // [inout] (ncol,nlay,ngpt) + const Float* tau_in, // (ncol,nlay,nbnd) + const Float* ssa_in, // (ncol,nlay,nbnd) + const int& nbnd, + const int* band_lims_gpoint); // (2,nbnd) + + void rte_inc_1scalar_by_nstream_bybnd( + const int& ncol, + const int& nlay, + const int& ngpt, + Float* tau_inout, // [inout] (ncol,nlay,ngpt) + const Float* tau_in, // (ncol,nlay,nbnd) + const Float* ssa_in, // (ncol,nlay,nbnd) + const int& nbnd, + const int* band_lims_gpoint); // (2,nbnd) + + void rte_inc_2stream_by_1scalar_bybnd( + const int& ncol, + const int& nlay, + const int& ngpt, + Float* tau_inout, // [inout] (ncol,nlay,ngpt) + Float* ssa_inout, // [inout] (ncol,nlay,ngpt) + const Float* tau_in, // (ncol,nlay,nbnd) + const int& nbnd, + const int* band_lims_gpoint); // (2,nbnd) + + void rte_inc_2stream_by_2stream_bybnd( + const int& ncol, + const int& nlay, + const int& ngpt, + Float* tau_inout, // [inout] (ncol,nlay,ngpt) + Float* ssa_inout, // [inout] (ncol,nlay,ngpt) + Float* g_inout, // [inout] (ncol,nlay,ngpt) + const Float* tau_in, // (ncol,nlay,nbnd) + const Float* ssa_in, // (ncol,nlay,nbnd) + const Float* g_in, // (ncol,nlay,nbnd) + const int& nbnd, + const int* band_lims_gpoint); // (2,nbnd) + + void rte_inc_2stream_by_nstream_bybnd( + const int& ncol, + const int& nlay, + const int& ngpt, + const int& nmom, + Float* tau_inout, // [inout] (ncol,nlay,ngpt) + Float* ssa_inout, // [inout] (ncol,nlay,ngpt) + Float* g_inout, // [inout] (ncol,nlay,ngpt) + const Float* tau_in, // (ncol,nlay,nbnd) + const Float* ssa_in, // (ncol,nlay,nbnd) + const Float* p_in, // (nmom,ncol,nlay,nbnd) + const int& nbnd, + const int* band_lims_gpoint); // (2,nbnd) + + void rte_inc_nstream_by_1scalar_bybnd( + const int& ncol, + const int& nlay, + const int& ngpt, + Float* tau_inout, // [inout] (ncol,nlay,ngpt) + Float* ssa_inout, // [inout] (ncol,nlay,ngpt) + const Float* tau_in, // (ncol,nlay,nbnd) + const int& nbnd, + const int* band_lims_gpoint); // (2,nbnd) + + void rte_inc_nstream_by_2stream_bybnd( + const int& ncol, + const int& nlay, + const int& ngpt, + const int& nmom1, + Float* tau_inout, // [inout] (ncol,nlay,ngpt) + Float* ssa_inout, // [inout] (ncol,nlay,ngpt) + Float* p_inout, + // [inout] (nomo,ncol,nlay,ngpt) + const Float* tau_in, // (ncol,nlay,nbnd) + const Float* ssa_in, // (ncol,nlay,nbnd) + const Float* g_in, // (ncol,nlay,nbnd) + const int& nbnd, + const int* band_lims_gpoint); // (2,nbnd) + + void rte_inc_nstream_by_nstream_bybnd( + const int& ncol, + const int& nlay, + const int& ngpt, + const int& nmom1, + const int& nmom2, + Float* tau_inout, // [inout] (ncol,nlay,ngpt) + Float* ssa_inout, // [inout] (ncol,nlay,ngpt) + Float* p_inout, + // [inout] (nomo,ncol,nlay,ngpt) + const Float* tau_in, // (ncol,nlay,nbnd) + const Float* ssa_in, // (ncol,nlay,nbnd) + const Float* p_in, // (nmom,ncol,nlay,nbnd) + const int& nbnd, + const int* band_lims_gpoint); // (2,nbnd) + + // + // OPTICAL PROPS - DELTA SCALING + // + void rte_delta_scale_2str_k( + const int& ncol, const int& nlay, const int& ngpt, + Float* tau_inout, // [inout] (ncol,nlev,ngpt) + Float* ssa_inout, // [inout] (ncol,nlev,ngpt) + Float* g_inout); // [inout] (ncol,nlev,ngpt) + + void rte_delta_scale_2str_f_k( + const int& ncol, const int& nlay, const int& ngpt, + Float* tau_inout, // [inout] (ncol,nlev,ngpt) + Float* ssa_inout, // [inout] (ncol,nlev,ngpt) + Float* g_inout, // [inout] (ncol,nlev,ngpt) + const Float* f); // (ncol,nlev,ngpt) + + // + // OPTICAL PROPS - SUBSET + // + void rte_extract_subset_dim1_3d( + const int& ncol, const int& nlay, const int& ngpt, + Float* array_in, // (ncol,nlay,ngpt) + const int& ncol_start, const int& ncol_end, + Float* array_out); // [out] (ncol_end-ncol_start+1,nlay,ngpt) + + void rte_extract_subset_dim2_4d( + const int& nmom, const int& ncol, const int& nlay, const int& ngpt, + const Float* array_in, // (nmom,ncol,nlay,ngpt) + const int& ncol_start, const int& ncol_end, + Float* array_out); // [out] (nmom,ncol_end-ncol_start+1,nlay,ngpt) + + void rte_extract_subset_absorption_tau( + const int& ncol, const int& nlay, const int& ngpt, + const Float* tau_in, // (ncol,nlay,ngpt) + const Float* ssa_in, // (ncol,nlay,ngpt) + const int& ncol_start, const int& ncol_end, + Float* tau_out); // [out] (ncol_end-ncol_start+1,nlay,ngpt) + + // + // Fluxes - reduction + // + void rte_sum_broadband( + const int& ncol, const int& nlev, const int& ngpt, + const Float* gpt_flux, // (ncol,nlev,ngpt) + Float* flux); // [out] (ncol,nlev) + + void rte_net_broadband_full( + const int& ncol, const int& nlev, const int& ngpt, + const Float* gpt_flux_dn, // (ncol,nlev,ngpt) + const Float* gpt_flux_up, // (ncol,nlev,ngpt) + Float* flux_net); // [out] (ncol,nlev) + + void rte_net_broadband_precalc( + const int& ncol, const int& nlev, + const Float* broadband_flux_dn, // (ncol, nlev) + const Float* broadband_flux_up, // (ncol, nlev) + Float* broadband_flux_net);//[out] (ncol, nlev) + + void rte_sum_byband( + const int& ncol, const int& nlev, const int& ngpt, const int& nbnd, + const int* band_lims, // dimension(2, nbnd) + const Float* gpt_flux, // (ncol, nlev, ngpt) + Float* bnd_flux); // [out] (ncol, nlev, ngpt) + + void rte_net_byband_full( + const int& ncol, const int& nlev, const int& ngpt, const int& nbnd, int* band_lims, + const Float* bnd_flux_dn, // (ncol,nlev,nbnd) + const Float* bnd_flux_up, // (ncol,nlev,nbnd) + Float* bnd_flux_net); // [out] (ncol,nlev) + // + // Array utilities + // + void zero_array_1D( + const int& ni, + Float* array); // [out] (ni) + + void zero_array_2D( + const int& ni, const int& nj, + Float* array); // [out] (ni, nj) + + void zero_array_3D( + const int& ni, const int& nj, const int& nk, + Float* array); // [out] (ni, nj, nk) + + void zero_array_4D( + const int& ni, const int& nj, const int& nk, const int& nl, + Float* array); // [out] (ni, nj, nk, nl) + +} diff --git a/rte-kernels/api/rte_types.h b/rte-kernels/api/rte_types.h new file mode 100644 index 000000000..fc645f203 --- /dev/null +++ b/rte-kernels/api/rte_types.h @@ -0,0 +1,34 @@ +/* This code is part of Radiative Transfer for Energetics (RTE) + +Contacts: Robert Pincus and Eli Mlawer +email: rrtmgp@aer.com + +Copyright 2024- + Trustees of Columbia University in the City of New York + 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 + +This header files C-compatible Boolean and floating point types (see mo_rte_type.F90 for the Fortran equivalent) + Adapted from code written by Chiel van Heerwaarden at Wageningen University and Research + +*/ +#pragma once + +#ifdef RTE_USE_CBOOL +using Bool = signed char; +#else +using Bool = int; +#endif + +#ifdef RTE_USE_SP +using Float = float; +const Float Float_epsilon = FLT_EPSILON; +#else +using Float = double; +const Float Float_epsilon = DBL_EPSILON; +#endif + + + diff --git a/tests/Makefile b/tests/Makefile index 1b90ed26c..a5060c623 100644 --- a/tests/Makefile +++ b/tests/Makefile @@ -1,3 +1,4 @@ +#!/usr/bin/env make # # Location of RTE+RRTMGP libraries, module files. # @@ -5,20 +6,18 @@ RRTMGP_BUILD = $(RRTMGP_ROOT)/build # # RRTMGP library, module files # -# LDFLAGS += -L$(RRTMGP_BUILD) -# LIBS += -lrrtmgp -lrte +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 -LDFLAGS += -L$(RRTMGP_BUILD) -LIBS += -lrte -lrrtmgp LIBS += -lnetcdff -lnetcdf -VPATH = .:$(RRTMGP_ROOT)/examples:$(RRTMGP_ROOT)/examples/rfmip-clear-sky:$(RRTMGP_ROOT)/examples/all-sky +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 @@ -80,7 +79,7 @@ rte_sw_solver_unit_tests : $(LIB_DEPS) mo_testing_utils.o rte_sw_solver_unit_te .PHONY: tests -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 From bbc14c20a2d9c6c1375128a24ba9978789cdb525 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Thu, 18 Apr 2024 07:32:42 -0400 Subject: [PATCH 17/20] Re-vectorize SW two-stream (#275) Adopt suggestions by @peterukk in #215. Makes expensive direct beam computations even when solar zenith angle is < 0, but masks those results out in layers where the direct beam does not reach. This will have practical impacts only in columns where the top of the column is in sunlight but the bottom not (or if users aren't removing columns entirely below the horizon). --- .github/workflows/containerized-ci.yml | 2 +- rte-kernels/mo_rte_solver_kernels.F90 | 113 +++++++++++++------------ tests/check_equivalence.F90 | 64 ++++++++------ 3 files changed, 99 insertions(+), 80 deletions(-) diff --git a/.github/workflows/containerized-ci.yml b/.github/workflows/containerized-ci.yml index 0f247ba05..2548e2940 100644 --- a/.github/workflows/containerized-ci.yml +++ b/.github/workflows/containerized-ci.yml @@ -21,7 +21,7 @@ jobs: 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 + fcflags: -m64 -g -traceback -heap-arrays -assume realloc_lhs -extend-source 132 -check bounds,uninit,pointers,stack -stand f08 -diag-disable=10448 # Set flags for Intel Fortran Compiler - fortran-compiler: ifx rte-kernels: default diff --git a/rte-kernels/mo_rte_solver_kernels.F90 b/rte-kernels/mo_rte_solver_kernels.F90 index f3308cc5e..e10a0cb36 100644 --- a/rte-kernels/mo_rte_solver_kernels.F90 +++ b/rte-kernels/mo_rte_solver_kernels.F90 @@ -1003,6 +1003,7 @@ pure subroutine sw_dif_and_source(ncol, nlay, top_at_1, mu0, sfc_albedo, & ! Ancillary variables real(wp), parameter :: min_k = 1.e4_wp * epsilon(1._wp) ! Suggestion from Chiel van Heerwaarden + real(wp), parameter :: min_mu0 = sqrt(epsilon(1._wp)) real(wp) :: k, exp_minusktau, k_mu, k_gamma3, k_gamma4 real(wp) :: RT_term, exp_minus2ktau real(wp) :: Rdir, Tdir, Tnoscat @@ -1022,6 +1023,7 @@ pure subroutine sw_dif_and_source(ncol, nlay, top_at_1, mu0, sfc_albedo, & dir_flux_trans => flux_dn_dir(:,lay_index ) end if + !$OMP SIMD do i = 1, ncol ! ! Scalars @@ -1029,7 +1031,6 @@ pure subroutine sw_dif_and_source(ncol, nlay, top_at_1, mu0, sfc_albedo, & tau_s = tau(i, lay_index) w0_s = w0 (i, lay_index) g_s = g (i, lay_index) - mu0_s = mu0(i, lay_index) ! ! Zdunkowski Practical Improved Flux Method "PIFM" ! (Zdunkowski et al., 1980; Contributions to Atmospheric Physics 53, 147-66) @@ -1059,63 +1060,69 @@ pure subroutine sw_dif_and_source(ncol, nlay, top_at_1, mu0, sfc_albedo, & ! ! On a round earth, where mu0 can increase with depth in the atmosphere, ! levels with mu0 <= 0 have no direct beam and hence no source for diffuse light + ! Compute transmission and reflection using a nominal value but mask out later ! - if(mu0_s > 0._wp) then - k_mu = k * mu0_s - ! - ! Equation 14, multiplying top and bottom by exp(-k*tau) - ! and rearranging to avoid div by 0. - ! - RT_term = w0_s * RT_term/merge(1._wp - k_mu*k_mu, & - epsilon(1._wp), & - abs(1._wp - k_mu*k_mu) >= epsilon(1._wp)) - ! - ! Zdunkowski Practical Improved Flux Method "PIFM" - ! (Zdunkowski et al., 1980; Contributions to Atmospheric Physics 53, 147-66) - ! - gamma3 = (2._wp - 3._wp * mu0_s * g_s ) * .25_wp - gamma4 = 1._wp - gamma3 - alpha1 = gamma1 * gamma4 + gamma2 * gamma3 ! Eq. 16 - alpha2 = gamma1 * gamma3 + gamma2 * gamma4 ! Eq. 17 + mu0_s = max(min_mu0, mu0(i, lay_index)) + k_mu = k * mu0_s + ! + ! Equation 14, multiplying top and bottom by exp(-k*tau) + ! and rearranging to avoid div by 0. + ! + RT_term = w0_s * RT_term/merge(1._wp - k_mu*k_mu, & + epsilon(1._wp), & + abs(1._wp - k_mu*k_mu) >= epsilon(1._wp)) + ! + ! Zdunkowski Practical Improved Flux Method "PIFM" + ! (Zdunkowski et al., 1980; Contributions to Atmospheric Physics 53, 147-66) + ! + gamma3 = (2._wp - 3._wp * mu0_s * g_s ) * .25_wp + gamma4 = 1._wp - gamma3 + alpha1 = gamma1 * gamma4 + gamma2 * gamma3 ! Eq. 16 + alpha2 = gamma1 * gamma3 + gamma2 * gamma4 ! Eq. 17 - ! - ! Transmittance of direct, unscattered beam. - ! - k_gamma3 = k * gamma3 - k_gamma4 = k * gamma4 - Tnoscat = exp(-tau_s/mu0_s) - Rdir = RT_term * & - ((1._wp - k_mu) * (alpha2 + k_gamma3) - & - (1._wp + k_mu) * (alpha2 - k_gamma3) * exp_minus2ktau - & - 2.0_wp * (k_gamma3 - alpha2 * k_mu) * exp_minusktau * Tnoscat) - ! - ! Equation 15, multiplying top and bottom by exp(-k*tau), - ! multiplying through by exp(-tau/mu0) to - ! prefer underflow to overflow - ! Omitting direct transmittance - ! - Tdir = -RT_term * & - ((1._wp + k_mu) * (alpha1 + k_gamma4) * Tnoscat - & - (1._wp - k_mu) * (alpha1 - k_gamma4) * exp_minus2ktau * Tnoscat - & - 2.0_wp * (k_gamma4 + alpha1 * k_mu) * exp_minusktau) - ! Final check that energy is not spuriously created, by recognizing that - ! the beam can either be reflected, penetrate unscattered to the base of a layer, - ! or penetrate through but be scattered on the way - the rest is absorbed - ! Makes the equations safer in single precision. Credit: Robin Hogan, Peter Ukkonen - Rdir = max(0.0_wp, min(Rdir, (1.0_wp - Tnoscat ) )) - Tdir = max(0.0_wp, min(Tdir, (1.0_wp - Tnoscat - Rdir) )) + ! + ! Transmittance of direct, unscattered beam. + ! + k_gamma3 = k * gamma3 + k_gamma4 = k * gamma4 + Tnoscat = exp(-tau_s/mu0_s) + Rdir = RT_term * & + ((1._wp - k_mu) * (alpha2 + k_gamma3) - & + (1._wp + k_mu) * (alpha2 - k_gamma3) * exp_minus2ktau - & + 2.0_wp * (k_gamma3 - alpha2 * k_mu) * exp_minusktau * Tnoscat) + ! + ! Equation 15, multiplying top and bottom by exp(-k*tau), + ! multiplying through by exp(-tau/mu0) to + ! prefer underflow to overflow + ! Omitting direct transmittance + ! + Tdir = -RT_term * & + ((1._wp + k_mu) * (alpha1 + k_gamma4) * Tnoscat - & + (1._wp - k_mu) * (alpha1 - k_gamma4) * exp_minus2ktau * Tnoscat - & + 2.0_wp * (k_gamma4 + alpha1 * k_mu) * exp_minusktau) + ! Final check that energy is not spuriously created, by recognizing that + ! the beam can either be reflected, penetrate unscattered to the base of a layer, + ! or penetrate through but be scattered on the way - the rest is absorbed + ! Makes the equations safer in single precision. Credit: Robin Hogan, Peter Ukkonen + Rdir = max(0.0_wp, min(Rdir, (1.0_wp - Tnoscat ) )) + Tdir = max(0.0_wp, min(Tdir, (1.0_wp - Tnoscat - Rdir) )) - source_up(i,lay_index) = Rdir * dir_flux_inc(i) - source_dn(i,lay_index) = Tdir * dir_flux_inc(i) - dir_flux_trans(i) = Tnoscat * dir_flux_inc(i) - else - source_up(i,lay_index) = 0._wp - source_dn(i,lay_index) = 0._wp - dir_flux_trans(i) = 0._wp - end if + source_up(i,lay_index) = Rdir * dir_flux_inc(i) + source_dn(i,lay_index) = Tdir * dir_flux_inc(i) + dir_flux_trans(i) = Tnoscat * dir_flux_inc(i) end do end do - source_sfc(:) = dir_flux_trans(:)*sfc_albedo(:) + ! + ! T and R for the direct beam are computed using nominal values even when the + ! sun is below the horizon (mu0 < 0); set those values back to zero + ! This won't be efficient if many nighttime columns are passed + ! + source_sfc(:) = merge(dir_flux_trans(:)*sfc_albedo(:), & + 0._wp, mu0(:,lay_index) > 0._wp) + where(mu0(:,:) <= 0._wp) + source_up(:,:) = 0._wp + source_dn(:,:) = 0._wp + end where end subroutine sw_dif_and_source ! --------------------------------------------------------------- diff --git a/tests/check_equivalence.F90 b/tests/check_equivalence.F90 index ed6022788..f08ca8a26 100644 --- a/tests/check_equivalence.F90 +++ b/tests/check_equivalence.F90 @@ -421,55 +421,67 @@ program rte_check_equivalence ! Threshold of 4x spacing() works in double precision ! call stop_on_err(gas_optics%gas_optics(p_lay, p_lev, & - t_lay, & - gas_concs, & - atmos, & - toa_flux)) + t_lay, & + gas_concs, & + atmos, & + 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, & 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 = 6._wp) .or. & - .not. allclose(tst_flux_dir,ref_flux_dir,tol = 8._wp)) & + 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)) & call report_err(" halving/doubling fails") + ! + ! Incremement with 0 optical depth + ! + call stop_on_err(gas_optics%gas_optics(p_lay, p_lev, & + t_lay, & + gas_concs, & + atmos, & + toa_flux)) call increment_with_1scl(atmos) call stop_on_err(rte_sw(atmos, top_at_1, & 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 = 6._wp) .or. & - .not. allclose(tst_flux_dir,ref_flux_dir,tol = 6._wp)) & + 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)) & call report_err(" Incrementing with 1scl fails") - call stop_on_err(gas_optics%gas_optics(p_lay, p_lev, & - t_lay, & - gas_concs, & - atmos, & - toa_flux)) - call increment_with_2str(atmos) - if(.not. allclose(tst_flux_up, ref_flux_up, tol = 8._wp) .or. & - .not. allclose(tst_flux_dn, ref_flux_dn, tol = 6._wp) .or. & - .not. allclose(tst_flux_dir,ref_flux_dir,tol = 6._wp)) & + call stop_on_err(gas_optics%gas_optics(p_lay, p_lev, & + t_lay, & + gas_concs, & + atmos, & + toa_flux)) + call increment_with_2str(atmos) + call stop_on_err(rte_sw(atmos, top_at_1, & + 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)) & call report_err(" Incrementing with 2str fails") call stop_on_err(gas_optics%gas_optics(p_lay, p_lev, & - t_lay, & - gas_concs, & - atmos, & - toa_flux)) + t_lay, & + gas_concs, & + atmos, & + toa_flux)) call increment_with_nstr(atmos) call stop_on_err(rte_sw(atmos, top_at_1, & 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 = 6._wp) .or. & - .not. allclose(tst_flux_dir,ref_flux_dir,tol = 6._wp)) & + 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)) & call report_err(" Incrementing with nstr fails") print *, " Incrementing" end if From 8cf510ce4d50e9ddb56dd30449a9590067194e58 Mon Sep 17 00:00:00 2001 From: Sergey Kosukhin Date: Thu, 18 Apr 2024 15:05:55 +0200 Subject: [PATCH 18/20] Skip GitLab CI for dependabot (#278) --- .github/workflows/gitlab-ci.yml | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/.github/workflows/gitlab-ci.yml b/.github/workflows/gitlab-ci.yml index 0c987d844..374def61a 100644 --- a/.github/workflows/gitlab-ci.yml +++ b/.github/workflows/gitlab-ci.yml @@ -21,7 +21,8 @@ jobs: if: | 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.head.repo.owner.login == github.repository_owner && + github.event.pull_request.user.login != 'dependabot' )) runs-on: ubuntu-latest outputs: ref-name: ${{ steps.g-push-rev.outputs.ref-name }} @@ -103,7 +104,8 @@ jobs: if: | 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.head.repo.owner.login == github.repository_owner && + github.event.pull_request.user.login != 'dependabot' )) runs-on: ubuntu-latest outputs: ref-name: ${{ steps.g-push-rev.outputs.ref-name }} From 29cc4b37384684e663cb3ad32cf9a3dca9e39cc1 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Thu, 18 Apr 2024 09:13:32 -0400 Subject: [PATCH 19/20] Different user for dependabot jobs --- .github/workflows/gitlab-ci.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/gitlab-ci.yml b/.github/workflows/gitlab-ci.yml index 374def61a..2a46e8557 100644 --- a/.github/workflows/gitlab-ci.yml +++ b/.github/workflows/gitlab-ci.yml @@ -22,7 +22,7 @@ jobs: 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' )) + github.event.pull_request.user.login != 'dependabot[bot]' )) runs-on: ubuntu-latest outputs: ref-name: ${{ steps.g-push-rev.outputs.ref-name }} @@ -105,7 +105,7 @@ jobs: 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' )) + github.event.pull_request.user.login != 'dependabot[bot]' )) runs-on: ubuntu-latest outputs: ref-name: ${{ steps.g-push-rev.outputs.ref-name }} From 9caac5ec3ceab9a295699a2abb58809a9dca87de Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Thu, 18 Apr 2024 09:24:04 -0400 Subject: [PATCH 20/20] Bump JamesIves/github-pages-deploy-action from 4.5.0 to 4.6.0 (#276) Bumps [JamesIves/github-pages-deploy-action](https://github.com/jamesives/github-pages-deploy-action) from 4.5.0 to 4.6.0. - [Release notes](https://github.com/jamesives/github-pages-deploy-action/releases) - [Commits](https://github.com/jamesives/github-pages-deploy-action/compare/v4.5.0...v4.6.0) --- updated-dependencies: - dependency-name: JamesIves/github-pages-deploy-action dependency-type: direct:production update-type: version-update:semver-minor ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/doc-deployment.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/doc-deployment.yml b/.github/workflows/doc-deployment.yml index 2ee1de23b..d948d3453 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.5.0 + uses: JamesIves/github-pages-deploy-action@v4.6.0 if: ${{ github.event_name == 'push' && github.ref == 'refs/heads/documentation' }} with: branch: gh-pages