From db7044f9c85b9675b0447d2162001f32b61e0e84 Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Mon, 13 Jan 2025 11:44:04 -0500 Subject: [PATCH 1/3] Reformatting only: remove tailing whitespace, normalize file endings, cma ke-format and cmake-lint, (Github) actionslint via pre-commit --- .github/workflows/continuous-integration.yml | 4 +- .github/workflows/gitlab-ci.yml | 2 +- .github/workflows/module_switcher | 1 - CITATION.cff | 12 +- Contributing.md | 12 +- README.md | 16 +- doc/class_diagrams/mo_optical_props.puml | 4 +- doc/jekyll_site/Gemfile | 1 - .../_posts/2022-06-02-Release-notes.md | 2 +- .../_posts/2023-11-27-v1.7-Release-notes.md | 12 +- doc/jekyll_site/explanations/index.md | 2 +- doc/jekyll_site/how-tos/build-and-test.md | 50 ++-- doc/jekyll_site/index.markdown | 10 +- doc/jekyll_site/reference/index.md | 2 +- doc/jekyll_site/tutorials/index.md | 2 +- examples/all-sky/README.md | 2 +- examples/all-sky/make_problem_size_loop.py | 2 +- .../all-sky/mo_load_aerosol_coefficients.F90 | 8 +- .../all-sky/mo_load_cloud_coefficients.F90 | 2 +- examples/all-sky/rrtmgp_allsky.F90 | 268 +++++++++--------- examples/compare-to-reference.py | 16 +- examples/rfmip-clear-sky/CMakeLists.txt | 3 +- examples/rfmip-clear-sky/Readme.md | 2 +- examples/rfmip-clear-sky/mo_rfmip_io.F90 | 2 +- extensions/mo_cloud_sampling.F90 | 2 +- extensions/mo_compute_bc.F90 | 2 +- extensions/mo_fluxes_byband.F90 | 4 +- extensions/mo_fluxes_bygpoint.F90 | 2 +- extensions/mo_heating_rates.F90 | 4 +- .../mo_zenith_angle_spherical_correction.F90 | 2 +- gas-optics/mo_gas_concentrations.F90 | 52 ++-- gas-optics/mo_gas_optics.F90 | 38 +-- gas-optics/mo_gas_optics_util_string.F90 | 4 +- .../mo_aerosol_optics_rrtmgp_merra.F90 | 50 ++-- rrtmgp-frontend/mo_cloud_optics_rrtmgp.F90 | 4 +- rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 | 48 ++-- .../accel/mo_gas_optics_rrtmgp_kernels.F90 | 6 +- .../api/mo_cloud_optics_rrtmgp_kernels.F90 | 3 +- .../api/mo_gas_optics_rrtmgp_kernels.F90 | 116 ++++---- rrtmgp-kernels/api/rrtmgp_kernels.h | 20 +- .../mo_cloud_optics_rrtmgp_kernels.F90 | 2 +- .../mo_gas_optics_rrtmgp_kernels.F90 | 122 ++++---- rte-frontend/mo_optical_props.F90 | 16 +- rte-frontend/mo_rte_lw.F90 | 14 +- rte-frontend/mo_rte_sw.F90 | 2 +- rte-frontend/mo_rte_util_array_validation.F90 | 2 +- .../accel/mo_optical_props_kernels.F90 | 38 +-- rte-kernels/accel/mo_rte_solver_kernels.F90 | 12 +- .../api/mo_fluxes_broadband_kernels.F90 | 4 +- rte-kernels/api/mo_optical_props_kernels.F90 | 36 +-- rte-kernels/api/rte_kernels.h | 88 +++--- rte-kernels/api/rte_types.h | 3 - rte-kernels/mo_rte_solver_kernels.F90 | 32 +-- tests/Readme.md | 2 +- tests/check_equivalence.F90 | 134 ++++----- tests/check_variants.F90 | 22 +- tests/mo_gas_optics_defs_rrtmgp.F90 | 6 +- tests/mo_rcemip_profiles.F90 | 2 +- tests/mo_testing_utils.F90 | 132 ++++----- tests/rte_lw_solver_unit_tests.F90 | 124 ++++---- tests/rte_optic_prop_unit_tests.F90 | 146 +++++----- tests/rte_sw_solver_unit_tests.F90 | 112 ++++---- tests/validation-plots.py | 4 +- 63 files changed, 922 insertions(+), 927 deletions(-) diff --git a/.github/workflows/continuous-integration.yml b/.github/workflows/continuous-integration.yml index d9e346e9e..baff11923 100644 --- a/.github/workflows/continuous-integration.yml +++ b/.github/workflows/continuous-integration.yml @@ -24,7 +24,7 @@ jobs: gfortran-version: [12, 13, 14] gfortran-from: [system, conda] fpmodel: [DP, SP] - exclude: + exclude: - os: ubuntu-24.04 gfortran-from: conda - os: macos-13 @@ -38,7 +38,7 @@ jobs: gfortran-version: 14 gfortran-from: conda - os: windows-2022 - include: + include: - os: ubuntu-24.04 gfortran-version: 13 gfortran-from: conda diff --git a/.github/workflows/gitlab-ci.yml b/.github/workflows/gitlab-ci.yml index feacb5f7f..4b08fe737 100644 --- a/.github/workflows/gitlab-ci.yml +++ b/.github/workflows/gitlab-ci.yml @@ -103,7 +103,7 @@ jobs: # lumi-init: if: | - false && + false && github.repository_owner == 'earth-system-radiation' && ( github.event_name != 'pull_request' || ( github.event.pull_request.head.repo.owner.login == github.repository_owner && diff --git a/.github/workflows/module_switcher b/.github/workflows/module_switcher index 3037308c0..caad97d57 100644 --- a/.github/workflows/module_switcher +++ b/.github/workflows/module_switcher @@ -143,4 +143,3 @@ switch_for_module () eval "$sfm_cmd" done } - diff --git a/CITATION.cff b/CITATION.cff index b574ec618..dfa642b9b 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -5,10 +5,10 @@ authors: given-names: Robert orcid: "https://orcid.org/0000-0002-0016-3470" - family-names: Iacono - given-names: Michael J. + given-names: Michael J. orcid: "https://orcid.org/0000-0002-9578-0649" - family-names: Alexeev - given-names: Dmitry + given-names: Dmitry orcid: "https://orcid.org/0000-0002-6425-2181" - family-names: Adamidis given-names: Panos @@ -19,12 +19,12 @@ authors: given-names: Matthew orcid: "https://orcid.org/0000-0003-4764-3348" - family-names: Pfister - given-names: Erik + given-names: Erik orcid: "http://orcid.org/0009-0002-2688-138X" - family-names: Polonsky given-names: Igor N. - family-names: Romero - given-names: Nicols A. + given-names: Nicols A. - family-names: Kosukhin given-names: Sergey S. - family-names: Wehe @@ -35,7 +35,7 @@ license: BSD-3-Clause date-released: "2023-11-27" version: 1.7 -abstract: "RTE+RRTMGP is a set of codes for computing radiative fluxes in planetary atmospheres. +abstract: "RTE+RRTMGP is a set of codes for computing radiative fluxes in planetary atmospheres. RRTMGP uses a k-distribution to provide an optical description (absorption and possibly Rayleigh optical depth) of the gaseous atmosphere, along with the relevant source functions, on a pre-determined spectral grid given temperatures, pressures, and gas concentration. RTE computes fluxes given spectrally-resolved optical descriptions and source functions. The fluxes are normally summarized or reduced via a user extensible class." identifiers: - description: "All versions" @@ -65,4 +65,4 @@ references: volume: 11 number: 10 pages: 3074-3089 - year: 2019 \ No newline at end of file + year: 2019 diff --git a/Contributing.md b/Contributing.md index 0531ed768..bbaa1bebc 100644 --- a/Contributing.md +++ b/Contributing.md @@ -4,16 +4,16 @@ Thanks for considering making a contribution to RTE+RRTMGP. The code in this repository is intended to work with compilers supporting the Fortran 2008 standard. It is also expected to run end-to-end on GPUs when compiled with OpenACC or OpenMP (though OpenMP is still unreliable). Commits are tested automatically against a range of compilers using Github Actions and also resources provided by the [Swiss Supercomputing Center](https://cscs.ch). The testing uses two general codes in `examples/`for which results are compared against existing implemetations, and custom codes in `tests/` intended to excercise all code options. -##### Did you find a bug? +##### Did you find a bug? -Please file an issue on the [Github page](https://github.com/RobertPincus/rte-rrtmgp/issues). Please include a minimal reproducer of the bug it at all possible. +Please file an issue on the [Github page](https://github.com/RobertPincus/rte-rrtmgp/issues). Please include a minimal reproducer of the bug it at all possible. ##### Did you write a patch that fixes a bug? -Please fork this repository, branch from `develop`, make your changes, and open a Github [pull request](https://github.com/RobertPincus/rte-rrtmgp/pulls) against branch `develop`. +Please fork this repository, branch from `develop`, make your changes, and open a Github [pull request](https://github.com/RobertPincus/rte-rrtmgp/pulls) against branch `develop`. -##### Did you add functionality? +##### Did you add functionality? -Please fork this repository, branch from `develop`, make your changes, and open a Github [pull request](https://github.com/RobertPincus/rte-rrtmgp/pulls) against branch `develop`, adding a new regression test or comparison against the reference in `tests/verification.py` or `tests/validation-plots.py` as appropriate. Add the test to the `tests` target in `tests/Makefile`. +Please fork this repository, branch from `develop`, make your changes, and open a Github [pull request](https://github.com/RobertPincus/rte-rrtmgp/pulls) against branch `develop`, adding a new regression test or comparison against the reference in `tests/verification.py` or `tests/validation-plots.py` as appropriate. Add the test to the `tests` target in `tests/Makefile`. -RTE+RRTMGP is intended to be a core that users can extend with custom code to suit their own needs. \ No newline at end of file +RTE+RRTMGP is intended to be a core that users can extend with custom code to suit their own needs. diff --git a/README.md b/README.md index dfbb39610..123374a51 100644 --- a/README.md +++ b/README.md @@ -3,7 +3,7 @@ [RTE+RRTMGP's GitHub Pages site](https://earth-system-radiation.github.io/rte-rrtmgp/) contains a mix of automatically-generated documentation and hand-written descriptions. The documentation is incomplete and evolving. Thanks to the folks at [Sourcery Institute](https://www.sourceryinstitute.org) -for help in setting this up. +for help in setting this up. For the moment the [Wiki](https://github.com/earth-system-radiation/rte-rrtmgp/wiki) may also be useful. @@ -20,22 +20,22 @@ RTE computes fluxes given spectrally-resolved optical descriptions and source fu A description of building RTE+RRTMGP with an ad hoc homemade system is described in the [documentation](https://earth-system-radiation.github.io/rte-rrtmgp/how-tos/). -See also the `autoconf` branch for a Gnu autotools build system. +See also the `autoconf` branch for a Gnu autotools build system. ## Examples Two examples are provided in `examples/`, one for clear skies and one including clouds. Directory `tests/` contains regression testing (e.g. to ensure that answers are independent of orientation) and unit testing (to be sure all the code paths are tested). See the README file and codes in each directory for further information. -## Citing the code +## Citing the code -Code releases are archived at Zenodo. All releases are available at -[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.3403172.svg)](https://doi.org/10.5281/zenodo.3403172). +Code releases are archived at Zenodo. All releases are available at +[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.3403172.svg)](https://doi.org/10.5281/zenodo.3403172). The current release is available at: [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.7521518.svg)](https://doi.org/10.5281/zenodo.10211873) Please cite the code using these DOIs and the information in the `CITATION.cff` file in addition to the reference [paper](https://doi.org/10.1029/2019MS001621) ## Acknowledgements -The development of RTE+RRTMGP has been funded in the US by the Office of Naval Research, NASA, NOAA, and the Department of Energy. We -are grateful for contributions from a range of collaborators at institutions including the Swiss Supercomputing Center, -the German Climate Computing Center, and Nvidia. \ No newline at end of file +The development of RTE+RRTMGP has been funded in the US by the Office of Naval Research, NASA, NOAA, and the Department of Energy. We +are grateful for contributions from a range of collaborators at institutions including the Swiss Supercomputing Center, +the German Climate Computing Center, and Nvidia. diff --git a/doc/class_diagrams/mo_optical_props.puml b/doc/class_diagrams/mo_optical_props.puml index 00993d4b9..0e5bbe7ac 100644 --- a/doc/class_diagrams/mo_optical_props.puml +++ b/doc/class_diagrams/mo_optical_props.puml @@ -5,7 +5,7 @@ hide empty members class stuff_t{ z_ : complex - defined_ : logical + defined_ : logical z() : complex defined() : logical stuff_t(z : complex) : stuff_t @@ -26,4 +26,4 @@ class characterizable_stuff_t{ characterizable_stuff_t *-down- stuff_t : aggregates characterizable_stuff_t .up.|> characterizable_t : implements -@enduml \ No newline at end of file +@enduml diff --git a/doc/jekyll_site/Gemfile b/doc/jekyll_site/Gemfile index 743828e57..1a57fb1be 100644 --- a/doc/jekyll_site/Gemfile +++ b/doc/jekyll_site/Gemfile @@ -27,4 +27,3 @@ end # Performance-booster for watching directories on Windows gem "wdm", "~> 0.1.1", :platforms => [:mingw, :x64_mingw, :mswin] - diff --git a/doc/jekyll_site/_posts/2022-06-02-Release-notes.md b/doc/jekyll_site/_posts/2022-06-02-Release-notes.md index 83391b14d..9eabdd6db 100644 --- a/doc/jekyll_site/_posts/2022-06-02-Release-notes.md +++ b/doc/jekyll_site/_posts/2022-06-02-Release-notes.md @@ -4,7 +4,7 @@ title: "2022-06-02: Release notes" categories: Release-notes --- -Commit [a4fe30c](https://github.com/earth-system-radiation/rte-rrtmgp/commit/a4fe30cf4dab2e5fd8d3ab6f11683a82ae584475) +Commit [a4fe30c](https://github.com/earth-system-radiation/rte-rrtmgp/commit/a4fe30cf4dab2e5fd8d3ab6f11683a82ae584475) to branch `main` makes the following changes: - Solar zenith angle can vary with height, allowing for calculations on a pseudo-spherical earth diff --git a/doc/jekyll_site/_posts/2023-11-27-v1.7-Release-notes.md b/doc/jekyll_site/_posts/2023-11-27-v1.7-Release-notes.md index ef4cc27a1..e5d947c1d 100644 --- a/doc/jekyll_site/_posts/2023-11-27-v1.7-Release-notes.md +++ b/doc/jekyll_site/_posts/2023-11-27-v1.7-Release-notes.md @@ -4,12 +4,12 @@ title: "v1.7 Release notes" categories: Release-notes --- -Commit [3ac0636](https://github.com/earth-system-radiation/rte-rrtmgp/commit/3ac0636b17d6a3f11e4085f91679393fceaa4e18) +Commit [3ac0636](https://github.com/earth-system-radiation/rte-rrtmgp/commit/3ac0636b17d6a3f11e4085f91679393fceaa4e18) to branch `main` makes the following changes: -- Libraries can be built in single precision by changes in `rte-kind/mo_rte_kind.F90`. Differences with respect to double precision are roughly 0.13 W/m2. -- A class for computing the optical properties of aerosols following the MERRA representation has been added. -- The repository is reorganized into `frontend` and `kernel` directories for `rte` and `rrtmgp`. Data has been moved to a separate [repository]((https://github.com/earth-system-radiation/rrtmgp-data/). -- Citation information has been added. +- Libraries can be built in single precision by changes in `rte-kind/mo_rte_kind.F90`. Differences with respect to double precision are roughly 0.13 W/m2. +- A class for computing the optical properties of aerosols following the MERRA representation has been added. +- The repository is reorganized into `frontend` and `kernel` directories for `rte` and `rrtmgp`. Data has been moved to a separate [repository]((https://github.com/earth-system-radiation/rrtmgp-data/). +- Citation information has been added. -As usual some bugs have been fixed, the use of OpenACC and OpenMP GPU offload directives continues to evolve, and the continous integration continues to be fine-tuned. +As usual some bugs have been fixed, the use of OpenACC and OpenMP GPU offload directives continues to evolve, and the continous integration continues to be fine-tuned. diff --git a/doc/jekyll_site/explanations/index.md b/doc/jekyll_site/explanations/index.md index 24ea3e65a..7b205be25 100644 --- a/doc/jekyll_site/explanations/index.md +++ b/doc/jekyll_site/explanations/index.md @@ -9,4 +9,4 @@ title: "Explanations" The spectral properties of the atmosphere and the source functions depend on electromagnetic wavelength (or frequency or wavenumber). RTE treats this spectral dependence by dividing the spectrum into one or more _bands_, each of which represents a continuous set of wavelengths/frequencies. Bands may be further sub-divided into _g-points_ (the language is borrowed from _k_-distributions). Each _g_-point is treated as a independent psudo-monchromatic calculation but there is no inherent mapping between _g_-points and wavelengths; the sum over _g_-points is the band-average value. -The bands defined by RRTMGP cover the full spectrum of radiation emitted by the Sun and Earth: these are _broadband_ calculations. In RRTMGP the bands are continuous so that the ending wavelength of one band is the starting wavelength of the next. +The bands defined by RRTMGP cover the full spectrum of radiation emitted by the Sun and Earth: these are _broadband_ calculations. In RRTMGP the bands are continuous so that the ending wavelength of one band is the starting wavelength of the next. diff --git a/doc/jekyll_site/how-tos/build-and-test.md b/doc/jekyll_site/how-tos/build-and-test.md index bec95a35b..bf9417d08 100644 --- a/doc/jekyll_site/how-tos/build-and-test.md +++ b/doc/jekyll_site/how-tos/build-and-test.md @@ -4,12 +4,12 @@ title: "How to build and run tests" --- How to build the libraries, tests, and examples, run the tests, and verify the results -## In a nutshell -RTE+RRTMGP uses `CMake`. In the root directory: -`cmake -S . -B build` will guide you through configuration options. +## In a nutshell +RTE+RRTMGP uses `CMake`. In the root directory: +`cmake -S . -B build` will guide you through configuration options. -Environment variables can also be set and passed to `CMake` as in this example, which -builds and runs the tests: +Environment variables can also be set and passed to `CMake` as in this example, which +builds and runs the tests: ` cmake -S . -B build \ @@ -26,44 +26,44 @@ cmake -S . -B build \ Evaluating the results of the tests requires `Python` and the packages described in `environment*.yml`. -## Building and testing using (Gnu) autotools +## Building and testing using (Gnu) autotools Sergey Kosukhin and his colleagues at the Max Planck Institute for Meteorology maintain the `autoconf` branch which adds Gnu `autotools` building to `main` branch. -## Supplying data +## Supplying data -Running the tests and verifying the results requires the RRTMGP data. Clone the -[data repository](https://github.com/earth-system-radiation/rrtmgp-data) or download the -[Zenodo archive](https://doi.org/10.5281/zenodo.7988260). Set the environment variable `RRTMGP_DATA` -to the root of this directory. +Running the tests and verifying the results requires the RRTMGP data. Clone the +[data repository](https://github.com/earth-system-radiation/rrtmgp-data) or download the +[Zenodo archive](https://doi.org/10.5281/zenodo.7988260). Set the environment variable `RRTMGP_DATA` +to the root of this directory. -## Example compiler flags +## Example compiler flags In these examples `FC` is the Fortran compilers using flags `FCFLAGS` -### Gnu Fortran +### Gnu Fortran (see also the [continuous integration](https://github.com/earth-system-radiation/rte-rrtmgp/blob/main/.github/workflows/continuous-integration.yml)) `FC`: `gfortran-10` or `gfortran-11` or `gfortran-12` #### Debugging flags -`FCFLAGS: "-ffree-line-length-none -m64 -std=f2008 -march=native -fbounds-check -finit-real=nan -DRTE_USE_CBOOL"` +`FCFLAGS: "-ffree-line-length-none -m64 -std=f2008 -march=native -fbounds-check -finit-real=nan -DRTE_USE_CBOOL"` #### Even stricter debugging flags -`FCFLAGS: "-ffree-line-length-none -m64 -std=f2008 -march=native -fbounds-check -fbacktrace -finit-real=nan -DRTE_USE_CBOOL -pedantic -g -Wall"` +`FCFLAGS: "-ffree-line-length-none -m64 -std=f2008 -march=native -fbounds-check -fbacktrace -finit-real=nan -DRTE_USE_CBOOL -pedantic -g -Wall"` -### Intel Fortran Classic +### Intel Fortran Classic (see also the [continuous integration](https://github.com/earth-system-radiation/rte-rrtmgp/blob/main/.github/workflows/containerized-ci.yml)) -`FC: ifort` +`FC: ifort` #### Debugging flags -`FCFLAGS: "-m64 -g -traceback -heap-arrays -assume realloc_lhs -extend-source 132 -check bounds,uninit,pointers,stack -stand f08"` -#### Optimization flags: +`FCFLAGS: "-m64 -g -traceback -heap-arrays -assume realloc_lhs -extend-source 132 -check bounds,uninit,pointers,stack -stand f08"` +#### Optimization flags: `FCFLAGS:"-m64 -O3 -g -traceback -heap-arrays -assume realloc_lhs -extend-source 132"` -### Intel Fortran +### Intel Fortran (LLVM, see also the [continuous integration](https://github.com/earth-system-radiation/rte-rrtmgp/blob/main/.github/workflows/containerized-ci.yml)) -`FC: ifort` +`FC: ifort` #### Debugging flags -`FCFLAGS: "-debug -traceback -heap-arrays -assume realloc_lhs -extend-source 132 -stand f08"` -#### Using OpenMP GPU offload +`FCFLAGS: "-debug -traceback -heap-arrays -assume realloc_lhs -extend-source 132 -stand f08"` +#### Using OpenMP GPU offload See [this open issue](https://github.com/earth-system-radiation/rte-rrtmgp/issues/194) ### NVFortran @@ -71,10 +71,10 @@ See [this open issue](https://github.com/earth-system-radiation/rte-rrtmgp/issue `FC: nvfortran` #### Debugging flags `FCFLAGS: "-g -Minfo -Mbounds -Mchkptr -Mstandard -Kieee -Mchkstk -Mallocatable=03 -Mpreprocess"` -#### Optimization flags: +#### Optimization flags: `FCFLAGS: "-O3 -fast -Minfo -Mallocatable=03 -Mpreprocess"` ### HPE CCE for GPU using OpenMP-acc: crayftn -- requires at least CCE 14.0.0 `FC: crayftn` #### Debugging flags (these appear to be insufficient during the link stage) -`FCFLAGS: "-hnoacc -homp -O0"` \ No newline at end of file +`FCFLAGS: "-hnoacc -homp -O0"` diff --git a/doc/jekyll_site/index.markdown b/doc/jekyll_site/index.markdown index 90986395f..7fe51dd13 100644 --- a/doc/jekyll_site/index.markdown +++ b/doc/jekyll_site/index.markdown @@ -7,9 +7,9 @@ title: "RTE+RRTMGP documentation" This is the documentation for RTE+RRTMGP, a set of codes for computing radiative fluxes in planetary atmospheres. RTE+RRTMGP is described in a [paper](https://doi.org/10.1029/2019MS001621) in -[Journal of Advances in Modeling Earth Systems](http://james.agu.org). -The code itself can be sited as -doi:[10.5281/zenodo.3403172](https://doi.org/10.5281/zenodo.3403172) or via the +[Journal of Advances in Modeling Earth Systems](http://james.agu.org). +The code itself can be sited as +doi:[10.5281/zenodo.3403172](https://doi.org/10.5281/zenodo.3403172) or via the DOI attached to each release. RRTMGP uses a k-distribution to provide an optical description (absorption and @@ -34,5 +34,5 @@ We are starting with the [reference documentation](./reference/index.html), auto-generated from the code itself. This is provided separately for RTE and RRTMGP and for the user-facing classes and underlying computational kernels. -We welcome contributions to the documentation via pull requests to the `documentation` branch -of the Github repository. \ No newline at end of file +We welcome contributions to the documentation via pull requests to the `documentation` branch +of the Github repository. diff --git a/doc/jekyll_site/reference/index.md b/doc/jekyll_site/reference/index.md index 58dd27f49..a05bbe775 100644 --- a/doc/jekyll_site/reference/index.md +++ b/doc/jekyll_site/reference/index.md @@ -14,7 +14,7 @@ RTE aspires to follow a set of coding conventions: and spectral quadrature point. - RTE and RRTMGP are agnostic to vertical ordering - Units are MKS -- Procedures (with the exception of testing code) do not perform I/O +- Procedures (with the exception of testing code) do not perform I/O ## Fortran user-facing class interfaces diff --git a/doc/jekyll_site/tutorials/index.md b/doc/jekyll_site/tutorials/index.md index 87da31a76..1a9cd9d79 100644 --- a/doc/jekyll_site/tutorials/index.md +++ b/doc/jekyll_site/tutorials/index.md @@ -10,7 +10,7 @@ A typical workflow for a clear-sky calculation is to 1. allocate memory 2. set gas concentration values 3. compute the optical properties of the gaseous atmosphere -4. compute radiative fluxes +4. compute radiative fluxes This repository contains all the pieces needed to perform a clear-sky calculation. An [example](https://github.com/RobertPincus/rte-rrtmgp/tree/master/examples/rfmip-clear-sky) is provided. diff --git a/examples/all-sky/README.md b/examples/all-sky/README.md index 3a5e9a257..72b700510 100644 --- a/examples/all-sky/README.md +++ b/examples/all-sky/README.md @@ -4,4 +4,4 @@ This example provides a modestly more realistic setting the clear-sky problem do The example uses the first of the Garand atmosphere used for developing RRTMGP, as described in the [paper](https://doi.org/10.1029/2019MS001621) documenting the code, repeats the column a user-specified number of times, computes the optical properties of an arbitrary cloud in each column, and computes the broadband fluxes. Fractional cloudiness is not considered, and the clouds are extensive but quite boring, with uniform condensate amount and particle size everywhere (though with different values for liquid and ice). -Note that this example is run, and the results checked automatically, when `make` is invoked in the root directory. +Note that this example is run, and the results checked automatically, when `make` is invoked in the root directory. diff --git a/examples/all-sky/make_problem_size_loop.py b/examples/all-sky/make_problem_size_loop.py index 8b4cdee98..587dd8c76 100644 --- a/examples/all-sky/make_problem_size_loop.py +++ b/examples/all-sky/make_problem_size_loop.py @@ -16,7 +16,7 @@ if __name__ == '__main__': parser = argparse.ArgumentParser( description="Description here ") - # Argument parseing described at + # Argument parseing described at # https://stackoverflow.com/questions/15753701/how-can-i-pass-a-list-as-a-command-line-argument-with-argparse parser.add_argument("-x", "--executable", type=str, diff --git a/examples/all-sky/mo_load_aerosol_coefficients.F90 b/examples/all-sky/mo_load_aerosol_coefficients.F90 index f39643fff..eede4c96e 100644 --- a/examples/all-sky/mo_load_aerosol_coefficients.F90 +++ b/examples/all-sky/mo_load_aerosol_coefficients.F90 @@ -5,7 +5,7 @@ module mo_load_aerosol_coefficients ty_optical_props_1scl, & ty_optical_props_2str, & ty_optical_props_nstr - use mo_aerosol_optics_rrtmgp_merra, & + use mo_aerosol_optics_rrtmgp_merra, & only: ty_aerosol_optics_rrtmgp_merra use mo_simple_netcdf, only: read_field, read_string, var_exists, get_dim_size, & write_field, create_dim, create_var @@ -23,7 +23,7 @@ module mo_load_aerosol_coefficients ! read aerosol optical property LUT coefficients from NetCDF file ! subroutine load_aero_lutcoeff(aerosol_spec, aero_coeff_file) - class(ty_aerosol_optics_rrtmgp_merra), & + class(ty_aerosol_optics_rrtmgp_merra), & intent(inout) :: aerosol_spec character(len=*), intent(in ) :: aero_coeff_file ! ----------------- @@ -95,8 +95,8 @@ end subroutine load_aero_lutcoeff !-------------------------------------------------------------------------------------------------------------------- subroutine read_aero_state(filename, p_lay, p_lev, t_lay, vmr_h2o) character(len=*), intent(in ) :: filename - real(wp), dimension(:,:), allocatable, intent(inout) :: p_lay ! layer pressure - real(wp), dimension(:,:), allocatable, intent(inout) :: p_lev ! level pressure + real(wp), dimension(:,:), allocatable, intent(inout) :: p_lay ! layer pressure + real(wp), dimension(:,:), allocatable, intent(inout) :: p_lev ! level pressure real(wp), dimension(:,:), allocatable, intent(inout) :: t_lay ! layer temperature real(wp), dimension(:,:), allocatable, intent(inout) :: vmr_h2o ! water volume mixing ratio diff --git a/examples/all-sky/mo_load_cloud_coefficients.F90 b/examples/all-sky/mo_load_cloud_coefficients.F90 index 10e40cafc..a79469bef 100644 --- a/examples/all-sky/mo_load_cloud_coefficients.F90 +++ b/examples/all-sky/mo_load_cloud_coefficients.F90 @@ -5,7 +5,7 @@ module mo_load_cloud_coefficients ty_optical_props_1scl, & ty_optical_props_2str, & ty_optical_props_nstr - use mo_cloud_optics_rrtmgp, & + use mo_cloud_optics_rrtmgp, & only: ty_cloud_optics_rrtmgp use mo_simple_netcdf, only: read_field, read_string, var_exists, get_dim_size, & write_field, create_dim, create_var diff --git a/examples/all-sky/rrtmgp_allsky.F90 b/examples/all-sky/rrtmgp_allsky.F90 index f891a194e..53bb547f4 100644 --- a/examples/all-sky/rrtmgp_allsky.F90 +++ b/examples/all-sky/rrtmgp_allsky.F90 @@ -1,5 +1,5 @@ program rte_rrtmgp_allsky - use, intrinsic :: iso_fortran_env, & + use, intrinsic :: iso_fortran_env, & only: output_unit use mo_rte_kind, only: wp, i8, wl use mo_optical_props, only: ty_optical_props, & @@ -38,7 +38,7 @@ program rte_rrtmgp_allsky real(wp), dimension(:), allocatable :: mu0 real(wp), dimension(:,:), allocatable :: sfc_alb_dir, sfc_alb_dif ! First dimension is band ! - ! Gas concentrations + ! Gas concentrations ! character(len=3), dimension(8), parameter :: & gas_names = ['h2o', 'co2', 'o3 ', 'n2o', 'co ', 'ch4', 'o2 ', 'n2 '] @@ -58,7 +58,7 @@ program rte_rrtmgp_allsky ! Aerosols ! logical :: cell_has_aerosols - integer, dimension(:,:), allocatable :: aero_type + integer, dimension(:,:), allocatable :: aero_type ! MERRA2/GOCART aerosol type real(wp), dimension(:,:), allocatable :: aero_size ! Aerosol size for dust and sea salt @@ -79,7 +79,7 @@ program rte_rrtmgp_allsky ! type(ty_gas_optics_rrtmgp) :: k_dist type(ty_cloud_optics_rrtmgp) :: cloud_optics - type(ty_aerosol_optics_rrtmgp_merra) & + type(ty_aerosol_optics_rrtmgp_merra) & :: aerosol_optics type(ty_gas_concs) :: gas_concs, gas_concs_garand, gas_concs_1col class(ty_optical_props_arry), & @@ -93,12 +93,12 @@ program rte_rrtmgp_allsky integer :: nbnd, ngpt integer :: icol, ilay, ibnd, iloop, igas - logical :: top_is_at_1 ! CCE OMP workaround + logical :: top_is_at_1 ! CCE OMP workaround character(len=8) :: char_input integer :: nUserArgs, nloops, ncol, nlay ! logical :: write_fluxes = .false. - logical :: do_clouds = .false., use_luts = .true. + logical :: do_clouds = .false., use_luts = .true. logical :: do_aerosols = .false. integer, parameter :: ngas = 8 @@ -115,7 +115,7 @@ program rte_rrtmgp_allsky ! Code ! ---------------------------------------------------------------------------------- ! - ! Parse command line: rrtmgp_allsky ncol nlay nreps kdist [clouds [aerosols]] + ! Parse command line: rrtmgp_allsky ncol nlay nreps kdist [clouds [aerosols]] ! nUserArgs = command_argument_count() @@ -136,14 +136,14 @@ program rte_rrtmgp_allsky call get_command_argument(4,output_file) call get_command_argument(5,k_dist_file) - if (nUserArgs >= 6) then + if (nUserArgs >= 6) then call get_command_argument(6,cloud_optics_file) - do_clouds = .true. - end if - if (nUserArgs >= 7) then + do_clouds = .true. + end if + if (nUserArgs >= 7) then call get_command_argument(7,aerosol_optics_file) - do_aerosols = .true. - end if + do_aerosols = .true. + end if if (nUserArgs > 7) print *, "Ignoring command line arguments beyond the first seven..." ! ----------------------------------------------------------------------------------- allocate(p_lay(ncol, nlay), t_lay(ncol, nlay), p_lev(ncol, nlay+1), t_lev(ncol, nlay+1)) @@ -153,20 +153,20 @@ program rte_rrtmgp_allsky call compute_profiles(300._wp, ncol, nlay, p_lay, t_lay, p_lev, t_lev, q, o3) call stop_on_err(gas_concs%init(gas_names)) - call stop_on_err(gas_concs%set_vmr("h2o", q )) - call stop_on_err(gas_concs%set_vmr("o3", o3)) - call stop_on_err(gas_concs%set_vmr("co2", 348.e-6_wp)) - call stop_on_err(gas_concs%set_vmr("ch4", 1650.e-9_wp)) - call stop_on_err(gas_concs%set_vmr("n2o", 306.e-9_wp)) - call stop_on_err(gas_concs%set_vmr("n2", 0.7808_wp)) - call stop_on_err(gas_concs%set_vmr("o2", 0.2095_wp)) - call stop_on_err(gas_concs%set_vmr("co", 0._wp)) + call stop_on_err(gas_concs%set_vmr("h2o", q )) + call stop_on_err(gas_concs%set_vmr("o3", o3)) + call stop_on_err(gas_concs%set_vmr("co2", 348.e-6_wp)) + call stop_on_err(gas_concs%set_vmr("ch4", 1650.e-9_wp)) + call stop_on_err(gas_concs%set_vmr("n2o", 306.e-9_wp)) + call stop_on_err(gas_concs%set_vmr("n2", 0.7808_wp)) + call stop_on_err(gas_concs%set_vmr("o2", 0.2095_wp)) + call stop_on_err(gas_concs%set_vmr("co", 0._wp)) ! ---------------------------------------------------------------------------- ! load data into classes call load_and_init(k_dist, k_dist_file, gas_concs) is_sw = k_dist%source_is_external() is_lw = .not. is_sw - if (do_clouds) then + if (do_clouds) then ! ! Should also try with Pade calculations ! call load_cld_padecoeff(cloud_optics, cloud_optics_file) @@ -179,12 +179,12 @@ program rte_rrtmgp_allsky call stop_on_err(cloud_optics%set_ice_roughness(2)) end if - if (do_aerosols) then + if (do_aerosols) then ! ! Load aerosol optics coefficients from lookup tables ! call load_aero_lutcoeff (aerosol_optics, aerosol_optics_file) - end if + end if ! ---------------------------------------------------------------------------- ! @@ -249,7 +249,7 @@ program rte_rrtmgp_allsky !$acc enter data create (t_sfc, emis_sfc) !$omp target enter data map(alloc:t_sfc, emis_sfc) ! Surface temperature - top_is_at_1 = atmos%top_is_at_1() ! CCE OMP workaround + top_is_at_1 = atmos%top_is_at_1() ! CCE OMP workaround !$acc kernels !$omp target t_sfc = t_lev(1, merge(nlay+1, 1, top_is_at_1)) @@ -294,12 +294,12 @@ program rte_rrtmgp_allsky ! ! Cloud optics ! - if(do_clouds) & + if(do_clouds) & call stop_on_err(cloud_optics%cloud_optics(lwp, iwp, rel, rei, clouds)) ! ! Aerosol optics ! - if(do_aerosols) & + if(do_aerosols) & call stop_on_err(aerosol_optics%aerosol_optics(aero_type, aero_size, & aero_mass, relhum, aerosols)) ! @@ -309,8 +309,8 @@ program rte_rrtmgp_allsky fluxes%flux_dn => flux_dn(:,:) if(is_lw) then ! - ! Should we allocate these once, rather than once per loop? They're big. - ! + ! Should we allocate these once, rather than once per loop? They're big. + ! !$acc data create( lw_sources, lw_sources%lay_source, lw_sources%lev_source) & !$acc create( lw_sources%sfc_source, lw_sources%sfc_source_Jac) !$omp target data map(alloc: lw_sources%lay_source, lw_sources%lev_source) & @@ -340,18 +340,18 @@ program rte_rrtmgp_allsky gas_concs, & atmos, & toa_flux)) - if(do_clouds) then + if(do_clouds) then call stop_on_err(clouds%delta_scale()) call stop_on_err(clouds%increment(atmos)) - end if - if(do_aerosols) then + end if + if(do_aerosols) then call stop_on_err(aerosols%delta_scale()) call stop_on_err(aerosols%increment(atmos)) - end if + end if call stop_on_err(rte_sw(atmos, mu0, toa_flux, & sfc_alb_dir, sfc_alb_dif, & fluxes)) - !$acc end data + !$acc end data !$omp end target data end if call system_clock(finish, clock_rate) @@ -361,20 +361,20 @@ program rte_rrtmgp_allsky call system_clock(finish_all, clock_rate) avg = sum( elapsed(merge(2,1,nloops>1):) ) / real(merge(nloops-1,nloops,nloops>1)) - mint = minval(elapsed) + mint = minval(elapsed) - ! What to print? - ! ncol, nlay, ngpt; are clouds used, are aerosols used; time per column, total, min; + ! What to print? + ! ncol, nlay, ngpt; are clouds used, are aerosols used; time per column, total, min; print *, " ncol nlay ngpt clouds aerosols time_per_col_ms nloops time_total_s time_min_s" - write(output_unit, '(3(i6, 1x), 6x, 2(i1, 8x), 1x, f7.3, 1x, i6, 2x, 2(4x,f7.3))') & - ncol, nlay, ngpt, merge(1,0,do_clouds), merge(1,0,do_aerosols), & + write(output_unit, '(3(i6, 1x), 6x, 2(i1, 8x), 1x, f7.3, 1x, i6, 2x, 2(4x,f7.3))') & + ncol, nlay, ngpt, merge(1,0,do_clouds), merge(1,0,do_aerosols), & avg/(real(ncol) * (1.0e-3*clock_rate)), nloops, sum(elapsed) / real(clock_rate), mint / real(clock_rate) call write_fluxes - ! - ! Memory for bounday conditions on the GPU was allocated with unstructured data dataments - ! (acc enter data). Deallocate it expliicity + ! + ! Memory for bounday conditions on the GPU was allocated with unstructured data dataments + ! (acc enter data). Deallocate it expliicity ! if(is_lw) then !$acc exit data delete( t_sfc, emis_sfc) @@ -383,9 +383,9 @@ program rte_rrtmgp_allsky !$acc exit data delete( sfc_alb_dir, sfc_alb_dif, mu0) !$omp target exit data map(release:sfc_alb_dir, sfc_alb_dif, mu0) end if - + ! - ! Clouds and aerosols also used enter data + ! Clouds and aerosols also used enter data ! if(do_clouds) then #ifndef _CRAYFTN @@ -402,9 +402,9 @@ program rte_rrtmgp_allsky !$acc exit data delete (clouds%tau, clouds%ssa, clouds%g, clouds) !$omp target exit data map(release:clouds%tau, clouds%ssa, clouds%g) end select - ! - ! Explicit finalization of cloud optical properties - not really necessary since memory - ! will be freed when the program ends, but useful for testing + ! + ! Explicit finalization of cloud optical properties - not really necessary since memory + ! will be freed when the program ends, but useful for testing ! call clouds%finalize end if @@ -419,9 +419,9 @@ program rte_rrtmgp_allsky !$acc exit data delete (aerosols%tau, aerosols%ssa, aerosols%g, aerosols) !$omp target exit data map(release:aerosols%tau, aerosols%ssa, aerosols%g) end select - ! - ! Explicit finalization of aerosol optical properties - not really necessary since memory - ! will be freed when the program ends, but useful for testing + ! + ! Explicit finalization of aerosol optical properties - not really necessary since memory + ! will be freed when the program ends, but useful for testing ! call aerosols%finalize end if @@ -429,29 +429,29 @@ program rte_rrtmgp_allsky ! k-distribution ! call k_dist%finalize - + if(.not. is_lw) then !$acc exit data delete( flux_dir) !$omp target exit data map(release:flux_dir) end if - ! fluxes - but not flux_dir, which used enter data - !$acc end data + ! fluxes - but not flux_dir, which used enter data + !$acc end data !$omp end target data ! p_lay etc - !$acc end data + !$acc end data !$omp end target data contains ! ---------------------------------------------------------------------------------- subroutine compute_profiles(SST, ncol, nlay, p_lay, t_lay, p_lev, t_lev, q_lay, o3) ! - ! Construct profiles of pressure, temperature, humidity, and ozone + ! Construct profiles of pressure, temperature, humidity, and ozone ! more or less following the RCEMIP protocol for a surface temperature of 300K ! more or less follows a Python implementation by Chiel van Heerwardeen - ! Extensions for future - variable SST and T profile, variable RH, lapse rate in stratosphere - ! will all access absorption coefficient data more realistically + ! Extensions for future - variable SST and T profile, variable RH, lapse rate in stratosphere + ! will all access absorption coefficient data more realistically ! - real(wp), intent(in ) :: SST + real(wp), intent(in ) :: SST integer, intent(in ) :: ncol, nlay real(wp), dimension(ncol, nlay ), intent(out) :: p_lay, t_lay, q_lay, o3 real(wp), dimension(ncol, nlay+1), intent(out) :: p_lev, t_lev @@ -462,76 +462,76 @@ subroutine compute_profiles(SST, ncol, nlay, p_lay, t_lay, p_lev, t_lev, q_lay, integer :: icol, ilay, i real(wp), parameter :: z_trop = 15000._wp, z_top = 70.e3_wp - ! Ozone profile - maybe only a single profile? - real(wp), parameter :: g1 = 3.6478_wp, g2 = 0.83209_wp, g3 = 11.3515_wp, o3_min = 1e-13_wp + ! Ozone profile - maybe only a single profile? + real(wp), parameter :: g1 = 3.6478_wp, g2 = 0.83209_wp, g3 = 11.3515_wp, o3_min = 1e-13_wp ! According to CvH RRTMGP in Single Precision will fail with lower ozone concentrations - real(wp), parameter :: g = 9.79764, Rd = 287.04, p0 = 101480. ! Surface pressure + real(wp), parameter :: g = 9.79764, Rd = 287.04, p0 = 101480. ! Surface pressure real(wp), parameter :: z_q1 = 4.0e3, z_q2 = 7.5e3, q_t = 1.e-8 real(wp), parameter :: gamma = 6.7e-3 - + real(wp), parameter :: q_0 = 0.01864 ! for 300 K SST. ! ------------------- Tv0 = (1. + 0.608*q_0) * SST ! ! Split resolution above and below RCE tropopause (15 km or about 125 hPa) ! - z_lev(:) = [0._wp, 2._wp* z_trop /nlay * [(i, i=1, nlay/2)], & + z_lev(:) = [0._wp, 2._wp* z_trop /nlay * [(i, i=1, nlay/2)], & z_trop + 2._wp * (z_top - z_trop)/nlay * [(i, i=1, nlay/2)]] z_lay(:) = 0.5_wp * (z_lev(1:nlay) + z_lev(2:nlay+1)) - - !$acc data copyin(z_lev, z_lay) + + !$acc data copyin(z_lev, z_lay) !$omp target data map(to:z_lev, z_lay) ! - ! The two loops are the same, except applied to layers and levels + ! The two loops are the same, except applied to layers and levels ! but nvfortran doesn't seems to support elemental procedures in OpenACC loops ! - !$acc parallel loop collapse(2) - !$omp target teams distribute parallel do simd collapse(2) - do ilay = 1, nlay - do icol = 1, ncol - z = z_lay(ilay) - if (z > z_trop) then + !$acc parallel loop collapse(2) + !$omp target teams distribute parallel do simd collapse(2) + do ilay = 1, nlay + do icol = 1, ncol + z = z_lay(ilay) + if (z > z_trop) then q = q_t T = SST - gamma*z_trop/(1. + 0.608*q_0) Tv = (1. + 0.608*q ) * T p = p0 * (Tv/Tv0)**(g/(Rd*gamma)) * exp( -((g*(z-z_trop))/(Rd*Tv)) ) - else + else q = q_0 * exp(-z/z_q1) * exp(-(z/z_q2)**2) - T = SST - gamma*z / (1. + 0.608*q) + T = SST - gamma*z / (1. + 0.608*q) Tv = (1. + 0.608*q ) * T p = p0 * (Tv/Tv0)**(g/(Rd*gamma)) - end if - p_lay(icol,ilay) = p + end if + p_lay(icol,ilay) = p t_lay(icol,ilay) = T q_lay(icol,ilay) = q p_hpa = p_lay(icol,ilay) / 100._wp - o3(icol, ilay) = max(o3_min, & + o3(icol, ilay) = max(o3_min, & g1 * p_hpa**g2 * exp(-p_hpa/g3) * 1.e-6_wp) end do - end do + end do - !$acc parallel loop collapse(2) - !$omp target teams distribute parallel do simd collapse(2) + !$acc parallel loop collapse(2) + !$omp target teams distribute parallel do simd collapse(2) do ilay = 1, nlay+1 - do icol = 1, ncol - z = z_lev(ilay) - if (z > z_trop) then + do icol = 1, ncol + z = z_lev(ilay) + if (z > z_trop) then q = q_t T = SST - gamma*z_trop/(1. + 0.608*q_0) Tv = (1. + 0.608*q ) * T p = p0 * (Tv/Tv0)**(g/(Rd*gamma)) * exp( -((g*(z-z_trop))/(Rd*Tv)) ) - else + else q = q_0 * exp(-z/z_q1) * exp(-(z/z_q2)**2) - T = SST - gamma*z / (1. + 0.608*q) + T = SST - gamma*z / (1. + 0.608*q) Tv = (1. + 0.608*q ) * T p = p0 * (Tv/Tv0)**(g/(Rd*gamma)) - end if - p_lev(icol,ilay) = p + end if + p_lev(icol,ilay) = p t_lev(icol,ilay) = T - end do - end do + end do + end do !$acc end data !$omp end target data end subroutine compute_profiles @@ -548,10 +548,10 @@ subroutine stop_on_err(error_msg) end subroutine stop_on_err ! -------------------------------------------------------------------------------------- ! - subroutine compute_clouds + subroutine compute_clouds real(wp) :: rel_val, rei_val - ! - ! Variable and memory allocation + ! + ! Variable and memory allocation ! if(is_sw) then allocate(ty_optical_props_2str::clouds) @@ -574,7 +574,7 @@ subroutine compute_clouds call stop_on_err("rte_rrtmgp_allsky: Don't recognize the kind of optical properties ") end select ! - ! Cloud physical properties + ! Cloud physical properties ! allocate(lwp(ncol,nlay), iwp(ncol,nlay), & rel(ncol,nlay), rei(ncol,nlay), cloud_mask(ncol,nlay)) @@ -605,16 +605,16 @@ subroutine compute_clouds end do !$acc exit data delete(cloud_mask) !$omp target exit data map(release:cloud_mask) - + end subroutine compute_clouds ! ! -------------------------------------------------------------------------------------- ! subroutine compute_aerosols real(wp), dimension(ncol,nlay) :: vmr_h2o ! h2o vmr - logical :: is_sulfate, is_dust, is_even_column - ! - ! Variable and memory allocation + logical :: is_sulfate, is_dust, is_even_column + ! + ! Variable and memory allocation ! if(is_sw) then allocate(ty_optical_props_2str::aerosols) @@ -642,8 +642,8 @@ subroutine compute_aerosols !$omp target data map(alloc:vmr_h2o) call stop_on_err(gas_concs%get_vmr("h2o",vmr_h2o)) ! - ! Aerosol properties - ! + ! Aerosol properties + ! allocate(aero_type(ncol,nlay), aero_size(ncol,nlay), & aero_mass(ncol,nlay), relhum (ncol,nlay)) !$acc enter data create( aero_type, aero_size, aero_mass, relhum) @@ -657,20 +657,20 @@ subroutine compute_aerosols ! put them in 1/2 of the columns ! ! - !$acc parallel loop collapse(2) copyin(p_lay) - !$omp target teams distribute parallel do simd collapse(2) map(to:p_lay) + !$acc parallel loop collapse(2) copyin(p_lay) + !$omp target teams distribute parallel do simd collapse(2) map(to:p_lay) do ilay=1,nlay do icol=1,ncol - is_sulfate = (p_lay(icol,ilay) > 50._wp * 100._wp .and. & + is_sulfate = (p_lay(icol,ilay) > 50._wp * 100._wp .and. & p_lay(icol,ilay) < 100._wp * 100._wp) - is_dust = (p_lay(icol,ilay) > 700._wp * 100._wp .and. & + is_dust = (p_lay(icol,ilay) > 700._wp * 100._wp .and. & p_lay(icol,ilay) < 900._wp * 100._wp) is_even_column = mod(icol, 2) /= 0 - if (is_even_column .and. is_sulfate) then + if (is_even_column .and. is_sulfate) then aero_type(icol,ilay) = merra_aero_sulf aero_size(icol,ilay) = 0.2_wp aero_mass(icol,ilay) = 1.e-6_wp - else if(is_even_column .and. is_dust) then + else if(is_even_column .and. is_dust) then ! Dust aerosol aero_type(icol,ilay) = merra_aero_dust aero_size(icol,ilay) = 0.5_wp @@ -699,7 +699,7 @@ subroutine get_relhum(ncol, nlay, p_lay, t_lay, vmr_h2o, relhum) real(wp), intent(inout) :: relhum(ncol,nlay) ! relative humidity (fraction, 0-1) - ! Local variables + ! Local variables integer :: i, k real(wp) :: mmr_h2o ! water mass mixing ratio @@ -715,8 +715,8 @@ subroutine get_relhum(ncol, nlay, p_lay, t_lay, vmr_h2o, relhum) ! Derive layer virtual temperature !$acc parallel loop collapse(2) copyin(p_lay, vmr_h2o, t_lay) copyout( relhum) - !$omp target teams distribute parallel do simd collapse(2) map(to:p_lay, vmr_h2o, t_lay) map(from:relhum) - do i = 1, ncol + !$omp target teams distribute parallel do simd collapse(2) map(to:p_lay, vmr_h2o, t_lay) map(from:relhum) + do i = 1, ncol do k = 1, nlay ! Convert h2o vmr to mmr mmr_h2o = vmr_h2o(i,k) * mwd @@ -730,18 +730,18 @@ subroutine get_relhum(ncol, nlay, p_lay, t_lay, vmr_h2o, relhum) enddo end subroutine get_relhum !-------------------------------------------------------------------------------------------------------------------- - subroutine write_fluxes + subroutine write_fluxes use netcdf use mo_simple_netcdf, only: write_field integer :: ncid, i, col_dim, lay_dim, lev_dim, varid real(wp) :: vmr(ncol, nlay) character(len=3) :: flux_prefix ! - ! Write fluxes - make this optional? + ! Write fluxes - make this optional? ! ! - ! Define dimensions + ! Define dimensions ! if(nf90_create(trim(output_file), NF90_CLOBBER, ncid) /= NF90_NOERR) & call stop_on_err("rrtmgp_allsky: can't create file " // trim(output_file)) @@ -752,11 +752,11 @@ subroutine write_fluxes call stop_on_err("rrtmgp_allsky: can't define lay dimension") if(nf90_def_dim(ncid, "lev", nlay+1, lev_dim) /= NF90_NOERR) & call stop_on_err("rrtmgp_allsky: can't define lev dimension") - + ! - ! Define variables + ! Define variables ! - ! State + ! State ! call create_var("p_lev", ncid, [col_dim, lev_dim]) call create_var("t_lev", ncid, [col_dim, lev_dim]) @@ -765,29 +765,29 @@ subroutine write_fluxes call create_var("h2o", ncid, [col_dim, lay_dim]) call create_var("o3", ncid, [col_dim, lay_dim]) - ! All the gases except h2o, o3 - write as attributes? Or not bother? + ! All the gases except h2o, o3 - write as attributes? Or not bother? - if(do_clouds) then + if(do_clouds) then call create_var("lwp", ncid, [col_dim, lay_dim]) call create_var("iwp", ncid, [col_dim, lay_dim]) call create_var("rel", ncid, [col_dim, lay_dim]) call create_var("rei", ncid, [col_dim, lay_dim]) - end if - if(do_aerosols) then + end if + if(do_aerosols) then if(nf90_def_var(ncid, "aero_type", NF90_SHORT, [col_dim, lay_dim], varid) /= NF90_NOERR) & call stop_on_err("create_var: can't define variable aero_type") call create_var("aero_size", ncid, [col_dim, lay_dim]) call create_var("aero_mass", ncid, [col_dim, lay_dim]) - end if + end if ! - ! Fluxes - definitions + ! Fluxes - definitions ! - if(is_sw) then + if(is_sw) then flux_prefix = "sw_" call create_var(flux_prefix // "flux_dir", ncid, [col_dim, lev_dim]) else - flux_prefix = "lw_" - end if + flux_prefix = "lw_" + end if call create_var(flux_prefix // "flux_up", ncid, [col_dim, lev_dim]) call create_var(flux_prefix //"flux_dn", ncid, [col_dim, lev_dim]) if(nf90_enddef(ncid) /= NF90_NOERR) & @@ -796,7 +796,7 @@ subroutine write_fluxes ! ! Write variables ! - ! State - writing + ! State - writing !$acc update host(p_lev, t_lev, p_lay, t_lay) !$omp target update from(p_lev, t_lev, p_lay, t_lay) call stop_on_err(write_field(ncid, "p_lev", p_lev)) @@ -809,43 +809,43 @@ subroutine write_fluxes call stop_on_err(gas_concs%get_vmr("o3", vmr)) call stop_on_err(write_field(ncid, "o3", vmr)) - if(do_clouds) then + if(do_clouds) then !$acc update host(lwp, iwp, rel, rei) !$omp target update from(lwp, iwp, rel, rei) call stop_on_err(write_field(ncid, "lwp", lwp)) call stop_on_err(write_field(ncid, "iwp", iwp)) call stop_on_err(write_field(ncid, "rel", rel)) call stop_on_err(write_field(ncid, "rei", rei)) - end if + end if - if(do_aerosols) then + if(do_aerosols) then !$acc update host(aero_size, aero_mass, aero_type) !$omp target update from(aero_size, aero_mass, aero_type) call stop_on_err(write_field(ncid, "aero_size", aero_size)) call stop_on_err(write_field(ncid, "aero_mass", aero_mass)) call stop_on_err(write_field(ncid, "aero_type", aero_type)) - end if + end if - ! Fluxes - writing + ! Fluxes - writing !$acc update host(flux_up, flux_dn) !$omp target update from(flux_up, flux_dn) call stop_on_err(write_field(ncid, flux_prefix // "flux_up", flux_up)) call stop_on_err(write_field(ncid, flux_prefix // "flux_dn", flux_dn)) - if(.not. is_lw) then + if(.not. is_lw) then !$acc update host(flux_dir) !$omp target update from(flux_dir) call stop_on_err(write_field(ncid, flux_prefix // "flux_dir", flux_dir)) - end if + end if - ! Close netCDF + ! Close netCDF if(nf90_close(ncid) /= NF90_NOERR) call stop_on_err("rrtmgp_allsky: error closing file??") end subroutine write_fluxes ! --------------------------------------------------------- subroutine create_var(name, ncid, dim_ids) use netcdf - character(len=*), intent(in) :: name + character(len=*), intent(in) :: name integer, intent(in) :: ncid - integer, dimension(:), intent(in) :: dim_ids + integer, dimension(:), intent(in) :: dim_ids integer :: varid diff --git a/examples/compare-to-reference.py b/examples/compare-to-reference.py index bb7d847a3..6cb1ce9b0 100644 --- a/examples/compare-to-reference.py +++ b/examples/compare-to-reference.py @@ -1,13 +1,13 @@ #! /usr/bin/env python # -# General purpose comparison script -- compare all variables in a set of files, -# write output if differences exceed some threshold, -# error exit if differences exceed a different threshold +# General purpose comparison script -- compare all variables in a set of files, +# write output if differences exceed some threshold, +# error exit if differences exceed a different threshold +# +# Currently thresholds are specified as absolute differences +# worth revising but could change development practice +# Thresholds come from environement variables when set? # -# Currently thresholds are specified as absolute differences -# worth revising but could change development practice -# Thresholds come from environement variables when set? -# # import sys @@ -62,7 +62,7 @@ raise Exception( v + ": some test values are missing. Now that is strange.") # - # Reporting + # Reporting # if not np.allclose(tst[v], ref[v], atol=args.report_threshold, rtol=0): diff = abs((tst - ref)[v].values) diff --git a/examples/rfmip-clear-sky/CMakeLists.txt b/examples/rfmip-clear-sky/CMakeLists.txt index 37c69177a..8026f3afc 100644 --- a/examples/rfmip-clear-sky/CMakeLists.txt +++ b/examples/rfmip-clear-sky/CMakeLists.txt @@ -21,7 +21,8 @@ set(inoutputs list( TRANSFORM inoutputs PREPEND ${RRTMGP_DATA}/examples/rfmip-clear-sky/inputs/ - OUTPUT_VARIABLE inputs + OUTPUT_VARIABLE + inputs ) # The tests write to the input files, therefore we copy them: diff --git a/examples/rfmip-clear-sky/Readme.md b/examples/rfmip-clear-sky/Readme.md index f3a82a5db..24520ae27 100644 --- a/examples/rfmip-clear-sky/Readme.md +++ b/examples/rfmip-clear-sky/Readme.md @@ -5,7 +5,7 @@ the [RTE+RRTMGP](https://github.com/RobertPincus/rte-rrtmgp) radiation parameter Note that this example is run, and the results checked automatically, when `make` is invoked in the root directory. -By default needed files (input conditions, output templates) are copied from a subdirectory of `$RRTMGP_DATA`. +By default needed files (input conditions, output templates) are copied from a subdirectory of `$RRTMGP_DATA`. A Python script `stage_files.py` may also used to download relevant files from the [RFMIP web site](https://www.earthsystemcog.org/projects/rfmip/resources/).This script invokes another Python script to create empty output files. diff --git a/examples/rfmip-clear-sky/mo_rfmip_io.F90 b/examples/rfmip-clear-sky/mo_rfmip_io.F90 index 718f75de4..1d6072f4f 100644 --- a/examples/rfmip-clear-sky/mo_rfmip_io.F90 +++ b/examples/rfmip-clear-sky/mo_rfmip_io.F90 @@ -350,7 +350,7 @@ subroutine read_and_block_gases_ty(fileName, blocksize, gas_names, names_in_file call stop_on_err(gas_conc_array(b)%init(gas_names)) else call stop_on_err(gas_conc_array(b)%init([gas_names, 'h2o ', 'o3 ', 'no2 '])) - end if + end if end do ! ! Which gases are known to the k-distribution and available in the files? diff --git a/extensions/mo_cloud_sampling.F90 b/extensions/mo_cloud_sampling.F90 index f2359b096..c47ef4364 100644 --- a/extensions/mo_cloud_sampling.F90 +++ b/extensions/mo_cloud_sampling.F90 @@ -195,7 +195,7 @@ end function sampled_mask_max_ran ! Generate a McICA-sampled cloud mask for exponential-random overlap. ! The overlap parameter overlap_param is defined between pairs of layers. ! For layer i, overlap_param(i) describes the overlap between cloud_frac(i) and cloud_frac(i+1). - ! It is a correlation coefficient in [-1,1]. E.g., + ! It is a correlation coefficient in [-1,1]. E.g., ! +1 gives perfect correlation or maximum cloud overlap between layers i & i+1; ! 0 gives no correlation or random cloud overlap between layers i & i+1; ! -1 gives perfect anticorrelation or minimum cloud overlap between layers i & i+1. diff --git a/extensions/mo_compute_bc.F90 b/extensions/mo_compute_bc.F90 index fbec1d974..9a7f4d02a 100644 --- a/extensions/mo_compute_bc.F90 +++ b/extensions/mo_compute_bc.F90 @@ -21,7 +21,7 @@ module mo_compute_bc ! ------------------------------------------------------------------------------------------------- use mo_rte_kind, only: wp, wl use mo_rte_config, only: check_extents - use mo_rte_util_array_validation, & + use mo_rte_util_array_validation, & only: extents_are use mo_source_functions, only: ty_source_func_lw use mo_gas_concentrations, only: ty_gas_concs diff --git a/extensions/mo_fluxes_byband.F90 b/extensions/mo_fluxes_byband.F90 index 7bf334b66..733c9cadd 100644 --- a/extensions/mo_fluxes_byband.F90 +++ b/extensions/mo_fluxes_byband.F90 @@ -18,7 +18,7 @@ module mo_fluxes_byband use mo_rte_kind, only: wp use mo_rte_config, only: check_extents - use mo_rte_util_array_validation, & + use mo_rte_util_array_validation, & only: extents_are use mo_fluxes, only: ty_fluxes, ty_fluxes_broadband use mo_optical_props, only: ty_optical_props @@ -151,7 +151,7 @@ function are_desired_byband(this) end function are_desired_byband ! ---------------------------------------------------------------------------- - ! Kernels (private to this module) + ! Kernels (private to this module) ! ---------------------------------------------------------------------------- ! ! Spectral reduction over all points diff --git a/extensions/mo_fluxes_bygpoint.F90 b/extensions/mo_fluxes_bygpoint.F90 index 20130f568..bdb35850d 100644 --- a/extensions/mo_fluxes_bygpoint.F90 +++ b/extensions/mo_fluxes_bygpoint.F90 @@ -16,7 +16,7 @@ ! module mo_fluxes_bygpoint use mo_rte_kind, only: wp - use mo_rte_util_array_validation, & + use mo_rte_util_array_validation, & only: extents_are use mo_fluxes, only: ty_fluxes use mo_optical_props, only: ty_optical_props diff --git a/extensions/mo_heating_rates.F90 b/extensions/mo_heating_rates.F90 index 9535ae4cd..90ebcee2e 100644 --- a/extensions/mo_heating_rates.F90 +++ b/extensions/mo_heating_rates.F90 @@ -16,9 +16,9 @@ module mo_heating_rates use mo_rte_kind, only: wp, wl use mo_rte_config, only: check_extents - use mo_rte_util_array_validation, & + use mo_rte_util_array_validation, & only: extents_are, any_vals_less_than - use mo_gas_optics_constants, & + use mo_gas_optics_constants, & only: cp_dry, grav ! Only needed for heating rate calculation implicit none private diff --git a/extensions/mo_zenith_angle_spherical_correction.F90 b/extensions/mo_zenith_angle_spherical_correction.F90 index 232af9ff6..134fccd01 100644 --- a/extensions/mo_zenith_angle_spherical_correction.F90 +++ b/extensions/mo_zenith_angle_spherical_correction.F90 @@ -17,7 +17,7 @@ module mo_zenith_angle_spherical_correction use mo_rte_kind, only: wp, wl use mo_rte_config, only: check_extents, check_values - use mo_rte_util_array_validation, & + use mo_rte_util_array_validation, & only: extents_are, any_vals_outside, any_vals_less_than implicit none private diff --git a/gas-optics/mo_gas_concentrations.F90 b/gas-optics/mo_gas_concentrations.F90 index 592baf57d..f83fb8ecf 100644 --- a/gas-optics/mo_gas_concentrations.F90 +++ b/gas-optics/mo_gas_concentrations.F90 @@ -16,7 +16,7 @@ ! !> Values may be provided as scalars, 1-dimensional profiles (nlay), or 2-D fields (ncol,nlay). !> `nlay` and `ncol` are determined from the input arrays; self-consistency is enforced. -!> No bounds are enforced on the sum of the mixing ratios. +!> No bounds are enforced on the sum of the mixing ratios. !> !> For example: !> ``` @@ -29,7 +29,7 @@ !> or as 2D fields. Values for all columns are returned although the entire collection !> can be subsetted in the column dimension !> -!> Subsets can be extracted in the column dimension. +!> Subsets can be extracted in the column dimension. !> !> Functions return strings. Non-empty strings indicate an error. !> @@ -38,7 +38,7 @@ module mo_gas_concentrations use mo_rte_kind, only: wp use mo_rte_config, only: check_values - use mo_rte_util_array_validation, & + use mo_rte_util_array_validation, & only: any_vals_outside implicit none integer, parameter, private :: GAS_NOT_IN_LIST = -1 @@ -74,11 +74,11 @@ module mo_gas_concentrations procedure, public :: reset generic, public :: set_vmr => set_vmr_scalar, & set_vmr_1d, & - set_vmr_2d !! ### Set concentration values + set_vmr_2d !! ### Set concentration values generic, public :: get_vmr => get_vmr_1d, & get_vmr_2d !! ### Get concentration values - generic, public :: get_subset => get_subset_range - !! ### Extract a subset of columns + generic, public :: get_subset => get_subset_range + !! ### Extract a subset of columns procedure, public :: get_num_gases procedure, public :: get_gas_names end type ty_gas_concs @@ -87,8 +87,8 @@ module mo_gas_concentrations !> ### Initialize the object function init(this, gas_names) result(error_msg) class(ty_gas_concs), intent(inout) :: this - character(len=*), dimension(:), intent(in ) :: gas_names !! names of all gases which might be provided - character(len=128) :: error_msg !! error string, empty when successful + character(len=*), dimension(:), intent(in ) :: gas_names !! names of all gases which might be provided + character(len=128) :: error_msg !! error string, empty when successful ! --------- integer :: i, j, ngas ! --------- @@ -125,13 +125,13 @@ function init(this, gas_names) result(error_msg) ! Set concentrations --- scalar, 1D, 2D ! ! ------------------------------------------------------------------------------------- - !> ### Set scalar concentrations + !> ### Set scalar concentrations function set_vmr_scalar(this, gas, w) result(error_msg) ! In OpenACC context scalar w always assumed to be on the CPU class(ty_gas_concs), intent(inout) :: this character(len=*), intent(in ) :: gas !! Name of the gas being provided - real(wp), intent(in ) :: w !! volume (molar) mixing ratio - character(len=128) :: error_msg !! error string, empty when successful + real(wp), intent(in ) :: w !! volume (molar) mixing ratio + character(len=128) :: error_msg !! error string, empty when successful ! --------- real(wp), dimension(:,:), pointer :: p integer :: igas @@ -177,14 +177,14 @@ function set_vmr_scalar(this, gas, w) result(error_msg) !$omp end target end function set_vmr_scalar ! ------------------------------------------------------------------------------------- - !> ### Set 1d (function of level) concentrations + !> ### Set 1d (function of level) concentrations function set_vmr_1d(this, gas, w) result(error_msg) ! In OpenACC context w assumed to be either on the CPU or on the GPU class(ty_gas_concs), intent(inout) :: this character(len=*), intent(in ) :: gas !! Name of the gas being provided real(wp), dimension(:), & - intent(in ) :: w !! volume (molar) mixing ratio - character(len=128) :: error_msg !! error string, empty when successful + intent(in ) :: w !! volume (molar) mixing ratio + character(len=128) :: error_msg !! error string, empty when successful ! --------- real(wp), dimension(:,:), pointer :: p integer :: igas @@ -239,15 +239,15 @@ function set_vmr_1d(this, gas, w) result(error_msg) !$acc exit data delete(w) end function set_vmr_1d ! ------------------------------------------------------------------------------------- - !> ### Set 2d concentrations + !> ### Set 2d concentrations function set_vmr_2d(this, gas, w) result(error_msg) ! In OpenACC context w assumed to be either on the CPU or on the GPU class(ty_gas_concs), intent(inout) :: this character(len=*), intent(in ) :: gas !! Name of the gas being provided real(wp), dimension(:,:), & - intent(in ) :: w !! volume (molar) mixing ratio - character(len=128) :: error_msg - !! error string, empty when successful + intent(in ) :: w !! volume (molar) mixing ratio + character(len=128) :: error_msg + !! error string, empty when successful ! --------- real(wp), dimension(:,:), pointer :: p integer :: igas @@ -317,8 +317,8 @@ end function set_vmr_2d function get_vmr_1d(this, gas, array) result(error_msg) class(ty_gas_concs) :: this character(len=*), intent(in ) :: gas !! Name of the gas - real(wp), dimension(:), intent(out) :: array !! Volume mixing ratio - character(len=128) :: error_msg !! Error string, empty if successful + real(wp), dimension(:), intent(out) :: array !! Volume mixing ratio + character(len=128) :: error_msg !! Error string, empty if successful ! --------------------- real(wp), dimension(:,:), pointer :: p integer :: igas @@ -374,8 +374,8 @@ end function get_vmr_1d function get_vmr_2d(this, gas, array) result(error_msg) class(ty_gas_concs) :: this character(len=*), intent(in ) :: gas !! Name of the gas - real(wp), dimension(:,:), intent(out) :: array !! Volume mixing ratio - character(len=128) :: error_msg !! Error string, empty if successful + real(wp), dimension(:,:), intent(out) :: array !! Volume mixing ratio + character(len=128) :: error_msg !! Error string, empty if successful ! --------------------- real(wp), dimension(:,:), pointer :: p integer :: icol, ilay, igas @@ -451,8 +451,8 @@ end function get_vmr_2d function get_subset_range(this, start, n, subset) result(error_msg) class(ty_gas_concs), intent(in ) :: this integer, intent(in ) :: start, n !! Index of first column, number of columns to extract - class(ty_gas_concs), intent(inout) :: subset !! Object to hold the subset of columns - character(len=128) :: error_msg !! Error string, empty if successful + class(ty_gas_concs), intent(inout) :: subset !! Object to hold the subset of columns + character(len=128) :: error_msg !! Error string, empty if successful ! --------------------- real(wp), dimension(:,:), pointer :: p1, p2 integer :: i @@ -564,11 +564,11 @@ end function get_gas_names ! Private procedures ! ! ------------------------------------------------------------------------------------- - !> Convert string to lower case + !> Convert string to lower case pure function lower_case( input_string ) result( output_string ) character(len=*), intent(in) :: input_string character(len=len(input_string)) :: output_string - + ! List of character for case conversion character(len=26), parameter :: LOWER_CASE_CHARS = 'abcdefghijklmnopqrstuvwxyz' character(len=26), parameter :: UPPER_CASE_CHARS = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ' diff --git a/gas-optics/mo_gas_optics.F90 b/gas-optics/mo_gas_optics.F90 index baa17d1e8..c40738487 100644 --- a/gas-optics/mo_gas_optics.F90 +++ b/gas-optics/mo_gas_optics.F90 @@ -10,26 +10,26 @@ ! BSD 3-clause license, see http://opensource.org/licenses/BSD-3-Clause ! ------------------------------------------------------------------------------------------------- !> -!> ## Generic Fortran interface for gas optics -!> -!> Defines an interface for gas optics parameterizations: +!> ## Generic Fortran interface for gas optics !> -!> - inquiry functions for the pressure and temperature limits -!> - inquiry functions for the source of radiation (internal or planetary radiation vs. +!> Defines an interface for gas optics parameterizations: +!> +!> - inquiry functions for the pressure and temperature limits +!> - inquiry functions for the source of radiation (internal or planetary radiation vs. !> external or stellar radiation) -!> - Method for computing gas optical optical properties and incident stellar radiation -!> given pressure, temperature, and gas concentrations -!> - Method for computing gas optical optical properties and internal Planck sources -!> given pressure, temperature, and gas concentrations -!> +!> - Method for computing gas optical optical properties and incident stellar radiation +!> given pressure, temperature, and gas concentrations +!> - Method for computing gas optical optical properties and internal Planck sources +!> given pressure, temperature, and gas concentrations +!> !> This (abstract) class is a sub-classes of `ty_optical_props` in the RTE module `mo_optical_props` -!> and inherits the procedures related to spectral discratization from that class. -!> Optical properties are returned in any variable of `ty_optical_props_arry` (that is, -!> an array of values with dimensions ncol, nlay, ngpt) in the same module. -!> Internal sources of radiation are provided in a variable of type -!> `ty_source_func_lw` in RTE module `ty_source_func_lw`. -!> The module also makes use of [[mo_gas_concentrations(module):ty_gas_concs(type)]] from -!> module [[mo_gas_concentrations]]. +!> and inherits the procedures related to spectral discratization from that class. +!> Optical properties are returned in any variable of `ty_optical_props_arry` (that is, +!> an array of values with dimensions ncol, nlay, ngpt) in the same module. +!> Internal sources of radiation are provided in a variable of type +!> `ty_source_func_lw` in RTE module `ty_source_func_lw`. +!> The module also makes use of [[mo_gas_concentrations(module):ty_gas_concs(type)]] from +!> module [[mo_gas_concentrations]]. ! ! ------------------------------------------------------------------------------------------------- module mo_gas_optics @@ -78,7 +78,7 @@ function gas_optics_ext_abstract(this, & tlay !! layer temperatures [K]; (ncol,nlay) type(ty_gas_concs), intent(in ) :: gas_desc !! Gas volume mixing ratios class(ty_optical_props_arry), & - intent(inout) :: optical_props !! + intent(inout) :: optical_props !! real(wp), dimension(:,:), intent( out) :: toa_src !! Incoming solar irradiance(ncol,ngpt) character(len=128) :: error_msg !! Empty if successful ! Optional inputs @@ -105,7 +105,7 @@ function gas_optics_int_abstract(this, & intent(inout) :: optical_props !! Optical properties class(ty_source_func_lw ), & intent(inout) :: sources !! Planck sources - character(len=128) :: error_msg !! Empty if successful + character(len=128) :: error_msg !! Empty if successful real(wp), dimension(:,:), intent(in ), & optional, target :: col_dry, & !! Column dry amount (molecules/cm^2); dim(ncol,nlay) tlev !! level temperatures [K]l (ncol,nlay+1) diff --git a/gas-optics/mo_gas_optics_util_string.F90 b/gas-optics/mo_gas_optics_util_string.F90 index d3c9a14e9..9fdf9fe51 100644 --- a/gas-optics/mo_gas_optics_util_string.F90 +++ b/gas-optics/mo_gas_optics_util_string.F90 @@ -28,7 +28,7 @@ module mo_gas_optics_util_string contains ! ------------------------------------------------------------------------------------------------- - !> Convert stringo to lower case + !> Convert stringo to lower case pure function lower_case( input_string ) result( output_string ) character(len=*), intent(in) :: input_string character(len=len(input_string)) :: output_string @@ -66,7 +66,7 @@ pure function string_in_array(string, array) end function string_in_array ! -------------------------------------------------------------------------------------- ! - !> what is the location of a string in an array + !> what is the location of a string in an array ! pure function string_loc_in_array(string, array) character(len=*), intent(in) :: string diff --git a/rrtmgp-frontend/mo_aerosol_optics_rrtmgp_merra.F90 b/rrtmgp-frontend/mo_aerosol_optics_rrtmgp_merra.F90 index 11c0a3715..0b752bc0e 100644 --- a/rrtmgp-frontend/mo_aerosol_optics_rrtmgp_merra.F90 +++ b/rrtmgp-frontend/mo_aerosol_optics_rrtmgp_merra.F90 @@ -16,22 +16,22 @@ ! one sulfate type, and both hydrophobic and hydrophilic black carbon and organic carbon. ! Input aerosol optical data are stored in look-up tables. ! -! References for the gocart interactive aerosols: +! References for the gocart interactive aerosols: ! Chin et al., jgr, 2000 (https://doi.org/10.1029/2000jd900384) ! Chin et al., jas, 2002 (https://doi.org/10.1175/1520-0469(2002)059<0461:TAOTFT>2.0.CO;2) ! Colarco et al., jgr, 2010 (https://doi.org/10.1029/2009jd012820) -! -! References for merra2 aerosol reanalysis: -! Randles et al., j. clim., 2017 (https://doi.org/10.1175/jcli-d-16-0609.1) +! +! References for merra2 aerosol reanalysis: +! Randles et al., j. clim., 2017 (https://doi.org/10.1175/jcli-d-16-0609.1) ! Buchard et al., j. clim., 2017 (https://doi.org/10.1175/jcli-d-16-0613.1) -! +! ! The class can be used as-is but is also intended as an example of how to extend the RTE framework ! ------------------------------------------------------------------------------------------------- module mo_aerosol_optics_rrtmgp_merra use mo_rte_kind, only: wp, wl use mo_rte_config, only: check_extents, check_values - use mo_rte_util_array_validation,& + use mo_rte_util_array_validation,& only: extents_are, any_vals_outside use mo_optical_props, only: ty_optical_props, & ty_optical_props_arry, & @@ -103,7 +103,7 @@ function load_lut(this, band_lims_wvn, & aero_ocar_tbl, aero_ocar_rh_tbl) & result(error_msg) - class(ty_aerosol_optics_rrtmgp_merra), & + class(ty_aerosol_optics_rrtmgp_merra), & intent(inout) :: this real(wp), dimension(:,:), intent(in ) :: band_lims_wvn ! spectral discretization ! Lookup table interpolation constants @@ -167,7 +167,7 @@ function load_lut(this, band_lims_wvn, & this%aero_bcar_rh_tbl(nrh, nval, nband), & this%aero_ocar_tbl(nval, nband), & this%aero_ocar_rh_tbl(nrh, nval, nband)) - + ! Copy LUT coefficients this%merra_aero_bin_lims = merra_aero_bin_lims this%aero_rh = aero_rh @@ -203,7 +203,7 @@ subroutine finalize(this) ! Lookup table aerosol optics interpolation arrays if(allocated(this%merra_aero_bin_lims)) then deallocate(this%merra_aero_bin_lims, this%aero_rh) - !$acc exit data delete( this%merra_aero_bin_lims, this%aero_rh) + !$acc exit data delete( this%merra_aero_bin_lims, this%aero_rh) !$omp target exit data map(release:this%merra_aero_bin_lims, this%aero_rh) end if @@ -214,8 +214,8 @@ subroutine finalize(this) !$acc delete(this%aero_ocar_tbl, this%aero_ocar_rh_tbl) & !$acc delete(this) !$omp target exit data map(release:this%aero_dust_tbl, this%aero_salt_tbl, this%aero_sulf_tbl) & - !$omp map(release:this%aero_bcar_tbl, this%aero_bcar_rh_tbl) & - !$omp map(release:this%aero_ocar_tbl, this%aero_ocar_rh_tbl) + !$omp map(release:this%aero_bcar_tbl, this%aero_bcar_rh_tbl) & + !$omp map(release:this%aero_ocar_tbl, this%aero_ocar_rh_tbl) deallocate(this%aero_dust_tbl, this%aero_salt_tbl, this%aero_sulf_tbl, & this%aero_bcar_tbl, this%aero_bcar_rh_tbl, & this%aero_ocar_tbl, this%aero_ocar_rh_tbl) @@ -235,7 +235,7 @@ function aerosol_optics(this, aero_type, aero_size, aero_mass, relhum, & class(ty_aerosol_optics_rrtmgp_merra), & intent(in ) :: this - integer, intent(in ) :: aero_type(:,:) ! MERRA2/GOCART aerosol type + integer, intent(in ) :: aero_type(:,:) ! MERRA2/GOCART aerosol type ! Dimensions: (ncol,nlay) ! 1 = merra_aero_dust (dust) ! 2 = merra_aero_salt (salt) @@ -322,13 +322,13 @@ function aerosol_optics(this, aero_type, aero_size, aero_mass, relhum, & if(error_msg /= "") return end if - !$acc data copyin(aero_type, aero_size, aero_mass, relhum) - !$omp target data map(to:aero_type, aero_size, aero_mass, relhum) + !$acc data copyin(aero_type, aero_size, aero_mass, relhum) + !$omp target data map(to:aero_type, aero_size, aero_mass, relhum) ! ! Aerosol mask; don't need aerosol optics if there's no aerosol ! !$acc data create(aeromsk) - !$omp target data map(alloc:aeromsk) + !$omp target data map(alloc:aeromsk) !$acc parallel loop default(present) collapse(2) !$omp target teams distribute parallel do simd collapse(2) do ilay = 1, nlay @@ -346,13 +346,13 @@ function aerosol_optics(this, aero_type, aero_size, aero_mass, relhum, & if(any_vals_outside(relhum, aeromsk, 0._wp, 1._wp)) & error_msg = 'aerosol optics: relative humidity fraction is out of bounds' end if - ! Release aerosol mask + ! Release aerosol mask !$acc end data - !$omp end target data + !$omp end target data if(error_msg == "") then !$acc data create(atau, ataussa, ataussag) - !$omp target data map(alloc:atau, ataussa, ataussag) + !$omp target data map(alloc:atau, ataussa, ataussag) ! ! ! ---------------------------------------- @@ -418,9 +418,9 @@ function aerosol_optics(this, aero_type, aero_size, aero_mass, relhum, & end select !$acc end data !$omp end target data - end if + end if !$acc end data - !$omp end target data + !$omp end target data end function aerosol_optics !-------------------------------------------------------------------------------------------------------------------- ! @@ -429,7 +429,7 @@ end function aerosol_optics !-------------------------------------------------------------------------------------------------------------------- ! ! For size dimension, select size bin appropriate for the requested aerosol size. - ! For rh dimension, linearly interpolate values from a lookup table with "nrh" + ! For rh dimension, linearly interpolate values from a lookup table with "nrh" ! unevenly-spaced elements "aero_rh". The last dimension for all tables is band. ! Returns zero where no aerosol is present. ! @@ -469,7 +469,7 @@ subroutine compute_all_from_table(ncol, nlay, npair, nval, nrh, nbin, nbnd, do ilay = 1,nlay do icol = 1, ncol ! Sequential loop to find size bin - do i=1,nbin + do i=1,nbin if (size(icol,ilay) .ge. merra_aero_bin_lims(1,i) .and. & size(icol,ilay) .le. merra_aero_bin_lims(2,i)) then ibin = i @@ -494,7 +494,7 @@ subroutine compute_all_from_table(ncol, nlay, npair, nval, nrh, nbin, nbnd, endif endif - ! Set aerosol optical properties where aerosol present. Use aerosol type array as the mask. + ! Set aerosol optical properties where aerosol present. Use aerosol type array as the mask. select case (itype) ! dust @@ -561,7 +561,7 @@ end subroutine compute_all_from_table ! ! Function for linearly interpolating MERRA aerosol optics tables in the rh dimension for ! a single parameter, aerosol type, spectral band, and size bin. Interpolation is performed - ! only where aerosol in present using aerosol type as the mask. + ! only where aerosol in present using aerosol type as the mask. ! function linear_interp_aero_table(table, index1, index2, weight) result(value) !$acc routine seq @@ -574,7 +574,7 @@ function linear_interp_aero_table(table, index1, index2, weight) result(value) real(wp) :: value value = table(index1) + weight * (table(index2) - table(index1)) - + end function linear_interp_aero_table ! ---------------------------------------------------------- logical function any_int_vals_outside_2D(array, checkMin, checkMax) diff --git a/rrtmgp-frontend/mo_cloud_optics_rrtmgp.F90 b/rrtmgp-frontend/mo_cloud_optics_rrtmgp.F90 index 07102b534..844e34267 100644 --- a/rrtmgp-frontend/mo_cloud_optics_rrtmgp.F90 +++ b/rrtmgp-frontend/mo_cloud_optics_rrtmgp.F90 @@ -21,14 +21,14 @@ module mo_cloud_optics_rrtmgp use mo_rte_kind, only: wp, wl use mo_rte_config, only: check_values, check_extents - use mo_rte_util_array_validation,& + use mo_rte_util_array_validation,& only: any_vals_less_than, any_vals_outside, extents_are use mo_optical_props, only: ty_optical_props, & ty_optical_props_arry, & ty_optical_props_1scl, & ty_optical_props_2str, & ty_optical_props_nstr - use mo_cloud_optics_rrtmgp_kernels, only: & + use mo_cloud_optics_rrtmgp_kernels, only: & compute_cld_from_table, compute_cld_from_pade implicit none private diff --git a/rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 b/rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 index 967d656a9..c4b020fb1 100644 --- a/rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 +++ b/rrtmgp-frontend/mo_gas_optics_rrtmgp.F90 @@ -10,7 +10,7 @@ ! BSD 3-clause license, see http://opensource.org/licenses/BSD-3-Clause ! ------------------------------------------------------------------------------------------------- ! -!> ## Class implementing the RRTMGP correlated-_k_ distribution +!> ## Class implementing the RRTMGP correlated-_k_ distribution !> !> Implements a class for computing spectrally-resolved gas optical properties and source functions !> given atmopsheric physical properties (profiles of temperature, pressure, and gas concentrations) @@ -26,11 +26,11 @@ module mo_gas_optics_rrtmgp use mo_rte_kind, only: wp, wl use mo_rte_config, only: check_extents, check_values use mo_rte_util_array, only: zero_array - use mo_rte_util_array_validation, & + use mo_rte_util_array_validation, & only: any_vals_less_than, any_vals_outside, extents_are use mo_optical_props, only: ty_optical_props use mo_source_functions, only: ty_source_func_lw - use mo_gas_optics_rrtmgp_kernels, & + use mo_gas_optics_rrtmgp_kernels, & only: interpolation, compute_tau_absorption, compute_tau_rayleigh, compute_Planck_source use mo_gas_optics_constants, only: avogad, m_dry, m_h2o, grav use mo_gas_optics_util_string, only: lower_case, string_in_array, string_loc_in_array @@ -238,7 +238,7 @@ function gas_optics_int(this, & intent(inout) :: optical_props !! Optical properties class(ty_source_func_lw ), & intent(inout) :: sources !! Planck sources - character(len=128) :: error_msg !! Empty if succssful + character(len=128) :: error_msg !! Empty if succssful ! Optional inputs real(wp), dimension(:,:), intent(in ), & optional, target :: col_dry, & !! Column dry amount; dim(ncol,nlay) @@ -311,7 +311,7 @@ function gas_optics_int(this, & ! 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 + if(present(tlev)) then error_msg = source(this, & ncol, nlay, nband, ngpt, optical_props%top_is_at_1(), & play, plev, tlay, tsfc, & @@ -320,14 +320,14 @@ function gas_optics_int(this, & tlev) !$acc exit data delete(tlev) !$omp target exit data map(release:tlev) - else + else error_msg = source(this, & ncol, nlay, nband, ngpt, optical_props%top_is_at_1(), & play, plev, tlay, tsfc, & jtemp, jpress, jeta, tropo, fmajor, & sources) - end if + end if !$acc exit data delete(tsfc) !$omp target exit data map(release:tsfc) !$acc exit data delete(jtemp, jpress, tropo, fmajor, jeta) @@ -337,7 +337,7 @@ end function gas_optics_int !------------------------------------------------------------------------------------------ ! !> Compute gas optical depth given temperature, pressure, and composition - !> Top-of-atmosphere stellar insolation is also reported + !> Top-of-atmosphere stellar insolation is also reported ! function gas_optics_ext(this, & play, plev, tlay, gas_desc, & ! mandatory inputs @@ -351,7 +351,7 @@ function gas_optics_ext(this, & type(ty_gas_concs), intent(in ) :: gas_desc !! Gas volume mixing ratios ! output class(ty_optical_props_arry), & - intent(inout) :: optical_props + intent(inout) :: optical_props real(wp), dimension(:,:), intent( out) :: toa_src !! Incoming solar irradiance(ncol,ngpt) character(len=128) :: error_msg !! Empty if successful @@ -775,7 +775,7 @@ function set_solar_variability(this, & real(wp), intent(in) :: mg_index !! facular brightening index (NRLSSI2 facular "Bremen" index) real(wp), intent(in) :: sb_index !! sunspot dimming index (NRLSSI2 sunspot "SPOT67" index) real(wp), optional, intent(in) :: tsi !! total solar irradiance - character(len=128) :: error_msg !! Empty if successful + character(len=128) :: error_msg !! Empty if successful ! ---------------------------------------------------------- integer :: igpt real(wp), parameter :: a_offset = 0.1495954_wp @@ -808,10 +808,10 @@ function set_tsi(this, tsi) result(error_msg) ! class(ty_gas_optics_rrtmgp), intent(inout) :: this real(wp), intent(in ) :: tsi !! user-specified total solar irradiance; - character(len=128) :: error_msg !! Empty if successful + character(len=128) :: error_msg !! Empty if successful real(wp) :: norm - integer :: igpt, length + integer :: igpt, length ! ---------------------------------------------------------- error_msg = "" if(tsi < 0._wp) then @@ -880,7 +880,7 @@ function source(this, & ! 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 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) & @@ -896,7 +896,7 @@ function source(this, & ! Interpolation and extrapolation at boundaries is weighted by pressure ! !$acc parallel loop gang vector - !$omp target teams distribute parallel do simd + !$omp target teams distribute parallel do simd do icol = 1, ncol tlev_arr(icol,1) = tlay(icol,1) & + (plev(icol,1)-play(icol,1))*(tlay(icol,2)-tlay(icol,1)) & @@ -905,7 +905,7 @@ function source(this, & + (plev(icol,nlay+1)-play(icol,nlay))*(tlay(icol,nlay)-tlay(icol,nlay-1)) & / (play(icol,nlay)-play(icol,nlay-1)) end do - !$acc parallel loop gang vector collapse(2) + !$acc parallel loop gang vector collapse(2) !$omp target teams distribute parallel do simd collapse(2) do ilay = 2, nlay do icol = 1, ncol @@ -1477,7 +1477,7 @@ end function get_press_max ! pure function get_temp_min(this) class(ty_gas_optics_rrtmgp), intent(in) :: this - real(wp) :: get_temp_min !! minimum temperature for which the k-dsitribution is valid + real(wp) :: get_temp_min !! minimum temperature for which the k-dsitribution is valid get_temp_min = this%temp_ref_min end function get_temp_min @@ -1556,25 +1556,25 @@ function compute_optimal_angles(this, optical_props, optimal_angles) result(err_ ! input class(ty_gas_optics_rrtmgp), intent(in ) :: this class(ty_optical_props_arry), intent(in ) :: optical_props !! Optical properties - real(wp), dimension(:,:), intent( out) :: optimal_angles !! Secant of optical transport angle + real(wp), dimension(:,:), intent( out) :: optimal_angles !! Secant of optical transport angle character(len=128) :: err_msg !! Empty if successful !---------------------------- integer :: ncol, nlay, ngpt integer :: icol, ilay, igpt, bnd real(wp) :: t, trans_total #if defined _CRAYFTN && _RELEASE_MAJOR == 14 && _RELEASE_MINOR == 0 && _RELEASE_PATCHLEVEL == 3 -# define CRAY_WORKAROUND -#endif -#ifdef CRAY_WORKAROUND - integer, allocatable :: bands(:) -#else +# define CRAY_WORKAROUND +#endif +#ifdef CRAY_WORKAROUND + integer, allocatable :: bands(:) +#else integer :: bands(optical_props%get_ngpt()) #endif !---------------------------- ncol = optical_props%get_ncol() nlay = optical_props%get_nlay() ngpt = optical_props%get_ngpt() -#ifdef CRAY_WORKAROUND +#ifdef CRAY_WORKAROUND allocate( bands(ngpt) ) ! In order to work with CCE 14 (it is also better software) #endif @@ -2000,7 +2000,7 @@ end subroutine create_gpoint_flavor !-------------------------------------------------------------------------------------------------------------------- ! ! Utility function to combine optical depths from gas absorption and Rayleigh scattering - ! It may be more efficient to combine scattering and absorption optical depths in place + ! It may be more efficient to combine scattering and absorption optical depths in place ! rather than storing and processing two large arrays ! subroutine combine_abs_and_rayleigh(tau, tau_rayleigh, optical_props) diff --git a/rrtmgp-kernels/accel/mo_gas_optics_rrtmgp_kernels.F90 b/rrtmgp-kernels/accel/mo_gas_optics_rrtmgp_kernels.F90 index 103850776..b63b54a9e 100644 --- a/rrtmgp-kernels/accel/mo_gas_optics_rrtmgp_kernels.F90 +++ b/rrtmgp-kernels/accel/mo_gas_optics_rrtmgp_kernels.F90 @@ -604,7 +604,7 @@ subroutine compute_Planck_source( & integer :: ilay, icol, igpt, ibnd, itropo, iflav integer :: gptS, gptE real(wp), dimension(2), parameter :: one = [1._wp, 1._wp] - real(wp) :: pfrac, pfrac_m1 ! Planck fraction in this layer and the one below + real(wp) :: pfrac, pfrac_m1 ! Planck fraction in this layer and the one below real(wp) :: planck_function_1, planck_function_2 ! ----------------- @@ -636,7 +636,7 @@ subroutine compute_Planck_source( & ! 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 + if (ilay == 1) then lev_src(icol,ilay, igpt) = pfrac * planck_function_1 else itropo = merge(1,2,tropo(icol,ilay-1)) !WS moved itropo inside loop for GPU @@ -645,7 +645,7 @@ subroutine compute_Planck_source( & 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 + 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 diff --git a/rrtmgp-kernels/api/mo_cloud_optics_rrtmgp_kernels.F90 b/rrtmgp-kernels/api/mo_cloud_optics_rrtmgp_kernels.F90 index ed6959fa7..f60ac05e0 100644 --- a/rrtmgp-kernels/api/mo_cloud_optics_rrtmgp_kernels.F90 +++ b/rrtmgp-kernels/api/mo_cloud_optics_rrtmgp_kernels.F90 @@ -3,7 +3,7 @@ module mo_cloud_optics_rrtmgp_kernels implicit none private public :: compute_cld_from_table, compute_cld_from_pade - interface + interface !--------------------------------------------------------------------------- ! ! Linearly interpolate values from a lookup table with "nsteps" evenly-spaced @@ -55,4 +55,3 @@ end subroutine compute_cld_from_pade end interface !--------------------------------------------------------------------------- end module mo_cloud_optics_rrtmgp_kernels - diff --git a/rrtmgp-kernels/api/mo_gas_optics_rrtmgp_kernels.F90 b/rrtmgp-kernels/api/mo_gas_optics_rrtmgp_kernels.F90 index 9272bb16e..854a0d0d7 100644 --- a/rrtmgp-kernels/api/mo_gas_optics_rrtmgp_kernels.F90 +++ b/rrtmgp-kernels/api/mo_gas_optics_rrtmgp_kernels.F90 @@ -5,7 +5,7 @@ module mo_gas_optics_rrtmgp_kernels private public :: interpolation, compute_tau_absorption, compute_tau_rayleigh, compute_Planck_source ! ------------------------------------------------------------------------------------------------------------------ - interface + interface subroutine interpolation( & ncol,nlay,ngas,nflav,neta, npres, ntemp, & flavor, & @@ -19,13 +19,13 @@ subroutine interpolation( & integer, intent(in) :: ncol,nlay !! physical domain size integer, intent(in) :: ngas,nflav,neta,npres,ntemp - !! k-distribution table dimensions + !! 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 + !! log of pressure dimension in RRTMGP tables real(wp), dimension(ntemp), intent(in) :: temp_ref - !! temperature dimension in RRTMGP tables + !! temperature dimension in RRTMGP tables real(wp), intent(in) :: press_ref_log_delta, & temp_ref_min, temp_ref_delta, & press_ref_trop_log @@ -37,14 +37,14 @@ subroutine interpolation( & 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 + !! input column gas amount - molecules/cm^2 ! outputs integer, dimension(ncol,nlay), intent(out) :: jtemp, jpress - !! temperature and pressure interpolation indexes + !! temperature and pressure interpolation indexes logical(wl), dimension(ncol,nlay), intent(out) :: tropo - !! use lower (or upper) atmosphere tables + !! use lower (or upper) atmosphere tables integer, dimension(2, ncol,nlay,nflav), intent(out) :: jeta - !! Index for binary species interpolation + !! 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 @@ -61,11 +61,11 @@ subroutine interpolation( & !! 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 + !! Interpolation weights in pressure, eta, strat/trop real(wp), dimension(2,2, ncol,nlay,nflav), intent(out) :: fminor - !! Interpolation fraction in eta, strat/trop + !! Interpolation fraction in eta, strat/trop end subroutine interpolation - end interface + end interface ! ------------------------------------------------------------------------------------------------------------------ interface subroutine compute_tau_absorption( & @@ -99,8 +99,8 @@ subroutine 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) :: 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 @@ -109,56 +109,56 @@ subroutine compute_tau_absorption( & 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 + !! beginning and ending g-point for each band real(wp), dimension(ntemp,neta,npres+1,ngpt), intent(in) :: kmajor - !! absorption coefficient table - major gases + !! absorption coefficient table - major gases real(wp), dimension(ntemp,neta,nminorklower), intent(in) :: kminor_lower - !! absorption coefficient table - minor gases, lower atmosphere + !! absorption coefficient table - minor gases, lower atmosphere real(wp), dimension(ntemp,neta,nminorkupper), intent(in) :: kminor_upper - !! absorption coefficient table - minor gases, upper atmosphere + !! 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 + !! 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? + !! 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 + 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( 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 + 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? + !! 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() + !! 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() + !! 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 + !! 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) + !! input column gas amount (molecules/cm^2) integer, dimension(2, ncol,nlay,nflav ), intent(in) :: jeta - !! interpolation indexes in eta - computed in interpolation() + !! interpolation indexes in eta - computed in interpolation() integer, dimension( ncol,nlay ), intent(in) :: jtemp - !! interpolation indexes in temperature - computed in interpolation() + !! interpolation indexes in temperature - computed in interpolation() integer, dimension( ncol,nlay ), intent(in) :: jpress - !! interpolation indexes in pressure - computed in interpolation() + !! interpolation indexes in pressure - computed in interpolation() ! --------------------- ! output - optical depth - real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau !! aborption optional depth + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau !! aborption optional depth end subroutine compute_tau_absorption - end interface + end interface ! ------------------------------------------------------------------------------------------------------------------ interface subroutine compute_tau_rayleigh(ncol,nlay,nbnd,ngpt, & @@ -170,36 +170,36 @@ subroutine compute_tau_rayleigh(ncol,nlay,nbnd,ngpt, & 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 + !! input dimensions integer, intent(in ) :: ngas,nflav,neta,npres,ntemp - !! table dimensions - integer, dimension(2,ngpt), intent(in ) :: gpoint_flavor + !! 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 + !! 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 + !! 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() + !! interpolation weights for major gases - computed in interpolation() integer, dimension(2, ncol,nlay,nflav), intent(in ) :: jeta - !! interpolation indexes in eta - computed in interpolation() + !! interpolation indexes in eta - computed in interpolation() logical(wl), dimension(ncol,nlay), intent(in ) :: tropo - !! use upper- or lower-atmospheric tables? + !! use upper- or lower-atmospheric tables? integer, dimension(ncol,nlay), intent(in ) :: jtemp - !! interpolation indexes in temperature - computed in interpolation() + !! interpolation indexes in temperature - computed in interpolation() ! outputs real(wp), dimension(ncol,nlay,ngpt), intent(out) :: tau_rayleigh - !! Rayleigh optical depth + !! Rayleigh optical depth end subroutine compute_tau_rayleigh end interface ! ------------------------------------------------------------------------------------------------------------------ - interface + interface subroutine compute_Planck_source( & ncol, nlay, nbnd, ngpt, & nflav, neta, npres, ntemp, nPlanckTemp,& @@ -210,35 +210,35 @@ subroutine compute_Planck_source( & 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 + !! input dimensions integer, intent(in) :: nflav, neta, npres, ntemp, nPlanckTemp - !! table dimensions + !! 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 + 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() + !! interpolation weights for major gases - computed in interpolation() integer, dimension(2, ncol,nlay,nflav), intent(in) :: jeta - !! interpolation indexes in eta - computed in interpolation() + !! interpolation indexes in eta - computed in interpolation() logical(wl), dimension( ncol,nlay), intent(in) :: tropo - !! use upper- or lower-atmospheric tables? + !! use upper- or lower-atmospheric tables? integer, dimension( ncol,nlay), intent(in) :: jtemp, jpress - !! interpolation indexes in temperature and pressure - computed in interpolation() + !! 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 + 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, 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 + 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 ! ------------------------------------------------------------------------------------------------------------------ diff --git a/rrtmgp-kernels/api/rrtmgp_kernels.h b/rrtmgp-kernels/api/rrtmgp_kernels.h index 6d8976c3b..fa03904aa 100644 --- a/rrtmgp-kernels/api/rrtmgp_kernels.h +++ b/rrtmgp-kernels/api/rrtmgp_kernels.h @@ -19,7 +19,7 @@ This header files defines the C bindings for the kernels used in RRTMGP extern "C" { - /* Gas optics kernels */ + /* Gas optics kernels */ void rrtmgp_interpolation( const int& ncol, const int& nlay, const int& ngas, const int& nflav, const int& neta, @@ -121,8 +121,8 @@ extern "C" Float* lev_src, // [out] (ncol,nlay+1,ngpt) Float* sfc_src_jac // [out] (ncol, ngpt) ); - - /* Cloud optics kernels */ + + /* Cloud optics kernels */ void rrtmgp_compute_tau_rayleigh( const int& ncol, const int& nlay, const int& nband, const int& ngpt, const int& ngas, const int& nflav, const int& neta, const int& npres, const int& ntemp, @@ -144,7 +144,7 @@ extern "C" const Bool* mask, // (ncol,nlay) const Float* lwp, // (ncol,nlay) const Float* re, // (ncol,nlay) - const Float& step_size, + const Float& step_size, const Float& offset, const Float* tau_table, // (nsteps, nbnd) const Float* ssa_table, // (nsteps, nbnd) @@ -162,12 +162,12 @@ extern "C" const Float* re_bounds_ext, // (nsizes+1) const Float* re_bounds_ssa, // (nsizes+1) const Float* re_bounds_asy, // (nsizes+1) - const int& m_ext, int& n_ext, - const Float* coeffs_ext, // (nbnd,nsizes,0:m_ext+n_ext) - const int& m_ssa, int& n_ssa, - const Float* coeffs_ssa, // (nbnd,nsizes,0:m_ssa+n_ssa) - const int& m_asy, int& n_asy, - const Float* coeffs_asy, // (nbnd,nsizes,0:m_asy+n_asy) + const int& m_ext, int& n_ext, + const Float* coeffs_ext, // (nbnd,nsizes,0:m_ext+n_ext) + const int& m_ssa, int& n_ssa, + const Float* coeffs_ssa, // (nbnd,nsizes,0:m_ssa+n_ssa) + const int& m_asy, int& n_asy, + const Float* coeffs_asy, // (nbnd,nsizes,0:m_asy+n_asy) Float* tau, // (ncol,nlay,nbnd) Float* taussa, // (ncol,nlay,nbnd) Float* taussag // (ncol,nlay,nbnd) diff --git a/rrtmgp-kernels/mo_cloud_optics_rrtmgp_kernels.F90 b/rrtmgp-kernels/mo_cloud_optics_rrtmgp_kernels.F90 index a5c9e295e..ed297784d 100644 --- a/rrtmgp-kernels/mo_cloud_optics_rrtmgp_kernels.F90 +++ b/rrtmgp-kernels/mo_cloud_optics_rrtmgp_kernels.F90 @@ -1,7 +1,7 @@ ! This code is part of ! RRTM for GCM Applications - Parallel (RRTMGP) ! -! Copyright 2024-, Atmospheric and Environmental Research, +! Copyright 2024-, Atmospheric and Environmental Research, ! Trustees of Columbia University. All right reserved. ! ! Use and duplication is permitted under the terms of the diff --git a/rrtmgp-kernels/mo_gas_optics_rrtmgp_kernels.F90 b/rrtmgp-kernels/mo_gas_optics_rrtmgp_kernels.F90 index df656859b..d5608969c 100644 --- a/rrtmgp-kernels/mo_gas_optics_rrtmgp_kernels.F90 +++ b/rrtmgp-kernels/mo_gas_optics_rrtmgp_kernels.F90 @@ -14,12 +14,12 @@ !> !> ## Numeric calculations for gas optics. Absorption and Rayleigh optical depths, Planck source functions. !> -!> - Interpolation coefficients are computed, then used in subsequent routines. -!> - All applications will call compute_tau_absorption(); -!> compute_tau_rayleigh() and/or compute_Planck_source() will be called depending on the -!> configuration of the k-distribution. -!> - The details of the interpolation scheme are not particaulrly important as long as arrays including -!> tables are passed consisently between kernels. +!> - Interpolation coefficients are computed, then used in subsequent routines. +!> - All applications will call compute_tau_absorption(); +!> compute_tau_rayleigh() and/or compute_Planck_source() will be called depending on the +!> configuration of the k-distribution. +!> - The details of the interpolation scheme are not particaulrly important as long as arrays including +!> tables are passed consisently between kernels. !> ! ------------------------------------------------------------------------------------------------- @@ -46,13 +46,13 @@ subroutine interpolation( & integer, intent(in) :: ncol,nlay !! physical domain size integer, intent(in) :: ngas,nflav,neta,npres,ntemp - !! k-distribution table dimensions + !! 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 + !! log of pressure dimension in RRTMGP tables real(wp), dimension(ntemp), intent(in) :: temp_ref - !! temperature dimension in RRTMGP tables + !! temperature dimension in RRTMGP tables real(wp), intent(in) :: press_ref_log_delta, & temp_ref_min, temp_ref_delta, & press_ref_trop_log @@ -64,20 +64,20 @@ subroutine interpolation( & 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 + !! input column gas amount - molecules/cm^2 ! outputs integer, dimension(ncol,nlay), intent(out) :: jtemp, jpress - !! temperature and pressure interpolation indexes + !! temperature and pressure interpolation indexes logical(wl), dimension(ncol,nlay), intent(out) :: tropo - !! use lower (or upper) atmosphere tables + !! use lower (or upper) atmosphere tables integer, dimension(2, ncol,nlay,nflav), intent(out) :: jeta !! Index for binary species interpolation real(wp), dimension(2, ncol,nlay,nflav), intent(out) :: col_mix !! combination of major species's column amounts (first index is strat/trop) real(wp), dimension(2,2,2,ncol,nlay,nflav), intent(out) :: fmajor - !! Interpolation weights in pressure, eta, strat/trop + !! Interpolation weights in pressure, eta, strat/trop real(wp), dimension(2,2, ncol,nlay,nflav), intent(out) :: fminor - !! Interpolation fraction in eta, strat/trop + !! Interpolation fraction in eta, strat/trop ! ----------------- ! local real(wp), dimension(ncol,nlay) :: ftemp, fpress ! interpolation fraction for temperature, pressure @@ -203,8 +203,8 @@ subroutine compute_tau_absorption( & tau) bind(C, name="rrtmgp_compute_tau_absorption") ! --------------------- ! input dimensions - integer, intent(in) :: ncol,nlay,nbnd,ngpt !! array sizes - integer, intent(in) :: ngas,nflav,neta,npres,ntemp !! tables sizes + 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 @@ -213,54 +213,54 @@ subroutine compute_tau_absorption( & 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 + !! beginning and ending g-point for each band real(wp), dimension(ntemp,neta,npres+1,ngpt), intent(in) :: kmajor - !! absorption coefficient table - major gases + !! absorption coefficient table - major gases real(wp), dimension(ntemp,neta,nminorklower), intent(in) :: kminor_lower - !! absorption coefficient table - minor gases, lower atmosphere + !! absorption coefficient table - minor gases, lower atmosphere real(wp), dimension(ntemp,neta,nminorkupper), intent(in) :: kminor_upper - !! absorption coefficient table - minor gases, upper atmosphere + !! 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 + !! 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? + !! 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 + 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( 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 + 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? + !! 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() + !! 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() + !! 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 + !! 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) + !! input column gas amount (molecules/cm^2) integer, dimension(2, ncol,nlay,nflav ), intent(in) :: jeta - !! interpolation indexes in eta - computed in interpolation() + !! interpolation indexes in eta - computed in interpolation() integer, dimension( ncol,nlay ), intent(in) :: jtemp - !! interpolation indexes in temperature - computed in interpolation() + !! interpolation indexes in temperature - computed in interpolation() integer, dimension( ncol,nlay ), intent(in) :: jpress - !! interpolation indexes in pressure - computed in interpolation() + !! interpolation indexes in pressure - computed in interpolation() ! --------------------- ! output - optical depth - real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau !! aborption optional depth + real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau !! aborption optional depth ! --------------------- ! Local variables ! @@ -510,32 +510,32 @@ subroutine compute_tau_rayleigh(ncol,nlay,nbnd,ngpt, & fminor,jeta,tropo,jtemp, & tau_rayleigh) bind(C, name="rrtmgp_compute_tau_rayleigh") integer, intent(in ) :: ncol,nlay,nbnd,ngpt - !! input dimensions + !! input dimensions integer, intent(in ) :: ngas,nflav,neta,npres,ntemp - !! table dimensions - integer, dimension(2,ngpt), intent(in ) :: gpoint_flavor + !! 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 + !! 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 + !! 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() + !! interpolation weights for major gases - computed in interpolation() integer, dimension(2, ncol,nlay,nflav), intent(in ) :: jeta - !! interpolation indexes in eta - computed in interpolation() + !! interpolation indexes in eta - computed in interpolation() logical(wl), dimension(ncol,nlay), intent(in ) :: tropo - !! use upper- or lower-atmospheric tables? + !! use upper- or lower-atmospheric tables? integer, dimension(ncol,nlay), intent(in ) :: jtemp - !! interpolation indexes in temperature - computed in interpolation() + !! interpolation indexes in temperature - computed in interpolation() ! outputs real(wp), dimension(ncol,nlay,ngpt), intent(out) :: tau_rayleigh - !! Rayleigh optical depth + !! Rayleigh optical depth ! ----------------- ! local variables real(wp) :: k(ngpt) ! rayleigh scattering coefficient @@ -571,35 +571,35 @@ subroutine compute_Planck_source( & 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") integer, intent(in) :: ncol, nlay, nbnd, ngpt - !! input dimensions + !! input dimensions integer, intent(in) :: nflav, neta, npres, ntemp, nPlanckTemp - !! table dimensions + !! 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 + 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() + !! interpolation weights for major gases - computed in interpolation() integer, dimension(2, ncol,nlay,nflav), intent(in) :: jeta - !! interpolation indexes in eta - computed in interpolation() + !! interpolation indexes in eta - computed in interpolation() logical(wl), dimension( ncol,nlay), intent(in) :: tropo - !! use upper- or lower-atmospheric tables? + !! use upper- or lower-atmospheric tables? integer, dimension( ncol,nlay), intent(in) :: jtemp, jpress - !! interpolation indexes in temperature and pressure - computed in interpolation() + !! 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), intent(in) :: temp_ref_min, totplnk_delta !! interpolation constants real(wp), dimension(ntemp,neta,npres+1,ngpt), intent(in) :: pfracin !! Fraction of the Planck function in each g-point - real(wp), dimension(nPlanckTemp,nbnd), intent(in) :: totplnk !! Total Planck function by band at each temperature + 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 emission from the surface + real(wp), dimension(ncol, ngpt), intent(out) :: sfc_src !! Planck emission from the surface real(wp), dimension(ncol,nlay, ngpt), intent(out) :: lay_src !! Planck emission from layer centers real(wp), dimension(ncol,nlay+1,ngpt), intent(out) :: lev_src !! Planck emission from layer boundaries - real(wp), dimension(ncol, ngpt), intent(out) :: sfc_source_Jac - !! Jacobian (derivative) of the surface Planck source with respect to surface temperature + real(wp), dimension(ncol, ngpt), intent(out) :: sfc_source_Jac + !! Jacobian (derivative) of the surface Planck source with respect to surface temperature ! ----------------- ! local real(wp), parameter :: delta_Tsurf = 1.0_wp @@ -694,8 +694,8 @@ subroutine compute_Planck_source( & end do do ilay = 2, nlay do icol = 1, ncol - lev_src(icol,ilay,igpt) = sqrt(pfrac(icol,ilay-1, igpt) * & - pfrac(icol,ilay, igpt)) & + lev_src(icol,ilay,igpt) = sqrt(pfrac(icol,ilay-1, igpt) * & + pfrac(icol,ilay, igpt)) & * planck_function(icol,ilay, ibnd) end do end do diff --git a/rte-frontend/mo_optical_props.F90 b/rte-frontend/mo_optical_props.F90 index 5f216a52e..f061c9bbc 100644 --- a/rte-frontend/mo_optical_props.F90 +++ b/rte-frontend/mo_optical_props.F90 @@ -33,7 +33,7 @@ !> !> These classes must be allocated before use. Initialization and allocation can be combined. !> The classes have a validate() function that checks all arrays for valid values (e.g. tau > 0.) -!> The vertical orientation can be specified via this%set_top_at_1() or obtained via this%top_at_1(). +!> The vertical orientation can be specified via this%set_top_at_1() or obtained via this%top_at_1(). !> !> Optical properties can be delta-scaled (though this is currently implemented only for two-stream arrays) !> @@ -54,7 +54,7 @@ module mo_optical_props use mo_rte_kind, only: wp, wl use mo_rte_config, only: check_extents, check_values - use mo_rte_util_array_validation, & + use mo_rte_util_array_validation, & only: any_vals_less_than, any_vals_outside, extents_are use mo_optical_props_kernels, only: & increment_1scalar_by_1scalar, increment_1scalar_by_2stream, increment_1scalar_by_nstream, & @@ -110,7 +110,7 @@ module mo_optical_props ! ------------------------------------------------------------------------------------------------- type, extends(ty_optical_props), abstract, public :: ty_optical_props_arry real(wp), dimension(:,:,:), allocatable :: tau !! optical depth (ncol, nlay, ngpt) - logical(wl), private :: top_at_1 ! No default - maybe uninitialized values will get caught? + logical(wl), private :: top_at_1 ! No default - maybe uninitialized values will get caught? contains procedure, public :: get_ncol procedure, public :: get_nlay @@ -120,7 +120,7 @@ module mo_optical_props procedure, public :: increment !> !> - !> + !> procedure, public :: top_is_at_1 procedure, public :: set_top_at_1 @@ -1074,8 +1074,8 @@ end function get_nmom ! ------------------------------------------------------------------------------------------ pure function top_is_at_1(this) ! - ! Vertical orientation - .true. if array index 1 is top of atmosphere - ! + ! Vertical orientation - .true. if array index 1 is top of atmosphere + ! class(ty_optical_props_arry), intent(in) :: this logical :: top_is_at_1 @@ -1084,8 +1084,8 @@ end function top_is_at_1 ! ------------------------------------------------------------------------------------------ subroutine set_top_at_1(this, top_at_1) ! - !> Set vertical orientation of class - .true. if array index 1 is top of atmosphere - ! + !> Set vertical orientation of class - .true. if array index 1 is top of atmosphere + ! class(ty_optical_props_arry), intent(inout) :: this logical, intent(in ) :: top_at_1 diff --git a/rte-frontend/mo_rte_lw.F90 b/rte-frontend/mo_rte_lw.F90 index d03eaf1e9..fdb9eed6e 100644 --- a/rte-frontend/mo_rte_lw.F90 +++ b/rte-frontend/mo_rte_lw.F90 @@ -48,7 +48,7 @@ module mo_rte_lw use mo_rte_kind, only: wp, wl use mo_rte_config, only: check_extents, check_values use mo_rte_util_array,only: zero_array - use mo_rte_util_array_validation, & + use mo_rte_util_array_validation, & only: any_vals_less_than, any_vals_outside, extents_are use mo_optical_props, only: ty_optical_props, & ty_optical_props_arry, ty_optical_props_1scl, ty_optical_props_2str, ty_optical_props_nstr @@ -129,11 +129,11 @@ function rte_lw(optical_props, & integer, parameter :: max_gauss_pts = 4 real(wp), parameter, & dimension(max_gauss_pts, max_gauss_pts) :: & - ! + ! ! Values provided are for mu = cos(theta); we require the inverse ! - gauss_Ds = 1._wp / & - RESHAPE([0.6096748751_wp, huge(1._wp) , huge(1._wp) , huge(1._wp), & + gauss_Ds = 1._wp / & + RESHAPE([0.6096748751_wp, huge(1._wp) , huge(1._wp) , huge(1._wp), & 0.2509907356_wp, 0.7908473988_wp, huge(1._wp) , huge(1._wp), & 0.1024922169_wp, 0.4417960320_wp, 0.8633751621_wp, huge(1._wp), & 0.0454586727_wp, 0.2322334416_wp, 0.5740198775_wp, 0.9030775973_wp], & @@ -351,7 +351,7 @@ function rte_lw(optical_props, & end do end if call lw_solver_noscat(ncol, nlay, ngpt, & - logical(optical_props%top_is_at_1(), wl), & + logical(optical_props%top_is_at_1(), wl), & n_quad_angs, secants, gauss_wts(1:n_quad_angs,n_quad_angs), & optical_props%tau, & sources%lay_source, & @@ -371,7 +371,7 @@ function rte_lw(optical_props, & ! ! two-stream calculation with scattering ! - call lw_solver_2stream(ncol, nlay, ngpt, & + call lw_solver_2stream(ncol, nlay, ngpt, & logical(optical_props%top_is_at_1(), wl), & optical_props%tau, optical_props%ssa, optical_props%g, & sources%lay_source, sources%lev_source, & @@ -395,7 +395,7 @@ function rte_lw(optical_props, & ! Re-scaled solution to account for scattering ! call lw_solver_noscat(ncol, nlay, ngpt, & - logical(optical_props%top_is_at_1(), wl), & + logical(optical_props%top_is_at_1(), wl), & n_quad_angs, & secants, gauss_wts(1:n_quad_angs,n_quad_angs), & optical_props%tau, & diff --git a/rte-frontend/mo_rte_sw.F90 b/rte-frontend/mo_rte_sw.F90 index 7a4b359c0..b271bd8e8 100644 --- a/rte-frontend/mo_rte_sw.F90 +++ b/rte-frontend/mo_rte_sw.F90 @@ -36,7 +36,7 @@ module mo_rte_sw use mo_rte_kind, only: wp, wl use mo_rte_config, only: check_extents, check_values use mo_rte_util_array,only: zero_array - use mo_rte_util_array_validation, & + use mo_rte_util_array_validation, & only: any_vals_less_than, any_vals_outside, extents_are use mo_optical_props, only: ty_optical_props, & ty_optical_props_arry, ty_optical_props_1scl, ty_optical_props_2str, ty_optical_props_nstr diff --git a/rte-frontend/mo_rte_util_array_validation.F90 b/rte-frontend/mo_rte_util_array_validation.F90 index 79584efe3..33ea9fb63 100644 --- a/rte-frontend/mo_rte_util_array_validation.F90 +++ b/rte-frontend/mo_rte_util_array_validation.F90 @@ -368,5 +368,5 @@ function extents_are_2d_int(array, n1, n2) extents_are_2d_int = (size(array,1) == n1 .and. & size(array,2) == n2 ) end function extents_are_2d_int - + end module mo_rte_util_array_validation diff --git a/rte-kernels/accel/mo_optical_props_kernels.F90 b/rte-kernels/accel/mo_optical_props_kernels.F90 index ccdbb531f..508668b48 100644 --- a/rte-kernels/accel/mo_optical_props_kernels.F90 +++ b/rte-kernels/accel/mo_optical_props_kernels.F90 @@ -142,13 +142,13 @@ subroutine increment_1scalar_by_1scalar(ncol, nlay, ngpt, & ! tau1 and tau2 might be the same array, thus we need to perform copy and copyin ! in separate steps. Otherwise, at the time of copyout of tau1 runtime may see - ! it as present (as present counter of tau2, having the same memory address) + ! it as present (as present counter of tau2, having the same memory address) ! is not necessarily decrement yet, and copyout action is not carried out ! (present_or_copyout semantic) !$acc data copy(tau1) !$acc data copyin(tau2) - !$acc parallel loop collapse(3) + !$acc parallel loop collapse(3) !$omp target teams distribute parallel do simd collapse(3) & !$omp& map(to:tau2) & !$omp& map(tofrom:tau1) @@ -178,7 +178,7 @@ subroutine increment_1scalar_by_2stream(ncol, nlay, ngpt, & ! tau1 and tau2 might be the same array, thus we need to perform copy and copyin ! in separate steps. Otherwise, at the time of copyout of tau1 runtime may see - ! it as present (as present counter of tau2, having the same memory address) + ! it as present (as present counter of tau2, having the same memory address) ! is not necessarily decrement yet, and copyout action is not carried out ! (present_or_copyout semantic) !$acc data copy(tau1) @@ -215,7 +215,7 @@ subroutine increment_1scalar_by_nstream(ncol, nlay, ngpt, & ! -------------- ! tau1 and tau2 might be the same array, thus we need to perform copy and copyin ! in separate steps. Otherwise, at the time of copyout of tau1 runtime may see - ! it as present (as present counter of tau2, having the same memory address) + ! it as present (as present counter of tau2, having the same memory address) ! is not necessarily decrement yet, and copyout action is not carried out ! (present_or_copyout semantic) !$acc data copy(tau1) @@ -254,7 +254,7 @@ subroutine increment_2stream_by_1scalar(ncol, nlay, ngpt, & ! -------------- ! tau1 and tau2 might be the same array, thus we need to perform copy and copyin ! in separate steps. Otherwise, at the time of copyout of tau1 runtime may see - ! it as present (as present counter of tau2, having the same memory address) + ! it as present (as present counter of tau2, having the same memory address) ! is not necessarily decrement yet, and copyout action is not carried out ! (present_or_copyout semantic) !$acc data copy(tau1, ssa1) @@ -296,7 +296,7 @@ subroutine increment_2stream_by_2stream(ncol, nlay, ngpt, & ! -------------- ! tau1 and tau2 might be the same array, thus we need to perform copy and copyin ! in separate steps. Otherwise, at the time of copyout of tau1 runtime may see - ! it as present (as present counter of tau2, having the same memory address) + ! it as present (as present counter of tau2, having the same memory address) ! is not necessarily decrement yet, and copyout action is not carried out ! (present_or_copyout semantic) !$acc data copy(tau1, ssa1, g1) @@ -348,7 +348,7 @@ subroutine increment_2stream_by_nstream(ncol, nlay, ngpt, nmom2, & ! -------------- ! tau1 and tau2 might be the same array, thus we need to perform copy and copyin ! in separate steps. Otherwise, at the time of copyout of tau1 runtime may see - ! it as present (as present counter of tau2, having the same memory address) + ! it as present (as present counter of tau2, having the same memory address) ! is not necessarily decrement yet, and copyout action is not carried out ! (present_or_copyout semantic) !$acc data copy(tau1, ssa1, g1) @@ -399,7 +399,7 @@ subroutine increment_nstream_by_1scalar(ncol, nlay, ngpt, & ! -------------- ! tau1 and tau2 might be the same array, thus we need to perform copy and copyin ! in separate steps. Otherwise, at the time of copyout of tau1 runtime may see - ! it as present (as present counter of tau2, having the same memory address) + ! it as present (as present counter of tau2, having the same memory address) ! is not necessarily decrement yet, and copyout action is not carried out ! (present_or_copyout semantic) !$acc data copy(tau1, ssa1) @@ -445,7 +445,7 @@ subroutine increment_nstream_by_2stream(ncol, nlay, ngpt, nmom1, & ! -------------- ! tau1 and tau2 might be the same array, thus we need to perform copy and copyin ! in separate steps. Otherwise, at the time of copyout of tau1 runtime may see - ! it as present (as present counter of tau2, having the same memory address) + ! it as present (as present counter of tau2, having the same memory address) ! is not necessarily decrement yet, and copyout action is not carried out ! (present_or_copyout semantic) !$acc data copy(tau1, ssa1, p1) @@ -506,7 +506,7 @@ subroutine increment_nstream_by_nstream(ncol, nlay, ngpt, nmom1, nmom2, & mom_lim = min(nmom1, nmom2) ! tau1 and tau2 might be the same array, thus we need to perform copy and copyin ! in separate steps. Otherwise, at the time of copyout of tau1 runtime may see - ! it as present (as present counter of tau2, having the same memory address) + ! it as present (as present counter of tau2, having the same memory address) ! is not necessarily decrement yet, and copyout action is not carried out ! (present_or_copyout semantic) !$acc data copy(tau1, ssa1, p1) @@ -561,7 +561,7 @@ subroutine inc_1scalar_by_1scalar_bybnd(ncol, nlay, ngpt, & ! tau1 and tau2 might be the same array, thus we need to perform copy and copyin ! in separate steps. Otherwise, at the time of copyout of tau1 runtime may see - ! it as present (as present counter of tau2, having the same memory address) + ! it as present (as present counter of tau2, having the same memory address) ! is not necessarily decrement yet, and copyout action is not carried out ! (present_or_copyout semantic) !$acc data copy(tau1) @@ -601,7 +601,7 @@ subroutine inc_1scalar_by_2stream_bybnd(ncol, nlay, ngpt, & ! tau1 and tau2 might be the same array, thus we need to perform copy and copyin ! in separate steps. Otherwise, at the time of copyout of tau1 runtime may see - ! it as present (as present counter of tau2, having the same memory address) + ! it as present (as present counter of tau2, having the same memory address) ! is not necessarily decrement yet, and copyout action is not carried out ! (present_or_copyout semantic) !$acc data copy(tau1) @@ -641,7 +641,7 @@ subroutine inc_1scalar_by_nstream_bybnd(ncol, nlay, ngpt, & ! tau1 and tau2 might be the same array, thus we need to perform copy and copyin ! in separate steps. Otherwise, at the time of copyout of tau1 runtime may see - ! it as present (as present counter of tau2, having the same memory address) + ! it as present (as present counter of tau2, having the same memory address) ! is not necessarily decrement yet, and copyout action is not carried out ! (present_or_copyout semantic) !$acc data copy(tau1) @@ -684,7 +684,7 @@ subroutine inc_2stream_by_1scalar_bybnd(ncol, nlay, ngpt, & ! tau1 and tau2 might be the same array, thus we need to perform copy and copyin ! in separate steps. Otherwise, at the time of copyout of tau1 runtime may see - ! it as present (as present counter of tau2, having the same memory address) + ! it as present (as present counter of tau2, having the same memory address) ! is not necessarily decrement yet, and copyout action is not carried out ! (present_or_copyout semantic) !$acc data copy(tau1, ssa1) @@ -729,7 +729,7 @@ subroutine inc_2stream_by_2stream_bybnd(ncol, nlay, ngpt, & ! tau1 and tau2 might be the same array, thus we need to perform copy and copyin ! in separate steps. Otherwise, at the time of copyout of tau1 runtime may see - ! it as present (as present counter of tau2, having the same memory address) + ! it as present (as present counter of tau2, having the same memory address) ! is not necessarily decrement yet, and copyout action is not carried out ! (present_or_copyout semantic) !$acc data copy(tau1, ssa1, g1) @@ -786,7 +786,7 @@ subroutine inc_2stream_by_nstream_bybnd(ncol, nlay, ngpt, nmom2, & ! tau1 and tau2 might be the same array, thus we need to perform copy and copyin ! in separate steps. Otherwise, at the time of copyout of tau1 runtime may see - ! it as present (as present counter of tau2, having the same memory address) + ! it as present (as present counter of tau2, having the same memory address) ! is not necessarily decrement yet, and copyout action is not carried out ! (present_or_copyout semantic) !$acc data copy(tau1, ssa1, g1) @@ -841,7 +841,7 @@ subroutine inc_nstream_by_1scalar_bybnd(ncol, nlay, ngpt, & ! tau1 and tau2 might be the same array, thus we need to perform copy and copyin ! in separate steps. Otherwise, at the time of copyout of tau1 runtime may see - ! it as present (as present counter of tau2, having the same memory address) + ! it as present (as present counter of tau2, having the same memory address) ! is not necessarily decrement yet, and copyout action is not carried out ! (present_or_copyout semantic) !$acc data copy(tau1, ssa1) @@ -891,7 +891,7 @@ subroutine inc_nstream_by_2stream_bybnd(ncol, nlay, ngpt, nmom1, & ! tau1 and tau2 might be the same array, thus we need to perform copy and copyin ! in separate steps. Otherwise, at the time of copyout of tau1 runtime may see - ! it as present (as present counter of tau2, having the same memory address) + ! it as present (as present counter of tau2, having the same memory address) ! is not necessarily decrement yet, and copyout action is not carried out ! (present_or_copyout semantic) !$acc data copy(tau1, ssa1, p1) @@ -955,7 +955,7 @@ subroutine inc_nstream_by_nstream_bybnd(ncol, nlay, ngpt, nmom1, nmom2, & mom_lim = min(nmom1, nmom2) ! tau1 and tau2 might be the same array, thus we need to perform copy and copyin ! in separate steps. Otherwise, at the time of copyout of tau1 runtime may see - ! it as present (as present counter of tau2, having the same memory address) + ! it as present (as present counter of tau2, having the same memory address) ! is not necessarily decrement yet, and copyout action is not carried out ! (present_or_copyout semantic) !$acc data copy(tau1, ssa1, p1) diff --git a/rte-kernels/accel/mo_rte_solver_kernels.F90 b/rte-kernels/accel/mo_rte_solver_kernels.F90 index 36d2366ec..3b122a5c3 100644 --- a/rte-kernels/accel/mo_rte_solver_kernels.F90 +++ b/rte-kernels/accel/mo_rte_solver_kernels.F90 @@ -173,8 +173,8 @@ subroutine lw_solver_noscat_oneangle(ncol, nlay, ngpt, top_at_1, D, weight, & tau_loc(icol,ilay,igpt) = tau(icol,ilay,igpt)*D(icol,igpt) trans (icol,ilay,igpt) = exp(-tau_loc(icol,ilay,igpt)) end if - call lw_source_noscat(top_at_1, & - lay_source(icol,ilay,igpt), lev_source(icol,ilay,igpt), & + call lw_source_noscat(top_at_1, & + lay_source(icol,ilay,igpt), lev_source(icol,ilay,igpt), & lev_source(icol,ilay+1,igpt), & tau_loc (icol,ilay,igpt), trans (icol,ilay,igpt), & source_dn (icol,ilay,igpt), source_up (icol,ilay,igpt)) @@ -710,7 +710,7 @@ subroutine lw_source_noscat(top_at_1, lay_source, lev_source, levp1_source, tau, ! Down at the bottom of the layer, up at the top ! -------------------------------- real(wp), parameter :: tau_thresh = sqrt(sqrt(epsilon(tau))) - real(wp) :: fact, source_inc, source_dec + real(wp) :: fact, source_inc, source_dec ! --------------------------------------------------------------- ! ! Weighting factor. Use 3rd order series expansion when rounding error (~tau^2) @@ -727,13 +727,13 @@ subroutine lw_source_noscat(top_at_1, lay_source, lev_source, levp1_source, tau, ! source_inc = (1._wp - trans) * levp1_source + 2._wp * fact * (lay_source - levp1_source) source_dec = (1._wp - trans) * lev_source + 2._wp * fact * (lay_source - lev_source) - if (top_at_1) then + if (top_at_1) then source_dn = source_inc source_up = source_dec else source_up = source_inc source_dn = source_dec - end if + end if end subroutine lw_source_noscat ! --------------------------------------------------------------- ! @@ -1125,7 +1125,7 @@ subroutine sw_dif_and_source(ncol, nlay, ngpt, top_at_1, mu0, sfc_albedo, & (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, + ! 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 ) )) diff --git a/rte-kernels/api/mo_fluxes_broadband_kernels.F90 b/rte-kernels/api/mo_fluxes_broadband_kernels.F90 index d64c1fa29..64804322a 100644 --- a/rte-kernels/api/mo_fluxes_broadband_kernels.F90 +++ b/rte-kernels/api/mo_fluxes_broadband_kernels.F90 @@ -38,8 +38,8 @@ end subroutine sum_broadband ! ---------------------------------------------------------------------------- !> !> Spectral reduction over all points for net flux - !> Overloaded - which routine is called depends on arguments - !> + !> Overloaded - which routine is called depends on arguments + !> interface net_broadband ! ---------------------------------------------------------------------------- !> diff --git a/rte-kernels/api/mo_optical_props_kernels.F90 b/rte-kernels/api/mo_optical_props_kernels.F90 index bb7e5116f..afe5aa81d 100644 --- a/rte-kernels/api/mo_optical_props_kernels.F90 +++ b/rte-kernels/api/mo_optical_props_kernels.F90 @@ -81,7 +81,7 @@ pure subroutine increment_1scalar_by_1scalar(ncol, nlay, ngpt, & 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 + end interface ! --------------------------------- !> increase absorption optical depth with extinction optical depth (2-stream form) interface @@ -93,7 +93,7 @@ pure subroutine increment_1scalar_by_2stream(ncol, nlay, ngpt, & 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 + end interface ! --------------------------------- !> increase absorption optical depth with extinction optical depth (n-stream form) interface @@ -105,7 +105,7 @@ pure subroutine increment_1scalar_by_nstream(ncol, nlay, ngpt, & 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 + end interface ! --------------------------------- ! --------------------------------- !> increment two-stream optical properties \(\tau, \omega_0, g\) with absorption optical depth @@ -118,7 +118,7 @@ pure subroutine increment_2stream_by_1scalar(ncol, nlay, ngpt, & 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 + end interface ! --------------------------------- !> increment two-stream optical properties \(\tau, \omega_0, g\) with a second set interface @@ -130,7 +130,7 @@ pure subroutine increment_2stream_by_2stream(ncol, nlay, ngpt, & 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 + end interface ! --------------------------------- !> increment two-stream optical properties \(\tau, \omega_0, g\) with _n_-stream interface @@ -144,7 +144,7 @@ pure subroutine increment_2stream_by_nstream(ncol, nlay, ngpt, nmom2, & 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 + end interface ! --------------------------------- ! --------------------------------- !> increment _n_-stream optical properties \(\tau, \omega_0, p\) with absorption optical depth @@ -157,7 +157,7 @@ pure subroutine increment_nstream_by_1scalar(ncol, nlay, ngpt, & 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 + end interface ! --------------------------------- !> increment _n_-stream optical properties \(\tau, \omega_0, p\) with two-stream values interface @@ -171,7 +171,7 @@ pure subroutine increment_nstream_by_2stream(ncol, nlay, ngpt, 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 + end interface ! --------------------------------- !> increment one set of _n_-stream optical properties with another set interface @@ -187,7 +187,7 @@ pure subroutine increment_nstream_by_nstream(ncol, nlay, ngpt, nmom1, nmom2, & 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 + end interface ! ------------------------------------------------------------------------------------------------- ! ! Incrementing when the second set of optical properties is defined at lower spectral resolution @@ -206,7 +206,7 @@ pure subroutine inc_1scalar_by_1scalar_bybnd(ncol, nlay, ngpt, & 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 + end interface ! --------------------------------- !> increase absorption optical depth defined on g-points with extinction optical depth (2-stream form) defined on bands interface @@ -220,7 +220,7 @@ pure subroutine inc_1scalar_by_2stream_bybnd(ncol, nlay, ngpt, & 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 + end interface ! --------------------------------- !> increase absorption optical depth defined on g-points with extinction optical depth (n-stream form) defined on bands interface @@ -234,7 +234,7 @@ pure subroutine inc_1scalar_by_nstream_bybnd(ncol, nlay, ngpt, & 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 + end interface ! --------------------------------- !> increment two-stream optical properties \(\tau, \omega_0, g\) defined on g-points with absorption optical depth defined on bands interface @@ -248,7 +248,7 @@ pure subroutine inc_2stream_by_1scalar_bybnd(ncol, nlay, ngpt, & 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 + end interface ! --------------------------------- !> increment 2-stream optical properties defined on g-points with another set defined on bands interface @@ -262,7 +262,7 @@ pure subroutine inc_2stream_by_2stream_bybnd(ncol, nlay, ngpt, & 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 + end interface ! --------------------------------- !> increment 2-stream optical properties defined on g-points with _n_-stream properties set defined on bands interface @@ -278,7 +278,7 @@ pure subroutine inc_2stream_by_nstream_bybnd(ncol, nlay, ngpt, 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 + end interface ! --------------------------------- ! --------------------------------- !> increment _n_-stream optical properties defined on g-points with absorption optical depth defined on bands @@ -293,7 +293,7 @@ pure subroutine inc_nstream_by_1scalar_bybnd(ncol, nlay, ngpt, & 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 + end interface ! --------------------------------- !> increment n-stream optical properties defined on g-points with 2-stream properties set defined on bands interface @@ -309,7 +309,7 @@ pure subroutine inc_nstream_by_2stream_bybnd(ncol, nlay, ngpt, nmom1, & 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 + end interface ! --------------------------------- !> increment _n_-stream optical properties defined on g-points with a second set defined on bands interface @@ -327,7 +327,7 @@ pure subroutine inc_nstream_by_nstream_bybnd(ncol, nlay, ngpt, nmom1, 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 + end interface ! ------------------------------------------------------------------------------------------------- ! ! Subsetting, meaning extracting some portion of the 3D domain diff --git a/rte-kernels/api/rte_kernels.h b/rte-kernels/api/rte_kernels.h index 71f86b81e..8debdcbc9 100644 --- a/rte-kernels/api/rte_kernels.h +++ b/rte-kernels/api/rte_kernels.h @@ -30,9 +30,9 @@ extern "C" Float* flux_dir); // [out] (ncol,nlay+1,ngpt) void rte_sw_solver_2stream( - const int& ncol, - const int& nlay, - const int& ngpt, + 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) @@ -52,10 +52,10 @@ extern "C" Float* broadband_dir); // [out] (ncol,nlay+1) void rte_lw_solver_noscat( - const int& ncol, - const int& nlay, + const int& ncol, + const int& nlay, const int& ngpt, - const Bool& top_at_1, + const Bool& top_at_1, const int& nmus, const Float* secants, // (nmus) const Float* weights, // (nmus) @@ -83,9 +83,9 @@ extern "C" void rte_lw_solver_2stream( - const int& ncol, - const int& nlay, - const int& ngpt, + 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) @@ -103,40 +103,40 @@ extern "C" // OPTICAL PROPS - INCREMENT // void rte_increment_1scalar_by_1scalar( - const int& ncol, - const int& nlay, + 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& 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& 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& 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& ncol, + const int& nlay, const int& ngpt, Float* tau_inout, // [inout] (ncol,nlay,ngpt) Float* ssa_inout, // [inout] (ncol,nlay,ngpt) @@ -156,16 +156,16 @@ extern "C" const Float* p_in); //(nmom,ncol,nlay,ngpt) void rte_increment_nstream_by_1scalar( - const int& ncol, - const int& nlay, + 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& ncol, + const int& nlay, const int& ngpt, const int& nmom1, Float* tau_inout, // [inout] (ncol,nlay,ngpt) @@ -176,10 +176,10 @@ extern "C" const Float* g_in); // (ncol,nlay,ngpt) void rte_increment_nstream_by_nstream( - const int& ncol, - const int& nlay, + const int& ncol, + const int& nlay, const int& ngpt, - const int& nmom1, + const int& nmom1, const int& nmom2, Float* tau_inout, // [inout] (ncol,nlay,ngpt) Float* ssa_inout, // [inout] (ncol,nlay,ngpt) @@ -192,8 +192,8 @@ extern "C" // OPTICAL PROPS - INCREMENT BYBND // void rte_inc_1scalar_by_1scalar_bybnd( - const int& ncol, - const int& nlay, + const int& ncol, + const int& nlay, const int& ngpt, Float* tau_inout,// [inout] (ncol,nlay,ngpt) const Float* tau_in, // (ncol,nlay,nbnd) @@ -201,8 +201,8 @@ extern "C" const int* band_lims_gpoint); // (2,nbnd) void rte_inc_1scalar_by_2stream_bybnd( - const int& ncol, - const int& nlay, + const int& ncol, + const int& nlay, const int& ngpt, Float* tau_inout, // [inout] (ncol,nlay,ngpt) const Float* tau_in, // (ncol,nlay,nbnd) @@ -231,8 +231,8 @@ extern "C" const int* band_lims_gpoint); // (2,nbnd) void rte_inc_2stream_by_2stream_bybnd( - const int& ncol, - const int& nlay, + const int& ncol, + const int& nlay, const int& ngpt, Float* tau_inout, // [inout] (ncol,nlay,ngpt) Float* ssa_inout, // [inout] (ncol,nlay,ngpt) @@ -244,9 +244,9 @@ extern "C" 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& ncol, + const int& nlay, + const int& ngpt, const int& nmom, Float* tau_inout, // [inout] (ncol,nlay,ngpt) Float* ssa_inout, // [inout] (ncol,nlay,ngpt) @@ -258,8 +258,8 @@ extern "C" const int* band_lims_gpoint); // (2,nbnd) void rte_inc_nstream_by_1scalar_bybnd( - const int& ncol, - const int& nlay, + const int& ncol, + const int& nlay, const int& ngpt, Float* tau_inout, // [inout] (ncol,nlay,ngpt) Float* ssa_inout, // [inout] (ncol,nlay,ngpt) @@ -268,9 +268,9 @@ extern "C" 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& ncol, + const int& nlay, + const int& ngpt, const int& nmom1, Float* tau_inout, // [inout] (ncol,nlay,ngpt) Float* ssa_inout, // [inout] (ncol,nlay,ngpt) @@ -283,10 +283,10 @@ extern "C" 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& 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) diff --git a/rte-kernels/api/rte_types.h b/rte-kernels/api/rte_types.h index fc645f203..c3c244fa6 100644 --- a/rte-kernels/api/rte_types.h +++ b/rte-kernels/api/rte_types.h @@ -29,6 +29,3 @@ const Float Float_epsilon = FLT_EPSILON; using Float = double; const Float Float_epsilon = DBL_EPSILON; #endif - - - diff --git a/rte-kernels/mo_rte_solver_kernels.F90 b/rte-kernels/mo_rte_solver_kernels.F90 index 6a1156224..128dd29ec 100644 --- a/rte-kernels/mo_rte_solver_kernels.F90 +++ b/rte-kernels/mo_rte_solver_kernels.F90 @@ -625,23 +625,23 @@ subroutine lw_source_noscat(ncol, nlay, top_at_1, lay_source, lev_source, tau, t tau, & ! Optical path (tau/mu) trans ! Transmissivity (exp(-tau)) real(wp), dimension(ncol, nlay+1), intent(in) :: lev_source ! Planck source at levels (layer edges) - real(wp), dimension(ncol, nlay ), target, & + real(wp), dimension(ncol, nlay ), target, & intent(out):: source_dn, source_up ! Source function at layer edges ! Down at the bottom of the layer, up at the top ! -------------------------------- - real(wp), dimension(:,:), pointer :: source_inc, source_dec + real(wp), dimension(:,:), pointer :: source_inc, source_dec integer :: icol, ilay real(wp) :: fact real(wp), parameter :: tau_thresh = sqrt(sqrt(epsilon(tau))) ! --------------------------------------------------------------- - if (top_at_1) then - source_inc => source_dn + if (top_at_1) then + source_inc => source_dn source_dec => source_up else source_inc => source_up source_dec => source_dn - end if + end if do ilay = 1, nlay do icol = 1, ncol ! @@ -664,13 +664,13 @@ subroutine lw_source_noscat(ncol, nlay, top_at_1, lay_source, lev_source, tau, t ! ! Even better - omit the layer Planck source (not working so well) ! - if(.false.) then + if(.false.) then source_inc(icol,ilay) = (1._wp - trans(icol,ilay)) * lev_source(icol,ilay+1) + & fact * (lev_source(icol,ilay ) - lev_source(icol,ilay+1)) source_dec(icol,ilay) = (1._wp - trans(icol,ilay)) * lev_source(icol,ilay ) + & fact * (lev_source(icol,ilay+1) - lev_source(icol,ilay )) - end if - end do + end if + end do end do end subroutine lw_source_noscat ! ------------------------------------------------------------------------------------------------- @@ -1060,9 +1060,9 @@ 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 + ! Compute transmission and reflection using a nominal value but mask out later ! - mu0_s = max(min_mu0, mu0(i, lay_index)) + mu0_s = max(min_mu0, mu0(i, lay_index)) k_mu = k * mu0_s ! ! Equation 14, multiplying top and bottom by exp(-k*tau) @@ -1101,7 +1101,7 @@ pure subroutine sw_dif_and_source(ncol, nlay, top_at_1, mu0, sfc_albedo, & (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, + ! 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 ) )) @@ -1113,16 +1113,16 @@ pure subroutine sw_dif_and_source(ncol, nlay, top_at_1, mu0, sfc_albedo, & end do end do ! - ! T and R for the direct beam are computed using nominal values even when the + ! 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_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 where end subroutine sw_dif_and_source ! --------------------------------------------------------------- diff --git a/tests/Readme.md b/tests/Readme.md index 8f03c7805..4e57329d3 100644 --- a/tests/Readme.md +++ b/tests/Readme.md @@ -7,7 +7,7 @@ This directory contains testing codes used in development. The Makefile target `clear_sky_regression` is used for verification, meaning ensuring the package produces results as expected (e.g. that the rescaling approximation used for LW scattering reduces to the no-scattering -solution in the absence of clouds), and for validation, meaning +solution in the absence of clouds), and for validation, meaning comparisons among different approaches (e.g. different numbers of streams). ## Feature testing diff --git a/tests/check_equivalence.F90 b/tests/check_equivalence.F90 index 86d064232..aae500d8a 100644 --- a/tests/check_equivalence.F90 +++ b/tests/check_equivalence.F90 @@ -26,8 +26,8 @@ program rte_check_equivalence use iso_fortran_env, only : error_unit ! ! Exercise various paths through RTE+RRTMGP code that should result in the same answer - ! Some sections test e.g. initialization and finalization - ! Others compare two sets of fluxes which should be the same, e.g. with respect to vertical ordering + ! Some sections test e.g. initialization and finalization + ! Others compare two sets of fluxes which should be the same, e.g. with respect to vertical ordering ! use mo_rte_kind, only: wp use mo_optical_props, only: ty_optical_props, & @@ -83,7 +83,7 @@ program rte_check_equivalence ! Output variables ! real(wp), dimension(:,:), target, & - allocatable :: ref_flux_up, ref_flux_dn, ref_flux_dir, & + allocatable :: ref_flux_up, ref_flux_dn, ref_flux_dir, & tst_flux_up, tst_flux_dn, tst_flux_dir, & heating_rate, jFluxUp ! @@ -117,7 +117,7 @@ program rte_check_equivalence ! ! Parse command line for any file names, block size ! - failed = .false. + failed = .false. nUserArgs = command_argument_count() 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..." @@ -202,14 +202,14 @@ program rte_check_equivalence ! ! 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), & - heating_rate(ncol, nlay)) - if(is_lw) then + 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), & + heating_rate(ncol, nlay)) + if(is_lw) then allocate(jFluxUp(ncol,nlay+1)) else allocate(ref_flux_dir(ncol,nlay+1), tst_flux_dir(ncol,nlay+1)) - end if + end if ! ---------------------------------------------------------------------------- ! @@ -217,7 +217,7 @@ program rte_check_equivalence ! if(is_lw) then ! - ! initialization, finalization of optical properties + ! initialization, finalization of optical properties ! call make_optical_props_1scl(gas_optics) call atmos%finalize() @@ -225,7 +225,7 @@ program rte_check_equivalence call atmos%set_name("gas only atmosphere") print *, " Intialized atmosphere twice" ! - ! Default calculation + ! Default calculation ! fluxes%flux_up => ref_flux_up(:,:) fluxes%flux_dn => ref_flux_dn(:,:) @@ -247,7 +247,7 @@ program rte_check_equivalence print *, " Computed heating rates" ! ------------------------------------------------------- ! - ! Net fluxes + ! Net fluxes ! nullify(fluxes%flux_up) nullify(fluxes%flux_dn) @@ -256,7 +256,7 @@ program rte_check_equivalence lw_sources, & sfc_emis, & fluxes)) - if(.not. allclose(fluxes%flux_net, ref_flux_dn-ref_flux_up) ) & + if(.not. allclose(fluxes%flux_net, ref_flux_dn-ref_flux_up) ) & call stop_on_err("Net fluxes don't match when computed alone") fluxes%flux_up => tst_flux_up(:,:) fluxes%flux_dn => tst_flux_dn(:,:) @@ -264,13 +264,13 @@ program rte_check_equivalence lw_sources, & sfc_emis, & fluxes)) - if(.not. allclose(fluxes%flux_net, ref_flux_dn-ref_flux_up) ) & + if(.not. allclose(fluxes%flux_net, ref_flux_dn-ref_flux_up) ) & call report_err("Net fluxes don't match when computed in tandem") - print *, " Net fluxes" + print *, " Net fluxes" nullify(fluxes%flux_net) ! ------------------------------------------------------- ! - ! Orientation invariance + ! Orientation invariance ! call lw_clear_sky_vr if(.not. allclose(tst_flux_up, ref_flux_up, tol=4._wp) .or. & @@ -279,25 +279,25 @@ program rte_check_equivalence print *, " Vertical orientation invariance" ! ------------------------------------------------------- ! - ! Subsets of atmospheric columns + ! Subsets of atmospheric columns ! call lw_clear_sky_subset - 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) .or. & + .not. allclose(tst_flux_dn, ref_flux_dn) ) & call report_err(" Doing problem in subsets fails") print *, " Subsetting invariance" ! ------------------------------------------------------- ! - ! Incrementing + ! Incrementing ! - atmos%tau(:,:,:) = 0.5_wp * atmos%tau(:,:,:) + atmos%tau(:,:,:) = 0.5_wp * atmos%tau(:,:,:) call stop_on_err(atmos%increment(atmos)) call stop_on_err(rte_lw(atmos, & lw_sources, & sfc_emis, & fluxes)) - if(.not. allclose(tst_flux_up, ref_flux_up) .or. & - .not. allclose(tst_flux_dn, ref_flux_dn) ) & + if(.not. allclose(tst_flux_up, ref_flux_up) .or. & + .not. allclose(tst_flux_dn, ref_flux_dn) ) & call report_err(" halving/doubling fails") call increment_with_1scl(atmos) @@ -305,8 +305,8 @@ program rte_check_equivalence lw_sources, & sfc_emis, & fluxes)) - 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) .or. & + .not. allclose(tst_flux_dn, ref_flux_dn) ) & call report_err(" Incrementing with 1scl fails") call increment_with_2str(atmos) @@ -314,8 +314,8 @@ program rte_check_equivalence lw_sources, & sfc_emis, & fluxes)) - 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) .or. & + .not. allclose(tst_flux_dn, ref_flux_dn) ) & call report_err(" Incrementing with 2str fails") call increment_with_nstr(atmos) @@ -323,24 +323,24 @@ program rte_check_equivalence lw_sources, & sfc_emis, & fluxes)) - 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) .or. & + .not. allclose(tst_flux_dn, ref_flux_dn) ) & call report_err(" Incrementing with nstr fails") print *, " Incrementing" ! ------------------------------------------------------- ! - ! Computing Jacobian shouldn't change net fluxes + ! Computing Jacobian shouldn't change net fluxes ! call stop_on_err(rte_lw(atmos, & lw_sources, & sfc_emis, & fluxes, & flux_up_Jac = jFluxUp)) - if(.not. allclose(tst_flux_up, ref_flux_up) .or. & - .not. allclose(tst_flux_dn, ref_flux_dn) ) & + if(.not. allclose(tst_flux_up, ref_flux_up) .or. & + .not. allclose(tst_flux_dn, ref_flux_dn) ) & call report_err(" Computing Jacobian changes fluxes") ! - ! Increase surface temperature by 1K and recompute fluxes + ! Increase surface temperature by 1K and recompute fluxes ! call stop_on_err(gas_optics%gas_optics(p_lay, p_lev, & t_lay, sfc_t + 1._wp, & @@ -353,7 +353,7 @@ program rte_check_equivalence sfc_emis, & fluxes)) ! - ! Comparision of fluxes with increased surface T aren't expected to match + ! Comparision of fluxes with increased surface T aren't expected to match ! fluxes + their Jacobian w.r.t. surface T exactly ! if (.not. allclose(tst_flux_up, ref_flux_up + jFluxUp, tol=30._wp)) then @@ -363,10 +363,10 @@ program rte_check_equivalence print *, " Jacobian" else ! - ! Shortwave + ! Shortwave ! ! - ! initialization, finalization of optical properties + ! initialization, finalization of optical properties ! call make_optical_props_2str(gas_optics) call atmos%finalize() @@ -374,7 +374,7 @@ program rte_check_equivalence print *, " Intialized atmosphere twice" ! - ! Default calculation + ! Default calculation ! fluxes%flux_up => ref_flux_up (:,:) fluxes%flux_dn => ref_flux_dn (:,:) @@ -394,12 +394,12 @@ program rte_check_equivalence ! ------------------------------------------------------- ! - ! Orientation invariance + ! Orientation invariance ! call sw_clear_sky_vr - if(.not. allclose(tst_flux_up, ref_flux_up, tol = 4._wp) .or. & - .not. allclose(tst_flux_dn, ref_flux_dn, tol = 4._wp) .or. & - .not. allclose(tst_flux_dir,ref_flux_dir,tol = 4._wp)) & + if(.not. allclose(tst_flux_up, ref_flux_up, tol = 4._wp) .or. & + .not. allclose(tst_flux_dn, ref_flux_dn, tol = 4._wp) .or. & + .not. allclose(tst_flux_dir,ref_flux_dir,tol = 4._wp)) & call report_err(" Vertical invariance failure") print *, " Vertical orientation invariance" ! ------------------------------------------------------- @@ -415,26 +415,26 @@ program rte_check_equivalence print *, " TSI invariance" ! ------------------------------------------------------- ! - ! Incrementing - ! Threshold of 4x spacing() works in double precision + ! Incrementing + ! 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)) - atmos%tau(:,:,:) = 0.5_wp * atmos%tau(:,:,:) + atmos%tau(:,:,:) = 0.5_wp * atmos%tau(:,:,:) call stop_on_err(atmos%increment(atmos)) call stop_on_err(rte_sw(atmos, mu0, toa_flux, & sfc_alb_dir, sfc_alb_dif, & fluxes)) - if(.not. allclose(tst_flux_up, ref_flux_up, tol = 8._wp) .or. & - .not. allclose(tst_flux_dn, ref_flux_dn, tol = 12._wp) .or. & - .not. allclose(tst_flux_dir,ref_flux_dir,tol = 12._wp)) & + if(.not. allclose(tst_flux_up, ref_flux_up, tol = 8._wp) .or. & + .not. allclose(tst_flux_dn, ref_flux_dn, tol = 12._wp) .or. & + .not. allclose(tst_flux_dir,ref_flux_dir,tol = 12._wp)) & call report_err(" halving/doubling fails") ! - ! Incremement with 0 optical depth + ! Incremement with 0 optical depth ! call stop_on_err(gas_optics%gas_optics(p_lay, p_lev, & t_lay, & @@ -446,9 +446,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 = 8._wp) .or. & - .not. allclose(tst_flux_dn, ref_flux_dn, tol = 12._wp) .or. & - .not. allclose(tst_flux_dir,ref_flux_dir,tol = 12._wp)) & + if(.not. allclose(tst_flux_up, ref_flux_up, tol = 8._wp) .or. & + .not. allclose(tst_flux_dn, ref_flux_dn, tol = 12._wp) .or. & + .not. allclose(tst_flux_dir,ref_flux_dir,tol = 12._wp)) & call report_err(" Incrementing with 1scl fails") call stop_on_err(gas_optics%gas_optics(p_lay, p_lev, & @@ -461,9 +461,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 = 8._wp) .or. & - .not. allclose(tst_flux_dn, ref_flux_dn, tol = 12._wp) .or. & - .not. allclose(tst_flux_dir,ref_flux_dir,tol = 12._wp)) & + if(.not. allclose(tst_flux_up, ref_flux_up, tol = 8._wp) .or. & + .not. allclose(tst_flux_dn, ref_flux_dn, tol = 12._wp) .or. & + .not. allclose(tst_flux_dir,ref_flux_dir,tol = 12._wp)) & call report_err(" Incrementing with 2str fails") call stop_on_err(gas_optics%gas_optics(p_lay, p_lev, & @@ -476,16 +476,16 @@ program rte_check_equivalence 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 = 12._wp) .or. & - .not. allclose(tst_flux_dir,ref_flux_dir,tol = 12._wp)) & + if(.not. allclose(tst_flux_up, ref_flux_up, tol = 8._wp) .or. & + .not. allclose(tst_flux_dn, ref_flux_dn, tol = 12._wp) .or. & + .not. allclose(tst_flux_dir,ref_flux_dir,tol = 12._wp)) & call report_err(" Incrementing with nstr fails") print *, " Incrementing" - end if + end if if(failed) error stop 1 contains ! ---------------------------------------------------------------------------- - ! Longwave + ! Longwave ! ---------------------------------------------------------------------------- ! Clear-sky longwave fluxes ! Reverse orientation in the vertical, compute, un-reverse @@ -571,7 +571,7 @@ subroutine lw_clear_sky_subset end subroutine lw_clear_sky_subset ! ---------------------------------------------------------------------------- - ! Shortwave + ! Shortwave ! ---------------------------------------------------------------------------- ! Shortwave - vertically reversed ! @@ -653,17 +653,17 @@ end subroutine sw_clear_sky_tsi ! ---------------------------------------------------------------------------- logical function allclose(array1, array2, tol) real(wp), dimension(:,:), intent(in) :: array1, array2 - real(wp), optional, intent(in) :: tol - - real(wp) :: tolerance - if (present(tol)) then - tolerance = tol + real(wp), optional, intent(in) :: tol + + real(wp) :: tolerance + if (present(tol)) then + tolerance = tol else tolerance = 2._wp - end if + end if allclose = all(abs(array1-array2) <= tolerance * spacing(array1)) - end function allclose + end function allclose ! ---------------------------------------------------------------------------- subroutine report_err(error_msg) use iso_fortran_env, only : error_unit diff --git a/tests/check_variants.F90 b/tests/check_variants.F90 index 385ca04a8..68f62ebdc 100644 --- a/tests/check_variants.F90 +++ b/tests/check_variants.F90 @@ -23,14 +23,14 @@ subroutine stop_on_err(error_msg) end subroutine stop_on_err ! ---------------------------------------------------------------------------------- ! -! Exercise a range of alternative approaches to RFMIP clear sky problem -! Accuracy is assessed relative to RFMIP submission with validation-plots.py -! Serves also to exercise various code paths -! Longwave: -! omiting level temperatures, use optimal angle, use three-angle integration, -! two-stream solution; reduced-resolution gas optics -! Shortwave: -! reduced-resolution gas optics +! Exercise a range of alternative approaches to RFMIP clear sky problem +! Accuracy is assessed relative to RFMIP submission with validation-plots.py +! Serves also to exercise various code paths +! Longwave: +! omiting level temperatures, use optimal angle, use three-angle integration, +! two-stream solution; reduced-resolution gas optics +! Shortwave: +! reduced-resolution gas optics ! program rte_clear_sky_regression use mo_rte_kind, only: wp @@ -206,8 +206,8 @@ program rte_clear_sky_regression allocate(flux_up(ncol,nlay+1), flux_dn(ncol,nlay+1), flux_dir(ncol,nlay+1)) ! ---------------------------------------------------------------------------- ! - ! Create output file and site, level, layer dimensions - ! + ! Create output file and site, level, layer dimensions + ! if(nf90_create(trim(output_file), NF90_CLOBBER, ncid) /= NF90_NOERR) & call stop_on_err("write_fluxes: can't create file " // trim(output_file)) if(nf90_def_dim(ncid, "site", ncol, dimid) /= NF90_NOERR) & @@ -279,7 +279,7 @@ subroutine lw_clear_sky_default atmos, & lw_sources, & tlev = t_lev)) - call stop_on_err(rte_lw(atmos, & + call stop_on_err(rte_lw(atmos, & lw_sources, & sfc_emis, & fluxes)) diff --git a/tests/mo_gas_optics_defs_rrtmgp.F90 b/tests/mo_gas_optics_defs_rrtmgp.F90 index b99422f94..67dd913b5 100644 --- a/tests/mo_gas_optics_defs_rrtmgp.F90 +++ b/tests/mo_gas_optics_defs_rrtmgp.F90 @@ -12,8 +12,8 @@ ! 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 +! 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 ! @@ -247,4 +247,4 @@ subroutine load_and_init(kdist, filename, available_gases) ncid = nf90_close(ncid) end subroutine load_and_init !-------------------------------------------------------------------------------------------------------------------- -end module mo_gas_optics_defs \ No newline at end of file +end module mo_gas_optics_defs diff --git a/tests/mo_rcemip_profiles.F90 b/tests/mo_rcemip_profiles.F90 index e09070e6c..036a2c240 100644 --- a/tests/mo_rcemip_profiles.F90 +++ b/tests/mo_rcemip_profiles.F90 @@ -17,7 +17,7 @@ ! ------------------------------------------------------------------------------ module mo_rcemip_profiles use mo_rte_kind, only: wp, wl - use mo_rte_util_array_validation, & + use mo_rte_util_array_validation, & only: extents_are use mo_gas_concentrations, only: ty_gas_concs implicit none diff --git a/tests/mo_testing_utils.F90 b/tests/mo_testing_utils.F90 index 95dce46c0..f663b3394 100644 --- a/tests/mo_testing_utils.F90 +++ b/tests/mo_testing_utils.F90 @@ -15,11 +15,11 @@ 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, & + 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 + private public :: allclose, ops_match, check_fluxes public :: stop_on_err, report_err public :: increment_with_1scl, increment_with_2str, increment_with_nstr, vr @@ -43,56 +43,56 @@ module mo_testing_utils ! ---------------------------------------------------------------------------- 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 + real(wp), optional, intent(in) :: tol + + real(wp) :: tolerance + if (present(tol)) then + tolerance = tol else tolerance = 2._wp - end if + 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 + real(wp), optional, intent(in) :: tol + + real(wp) :: tolerance + if (present(tol)) then + tolerance = tol else tolerance = 2._wp - end if + 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 + real(wp), optional, intent(in) :: tol + + real(wp) :: tolerance + if (present(tol)) then + tolerance = tol else tolerance = 2._wp - end if + 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 + real(wp), optional, intent(in) :: tol + + real(wp) :: tolerance + if (present(tol)) then + tolerance = tol else tolerance = 2._wp - end if + end if allclose_4= all(abs(tst_array-ref_array) <= tolerance * spacing(ref_array)) end function allclose_4 @@ -103,25 +103,25 @@ end function allclose_4 ! ---------------------------------------------------------------------------- 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 + 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 + real(wp), optional, intent(in) :: tol - ops_match_2str = allclose(tst_values%tau, ref_values%tau, tol) .and. & + 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 + real(wp), optional, intent(in) :: tol - ops_match_nstr = allclose(tst_values%tau, ref_values%tau, tol) .and. & + 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 @@ -152,48 +152,48 @@ subroutine stop_on_err(error_msg) end if end subroutine stop_on_err ! ---------------------------------------------------------------------------- - subroutine check_fluxes_1pair(flux_1, flux_2, status, message) + 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: ", & + 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 if end subroutine check_fluxes_1pair ! ---------------------------------------------------------------------------- - subroutine check_fluxes_2pair(flux_1, flux_2, flux_3, flux_4, status, message) + 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: ", & + 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 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? + ! 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 + ! Local variable type(ty_optical_props_1scl) :: transparent - integer :: ncol, nlay, ngpt + integer :: ncol, nlay, ngpt ncol = atmos%get_ncol() nlay = atmos%get_nlay() ngpt = atmos%get_ngpt() @@ -202,15 +202,15 @@ subroutine increment_with_1scl(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 + 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 + ! Local variable + type(ty_optical_props_2str) :: transparent + integer :: ncol, nlay, ngpt ncol = atmos%get_ncol() nlay = atmos%get_nlay() ngpt = atmos%get_ngpt() @@ -221,16 +221,16 @@ subroutine increment_with_2str(atmos) 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 + 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 + ! Local variable + type(ty_optical_props_nstr) :: transparent integer, parameter :: nmom = 4 - integer :: ncol, nlay, ngpt + integer :: ncol, nlay, ngpt ncol = atmos%get_ncol() nlay = atmos%get_nlay() ngpt = atmos%get_ngpt() @@ -241,15 +241,15 @@ subroutine increment_with_nstr(atmos) 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 + 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, & + type(ty_source_func_lw), optional, & intent(inout) :: sources integer :: nlay @@ -266,12 +266,12 @@ subroutine vr(atmos, sources) type is (ty_optical_props_nstr) atmos%ssa(:,:,:) = atmos%ssa(:,nlay:1:-1,:) atmos%p(:,:,:,:) = atmos%p(:,:,nlay:1:-1,:) - end select + end select - if(present(sources)) then + 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 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 index ed42440e4..f49a45c7c 100644 --- a/tests/rte_lw_solver_unit_tests.F90 +++ b/tests/rte_lw_solver_unit_tests.F90 @@ -14,9 +14,9 @@ program rte_lw_solver_unit_tests ! ! Exercise various paths through RTE LW solvers - ! Tests are run on an idealized problem (radiative equilibirum) + ! 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 + ! 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, & @@ -27,29 +27,29 @@ program rte_lw_solver_unit_tests 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, & + 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 + ! 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 + ! 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 + ! Answers are invariant to ! Extracting subsets ! Vertical orientation - ! Adding transparent optical properties - ! Longwave specific tests: + ! 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: + ! Other possibilites: ! Vertical discretization? Maybe just check boundary fluxes - ! Test the application of the boundary condition? + ! Test the application of the boundary condition? real(wp), parameter :: pi = acos(-1._wp) integer, parameter :: ncol = 8, nlay = 16 @@ -57,37 +57,37 @@ program rte_lw_solver_unit_tests ! ! Longwave tests - gray radiative equilibrium ! - real(wp), parameter :: sigma = 5.670374419e-8_wp, & ! Stefan-Boltzmann constant + real(wp), parameter :: sigma = 5.670374419e-8_wp, & ! Stefan-Boltzmann constant D = 1._wp/0.6096748751_wp ! Diffusivity angle, from single-angle RRTMGP solver - real(wp), dimension( ncol), parameter :: sfc_t = [(285._wp, icol = 1,ncol/2), & + 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 + 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_optical_props_1scl) :: lw_atmos type(ty_source_func_lw) :: lw_sources - type(ty_optical_props_2str) :: sw_atmos + 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, & + ref_flux_up, ref_flux_dn, ref_flux_net, & + tst_flux_up, tst_flux_dn, tst_flux_net, & jFluxUp - logical :: passed + logical :: passed ! ------------------------------------------------------------------------------------------------------ - top_at_1 = .true. + top_at_1 = .true. ! ------------------------------------------------------------------------------------------------------ - ! + ! ! Longwave tests ! ! ------------------------------------------------------------------------------------------------------ print *, "RTE LW solver unit tests" ! - ! Gray radiative equillibrium + ! 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) @@ -112,7 +112,7 @@ program rte_lw_solver_unit_tests 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 + ! Compute only net fluxes ! nullify(fluxes%flux_up) nullify(fluxes%flux_dn) @@ -122,24 +122,24 @@ program rte_lw_solver_unit_tests 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 + ! 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, & lw_sources, sfc_emis, & fluxes)) - call check_fluxes(ref_flux_net, tst_flux_dn-tst_flux_up, & + 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 + ! 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, & + tst_flux_dn, ref_flux_dn, & passed, "Doing problem in subsets fails") ! ------------------------------------------------------- @@ -153,18 +153,18 @@ program rte_lw_solver_unit_tests lw_sources, sfc_emis, & fluxes)) ! - ! Seems like these fluxes should be bitwise identical regardless of orientation + ! 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. & + 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. + 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 + ! 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) @@ -173,7 +173,7 @@ program rte_lw_solver_unit_tests sfc_emis, & fluxes, & flux_up_Jac = jFluxUp)) - call check_fluxes(tst_flux_up, ref_flux_up, tst_flux_dn, ref_flux_dn, & + 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 @@ -185,14 +185,14 @@ program rte_lw_solver_unit_tests sfc_emis, & fluxes)) ! - ! Comparision of fluxes with increased surface T aren't expected to match + ! 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, "%")', & + 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 + ! 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) @@ -208,7 +208,7 @@ program rte_lw_solver_unit_tests sfc_emis, & fluxes, & flux_up_Jac = jFluxUp)) - call check_fluxes(tst_flux_up, ref_flux_up, tst_flux_dn, ref_flux_dn, & + 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() ! ------------------------------------------------------------------------------------ @@ -221,7 +221,7 @@ program rte_lw_solver_unit_tests sfc_emis, & fluxes, & lw_Ds = lw_Ds)) - call check_fluxes(tst_flux_up, ref_flux_up, tst_flux_dn, ref_flux_dn, & + call check_fluxes(tst_flux_up, ref_flux_up, tst_flux_dn, ref_flux_dn, & passed, "Specifying diffusivity angle D fails") ! ------------------------------------------------------------------------------------ @@ -231,18 +231,18 @@ program rte_lw_solver_unit_tests print * if(.not. passed) error stop 1 ! ------------------------------------------------------------------------------------ -contains +contains ! ------------------------------------------------------------------------------------ ! - ! Define an atmosphere in gray radiative equillibrium + ! 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 + 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 + 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 @@ -251,13 +251,13 @@ subroutine gray_rad_equil(sfc_t, total_tau, nlay, top_at_1, atmos, sources) ! ! 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]), & + 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 + ! Divide optical depth evenly among layers ! atmos%tau(1:ncol,1:nlay,1) = spread(total_tau(1:ncol)/real(nlay, wp), dim=2, ncopies=nlay) call atmos%set_top_at_1(top_at_1) @@ -276,12 +276,12 @@ subroutine gray_rad_equil(sfc_t, total_tau, nlay, top_at_1, atmos, sources) 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(:) * & + 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%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 @@ -290,12 +290,12 @@ subroutine gray_rad_equil(sfc_t, total_tau, nlay, top_at_1, atmos, sources) ! 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 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... + ! 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 @@ -304,41 +304,41 @@ function check_gray_rad_equil(sfc_T, lw_tau, top_at_1, up_flux, net_flux) logical :: check_gray_rad_equil logical :: passed - integer :: toa + integer :: toa ! ------------------------------ - check_gray_rad_equil = .true. + check_gray_rad_equil = .true. toa = merge(1, size(up_flux, 2), top_at_1) ! - ! Check top-of-atmosphere energy balance + ! Check top-of-atmosphere energy balance ! if(.not. allclose(up_flux(:,toa), & gray_rad_equil_olr(sfc_t, lw_tau), tol=8._wp)) then call report_err("OLR is not consistent with gray radiative equilibrium") check_gray_rad_equil = .false. - end if + end if ! ! Check that net fluxes are constant with height - ! Fairly relaxed threshold w.r.t. spacing() because net flux is small relative to + ! 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(:,:), & + if(.not. allclose(net_flux(:,:), & spread(net_flux(:,1), dim=2, ncopies=size(net_flux,2)), & - tol = 100._wp)) then + tol = 100._wp)) then call report_err("Net flux not constant with tau in gray radiative equilibrium") check_gray_rad_equil = .false. - end if + 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 + ! 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(:)) + gray_rad_equil_olr(:) = (2._wp * sigma * T(:)**4)/(2 + D * tau(:)) end function gray_rad_equil_olr ! ------------------------------------------------------------------------------------ ! @@ -358,7 +358,7 @@ subroutine clear_sky_subset(atmos, sources, sfc_emis, 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, & + real(wp), dimension(atmos%get_ncol()/2, & atmos%get_nlay()+1), target & :: up, dn integer :: i, colS, colE @@ -382,4 +382,4 @@ subroutine clear_sky_subset(atmos, sources, sfc_emis, flux_up, flux_dn) 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 +end program rte_lw_solver_unit_tests diff --git a/tests/rte_optic_prop_unit_tests.F90 b/tests/rte_optic_prop_unit_tests.F90 index 4fea31ff5..4a1d28cf3 100644 --- a/tests/rte_optic_prop_unit_tests.F90 +++ b/tests/rte_optic_prop_unit_tests.F90 @@ -17,7 +17,7 @@ program optical_prop_unit_tests ! 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, & + 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, & @@ -29,25 +29,25 @@ program optical_prop_unit_tests 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), 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]), & + 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))', & + 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 + ! 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 + ! 2- and n-stream optical properties ! call stop_on_err(ref_2str%alloc_2str(ncol, nlay, ref_1scl)) ref_2str%tau = ref_1scl%tau @@ -57,125 +57,125 @@ program optical_prop_unit_tests 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 + ! Henyey-Greenstein phase function do imom = 1, nmom ref_nstr%p(imom,:,:,:) = g**imom - end do + end do - passed = .true. + passed = .true. ! ---------------------------------------------------------------------------- ! - ! Incrementing with transparent (tau=0) sets of optical properties + ! Incrementing with transparent (tau=0) sets of optical properties ! ! ---------------------------------------------------------------------------- print *, " Incrementing 1scl" ! - ! Increment 1scl + ! Increment 1scl ! call make_copy_1scl call increment_with_1scl(tst_1scl) - if(.not. ops_match(tst_1scl, ref_1scl)) then + if(.not. ops_match(tst_1scl, ref_1scl)) then call report_err("1scl+1scl fails") - passed = .false. - end if + passed = .false. + end if call make_copy_1scl call increment_with_2str(tst_1scl) - if(.not. ops_match(tst_1scl, ref_1scl)) then + if(.not. ops_match(tst_1scl, ref_1scl)) then call report_err("1scl+2str fails") - passed = .false. - end if + passed = .false. + end if call make_copy_1scl call increment_with_nstr(tst_1scl) - if(.not. ops_match(tst_1scl, ref_1scl)) then + if(.not. ops_match(tst_1scl, ref_1scl)) then call report_err("1scl+nstr fails") - passed = .false. - end if + passed = .false. + end if call tst_1scl%finalize() ! ---------------------------------------------------------------------------- print *, " Incrementing 2str" ! - ! Increment 2str + ! Increment 2str ! call make_copy_2str call increment_with_1scl(tst_2str) - if(.not. ops_match(tst_2str, ref_2str)) then + if(.not. ops_match(tst_2str, ref_2str)) then call report_err("2str+1scl fails") - passed = .false. - end if + passed = .false. + end if call make_copy_2str call increment_with_2str(tst_2str) - if(.not. ops_match(tst_2str, ref_2str)) then + if(.not. ops_match(tst_2str, ref_2str)) then call report_err("2str+2str fails") - passed = .false. - end if + passed = .false. + end if call make_copy_2str call increment_with_nstr(tst_2str) - if(.not. ops_match(tst_2str, ref_2str)) then + if(.not. ops_match(tst_2str, ref_2str)) then call report_err("2str+nstr fails") - passed = .false. - end if + passed = .false. + end if call tst_2str%finalize() ! ---------------------------------------------------------------------------- print *, " Incrementing nstr" ! - ! Increment nstr + ! Increment nstr ! call make_copy_nstr call increment_with_1scl(tst_nstr) - if(.not. ops_match(tst_nstr, ref_nstr)) then + if(.not. ops_match(tst_nstr, ref_nstr)) then call report_err("nstr+1scl fails") - passed = .false. - end if + passed = .false. + end if call make_copy_nstr call increment_with_2str(tst_nstr) - if(.not. ops_match(tst_nstr, ref_nstr)) then + if(.not. ops_match(tst_nstr, ref_nstr)) then call report_err("nstr+2str fails") - passed = .false. - end if + passed = .false. + end if call make_copy_nstr call increment_with_nstr(tst_nstr) - if(.not. ops_match(tst_nstr, ref_nstr)) then + if(.not. ops_match(tst_nstr, ref_nstr)) then call report_err("nstr+nstr fails") - passed = .false. - end if + passed = .false. + end if call tst_2str%finalize() ! ---------------------------------------------------------------------------- print *, " Halving/doubling optical thickness" ! - ! Adding two media of half optical thickness to recover original values + ! Adding two media of half optical thickness to recover original values ! call make_copy_1scl - tst_1scl%tau = 0.5_wp * tst_1scl%tau + 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 + if(.not. ops_match(tst_1scl, ref_1scl)) then call report_err("1scl half/double fails") - passed = .false. - end if + passed = .false. + end if call make_copy_2str - tst_2str%tau = 0.5_wp * tst_2str%tau + 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 + if(.not. ops_match(tst_2str, ref_2str)) then call report_err("2str half/double fails") - passed = .false. - end if + passed = .false. + end if call make_copy_nstr - tst_nstr%tau = 0.5_wp * tst_nstr%tau + 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 + if(.not. ops_match(tst_nstr, ref_nstr)) then call report_err("nstr half/double fails") - passed = .false. - end if + passed = .false. + end if ! ---------------------------------------------------------------------------- print *, " Delta scaling" ! @@ -183,49 +183,49 @@ program optical_prop_unit_tests ! 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 + if(.not. ops_match(tst_2str, ref_2str)) then call report_err("2str delta-scaling with f=0 fails") - passed = .false. - end if + passed = .false. + end if ! ---------------------------------------------------------------------------- if (.not. passed) call stop_on_err("Optical props unit tests fail") print *, "Optical properties unit testing finished" print * ! ---------------------------------------------------------------------------- -contains +contains ! ---------------------------------------------------------------------------- ! - ! Make copies of the existing optical depth + ! Make copies of the existing optical depth ! - subroutine make_copy_1scl + subroutine make_copy_1scl call tst_1scl%finalize() - call stop_on_err(tst_1scl%alloc_1scl(ref_1scl%get_ncol(), ref_1scl%get_nlay(), & + 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 + end subroutine make_copy_1scl ! ---------------------------------------------------------------------------- ! - ! Make copies of the existing optical depth + ! 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(), & + 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 + tst_2str%ssa = ref_2str%ssa + tst_2str%g = ref_2str%g + end subroutine make_copy_2str ! ---------------------------------------------------------------------------- ! - ! Make copies of the existing optical depth + ! 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(), & + 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 + 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 +end program optical_prop_unit_tests diff --git a/tests/rte_sw_solver_unit_tests.F90 b/tests/rte_sw_solver_unit_tests.F90 index 610e20d3a..bc292c34e 100644 --- a/tests/rte_sw_solver_unit_tests.F90 +++ b/tests/rte_sw_solver_unit_tests.F90 @@ -15,8 +15,8 @@ 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 + ! 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, & @@ -26,7 +26,7 @@ program rte_sw_solver_unit_tests 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, & + vr, & increment_with_1scl, increment_with_2str, increment_with_nstr implicit none ! ---------------------------------------------------------------------------------- @@ -34,20 +34,20 @@ program rte_sw_solver_unit_tests ! 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 + ! 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 + ! Answers are invariant to ! Extracting subsets ! Vertical orientation - ! Adding transparent optical properties - ! Shortwave specific tests: + ! Adding transparent optical properties + ! Shortwave specific tests: ! Solutions are linear in TOA flux ! - ! Other possibilites: + ! Other possibilites: ! Vertical discretization? Maybe just check boundary fluxes - ! Test the application of the boundary condition? + ! Test the application of the boundary condition? real(wp), parameter :: pi = acos(-1._wp) integer, parameter :: ncol = 8, nlay = 16 @@ -55,11 +55,11 @@ program rte_sw_solver_unit_tests integer :: icol, ilay, imu0 ! - ! Shorteave tests - thin atmosphere + ! 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 - & + 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 @@ -67,11 +67,11 @@ program rte_sw_solver_unit_tests real(wp), parameter :: factor = 2._wp ! for checking linearity real(wp), dimension(ncol) :: mu0_arr - type(ty_optical_props_2str) :: atmos + type(ty_optical_props_2str) :: atmos type(ty_fluxes_broadband) :: fluxes - logical, parameter :: top_at_1 = .true. + logical, parameter :: top_at_1 = .true. real(wp), dimension(ncol,nlay+1), target :: & - ref_flux_up, ref_flux_dn, ref_flux_dir, ref_flux_net, & + 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: @@ -79,7 +79,7 @@ program rte_sw_solver_unit_tests #endif real(wp), dimension(:), pointer :: sfc - logical :: passed + logical :: passed ! ------------------------------------------------------------------------------------ ! @@ -89,8 +89,8 @@ program rte_sw_solver_unit_tests 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]), & + 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 atmos%set_top_at_1(top_at_1) @@ -104,7 +104,7 @@ program rte_sw_solver_unit_tests fluxes%flux_dn_dir => ref_flux_dir(:,:) fluxes%flux_net => ref_flux_net(:,:) call stop_on_err(rte_sw(atmos, & - mu0_arr, & + mu0_arr, & toa_flux, & sfc_albedo, sfc_albedo, & fluxes)) @@ -112,24 +112,24 @@ program rte_sw_solver_unit_tests ! 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, & + ! 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(atmos%top_is_at_1()) then + if(atmos%top_is_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), & + 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 + end if ! ------------------------------------------------------------------------------------ ! ! Net fluxes on- vs off-line @@ -138,37 +138,37 @@ program rte_sw_solver_unit_tests 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 + ! Compute only net fluxes ! nullify(fluxes%flux_up) nullify(fluxes%flux_dn) call stop_on_err(rte_sw(atmos, & - mu0_arr, & + 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 + ! 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, & - mu0_arr, & + mu0_arr, & toa_flux, & sfc_albedo, sfc_albedo, & fluxes)) - call check_fluxes(ref_flux_net, tst_flux_dn-tst_flux_up, & + 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 + ! 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, & + tst_flux_dn, ref_flux_dn, & passed, "Doing problem in subsets fails") ! ------------------------------------------------------- ! @@ -178,7 +178,7 @@ program rte_sw_solver_unit_tests call thin_scattering(tau, ssa, g, nlay, atmos) call vr(atmos) call stop_on_err(rte_sw(atmos, & - mu0_arr, & + mu0_arr, & toa_flux, & sfc_albedo, sfc_albedo, & fluxes)) @@ -189,8 +189,8 @@ program rte_sw_solver_unit_tests 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, & + 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) @@ -201,35 +201,35 @@ program rte_sw_solver_unit_tests print *, " Linear in TOA flux" call thin_scattering(tau, ssa, g, nlay, atmos) call stop_on_err(rte_sw(atmos, & - mu0_arr, & + 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, & + 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 + end do ! Done ! print *, "RTE SW solver unit tests done" print * if(.not. passed) error stop 1 ! ------------------------------------------------------------------------------------ -contains +contains ! ------------------------------------------------------------------------------------ ! - ! Define an atmosphere in gray radiative equillibrium + ! 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, 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)) & + real(wp), dimension(size(tau)*size(ssa)*size(g)) & :: temp ntau = size(tau); nssa = size(ssa); ng = size(g) @@ -238,18 +238,18 @@ subroutine thin_scattering(tau, ssa, g, nlay, atmos) ! ! 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]), & + 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 + ! 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 + ! ... 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) @@ -264,21 +264,21 @@ subroutine thin_scattering(tau, ssa, g, nlay, atmos) print * end if ! - ! Delta-scale + ! 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 + 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 - + real(wp), dimension(atmos%get_ncol()) :: gamma3, R, T ! Reflectance, transmittance + check_thin_scattering = .true. ! ! Solutions for small tau @@ -297,7 +297,7 @@ function check_thin_scattering(atmos, mu0, top_at_1, ref_flux_up, ref_flux_dn, 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 if end function check_thin_scattering ! ------------------------------------------------------------------------------------ ! @@ -316,7 +316,7 @@ subroutine clear_sky_subset(atmos, mu0, toa_flux, sfc_albedo, 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, & + real(wp), dimension(atmos%get_ncol()/2, & atmos%get_nlay()+1), target & :: up, dn integer :: i, colS, colE @@ -332,7 +332,7 @@ subroutine clear_sky_subset(atmos, mu0, toa_flux, sfc_albedo, flux_up, flux_dn) colE = i * ncol/2 call stop_on_err(atmos%get_subset(colS, ncol/2, atmos_subset)) call stop_on_err(rte_sw(atmos_subset, & - mu0(colS:colE), & + mu0(colS:colE), & toa_flux(colS:colE,:), & sfc_albedo(:,colS:colE), sfc_albedo(:,colS:colE), & fluxes)) @@ -342,4 +342,4 @@ subroutine clear_sky_subset(atmos, mu0, toa_flux, sfc_albedo, flux_up, flux_dn) end subroutine clear_sky_subset ! ------------------------------------------------------------------------------------ -end program rte_sw_solver_unit_tests \ No newline at end of file +end program rte_sw_solver_unit_tests diff --git a/tests/validation-plots.py b/tests/validation-plots.py index 9564e9216..01f35b2fb 100755 --- a/tests/validation-plots.py +++ b/tests/validation-plots.py @@ -73,13 +73,13 @@ def main(): ) parser.add_argument( "--lw_vars_file", - help="Path to the LW results file", + help="Path to the LW results file", default="lw_flux_variants.nc" ) parser.add_argument( "--sw_vars_file", - help="Path to the SW results file", + help="Path to the SW results file", default="sw_flux_variants.nc" ) From cc0f75b44d536fd10fa0556474926e58bfbb3f4f Mon Sep 17 00:00:00 2001 From: Robert Pincus Date: Tue, 14 Jan 2025 07:26:38 -0500 Subject: [PATCH 2/3] Updating Cmake-format version --- examples/rfmip-clear-sky/CMakeLists.txt | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/examples/rfmip-clear-sky/CMakeLists.txt b/examples/rfmip-clear-sky/CMakeLists.txt index 8026f3afc..37c69177a 100644 --- a/examples/rfmip-clear-sky/CMakeLists.txt +++ b/examples/rfmip-clear-sky/CMakeLists.txt @@ -21,8 +21,7 @@ set(inoutputs list( TRANSFORM inoutputs PREPEND ${RRTMGP_DATA}/examples/rfmip-clear-sky/inputs/ - OUTPUT_VARIABLE - inputs + OUTPUT_VARIABLE inputs ) # The tests write to the input files, therefore we copy them: From 71520f8bcc8af9611e2a0f70eb3278fcc3064343 Mon Sep 17 00:00:00 2001 From: Sergey Kosukhin Date: Tue, 14 Jan 2025 16:43:35 +0100 Subject: [PATCH 3/3] Format Python code --- examples/all-sky/make_problem_size_loop.py | 99 ++++---- examples/all-sky/run-allsky-example.py | 70 ++++-- examples/compare-to-reference.py | 95 ++++---- .../rfmip-clear-sky/run-rfmip-examples.py | 57 +++-- tests/validation-plots.py | 212 +++++++++++------- 5 files changed, 332 insertions(+), 201 deletions(-) diff --git a/examples/all-sky/make_problem_size_loop.py b/examples/all-sky/make_problem_size_loop.py index 587dd8c76..8b72e788b 100644 --- a/examples/all-sky/make_problem_size_loop.py +++ b/examples/all-sky/make_problem_size_loop.py @@ -13,52 +13,73 @@ # output name -if __name__ == '__main__': - parser = argparse.ArgumentParser( - description="Description here ") +if __name__ == "__main__": + parser = argparse.ArgumentParser(description="Description here ") # Argument parseing described at # https://stackoverflow.com/questions/15753701/how-can-i-pass-a-list-as-a-command-line-argument-with-argparse - parser.add_argument("-x", "--executable", - type=str, - default="./rrtmgp_allsky", - help="Path to exectuable") - parser.add_argument("-c", "--ncol", - type=lambda items: [int(i) for i in - list(csv.reader([items]))[0]], - default="2,4,8,16", - help="Number of columns (multiple e.g. 2,4,8,16)") - parser.add_argument("-l", "--nlay", - type=lambda items: [int(i) for i in - list(csv.reader([items]))[0]], - default="32, 64, 96", - help="Number of layers (multiple e.g. 32,64.96)") - parser.add_argument("-i", "--nloops", - type=int, default=1, - help="Number of loops (same for all ncol)") - parser.add_argument("-o", "--output_file", - type=str, - default="rrtmgp-allsky-output.nc", - help="Path to output file") - parser.add_argument("-k", "--k_distribution", - type=str, - required=True, - help="Path to gas optics file [required]") - parser.add_argument("-cl", "--cloud_optics", - type=str, default="", - help="Path to cloud optics file") - parser.add_argument("-a", "--aerosol_optics", - type=str, default="", - help="Path to aerosol optics file") + parser.add_argument( + "-x", + "--executable", + type=str, + default="./rrtmgp_allsky", + help="Path to exectuable", + ) + parser.add_argument( + "-c", + "--ncol", + type=lambda items: [int(i) for i in list(csv.reader([items]))[0]], + default="2,4,8,16", + help="Number of columns (multiple e.g. 2,4,8,16)", + ) + parser.add_argument( + "-l", + "--nlay", + type=lambda items: [int(i) for i in list(csv.reader([items]))[0]], + default="32, 64, 96", + help="Number of layers (multiple e.g. 32,64.96)", + ) + parser.add_argument( + "-i", + "--nloops", + type=int, + default=1, + help="Number of loops (same for all ncol)", + ) + parser.add_argument( + "-o", + "--output_file", + type=str, + default="rrtmgp-allsky-output.nc", + help="Path to output file", + ) + parser.add_argument( + "-k", + "--k_distribution", + type=str, + required=True, + help="Path to gas optics file [required]", + ) + parser.add_argument( + "-cl", "--cloud_optics", type=str, default="", help="Path to cloud optics file" + ) + parser.add_argument( + "-a", + "--aerosol_optics", + type=str, + default="", + help="Path to aerosol optics file", + ) args = parser.parse_args() # Can't supply aerosols without clouds if args.cloud_optics == "" and args.aerosol_optics != "": - raise AssertionError( - "Need to supply cloud optics if providing aerosol optics") + raise AssertionError("Need to supply cloud optics if providing aerosol optics") # Every combo of ncol, nlay for l in args.nlay: for i in args.ncol: - print(f"{args.executable} {i:6d} {l:4d} {args.nloops:3d} " + - f"{args.output_file} {args.k_distribution} " + - f"{args.cloud_optics} {args.aerosol_optics} ") + print( + f"{args.executable} {i:6d} {l:4d} {args.nloops:3d} " + + f"{args.output_file} {args.k_distribution} " + + f"{args.cloud_optics} {args.aerosol_optics} " + ) diff --git a/examples/all-sky/run-allsky-example.py b/examples/all-sky/run-allsky-example.py index dd9b95a25..1a911bb08 100644 --- a/examples/all-sky/run-allsky-example.py +++ b/examples/all-sky/run-allsky-example.py @@ -20,37 +20,63 @@ # In the local directory all_sky_exe_name = os.path.join(all_sky_dir, "rrtmgp_allsky") -if __name__ == '__main__': +if __name__ == "__main__": parser = argparse.ArgumentParser( - description="Runs all-sky examples, resetting output.") - parser.add_argument("--run_command", type=str, default="", - help="Prefix ('jsrun' etc.) for running commands. " - "Use quote marks to enclose multi-part commands.") - parser.add_argument("--ncol", type=int, default=24, - help="Number of columns to compute") - parser.add_argument("--nlay", type=int, default=72, - help="Number of layers " - "(every one will have the same clouds)") - parser.add_argument("--nloops", type=int, default=1, - help="Number of times to compute 'nloops' " - "cloudy columns") + description="Runs all-sky examples, resetting output." + ) + parser.add_argument( + "--run_command", + type=str, + default="", + help="Prefix ('jsrun' etc.) for running commands. " + "Use quote marks to enclose multi-part commands.", + ) + parser.add_argument( + "--ncol", type=int, default=24, help="Number of columns to compute" + ) + parser.add_argument( + "--nlay", + type=int, + default=72, + help="Number of layers " "(every one will have the same clouds)", + ) + parser.add_argument( + "--nloops", + type=int, + default=1, + help="Number of times to compute 'nloops' " "cloudy columns", + ) args = parser.parse_args() - ncol_str = '{0:5d}'.format(args.ncol) - nlay_str = '{0:5d}'.format(args.nlay) - nloops_str = '{0:5d}'.format(args.nloops) + ncol_str = "{0:5d}".format(args.ncol) + nlay_str = "{0:5d}".format(args.nlay) + nloops_str = "{0:5d}".format(args.nloops) if args.run_command: print("using the run command") all_sky_exe_name = args.run_command + " " + all_sky_exe_name os.chdir(all_sky_dir) # Remove cloudy-sky fluxes from the file containing the atmospheric profiles subprocess.run( - [all_sky_exe_name, ncol_str, nlay_str, nloops_str, - "rrtmgp-allsky-lw-no-aerosols.nc", lw_gas_coeffs_file, - lw_clouds_coeff_file]) + [ + all_sky_exe_name, + ncol_str, + nlay_str, + nloops_str, + "rrtmgp-allsky-lw-no-aerosols.nc", + lw_gas_coeffs_file, + lw_clouds_coeff_file, + ] + ) subprocess.run( - [all_sky_exe_name, ncol_str, nlay_str, nloops_str, - "rrtmgp-allsky-sw-no-aerosols.nc", sw_gas_coeffs_file, - sw_clouds_coeff_file]) + [ + all_sky_exe_name, + ncol_str, + nlay_str, + nloops_str, + "rrtmgp-allsky-sw-no-aerosols.nc", + sw_gas_coeffs_file, + sw_clouds_coeff_file, + ] + ) # end main diff --git a/examples/compare-to-reference.py b/examples/compare-to-reference.py index 6cb1ce9b0..fcd094eef 100644 --- a/examples/compare-to-reference.py +++ b/examples/compare-to-reference.py @@ -9,58 +9,74 @@ # Thresholds come from environement variables when set? # # -import sys - import argparse -import numpy as np import os +import sys import warnings + +import numpy as np import xarray as xr # # Comparing reference and test results # -if __name__ == '__main__': +if __name__ == "__main__": warnings.simplefilter("ignore", xr.SerializationWarning) parser = argparse.ArgumentParser( - description="Compares example output to file in reference directory") - parser.add_argument("--ref_dir", type=str, - help="Directory containing reference files") - parser.add_argument("--tst_dir", type=str, default=".", - help="Directory contining test values") - parser.add_argument("--file_names", type=str, nargs='+', default=[], - help="Name[s] of files to compare") - parser.add_argument("--variables", type=str, nargs='+', default=None, - help="Name[s] of files to compare") - parser.add_argument("--report_threshold", type=float, - default=os.getenv("REPORTING_THRESHOLD", 0.), - help="Threshold for reporting differences") - parser.add_argument("--failure_threshold", type=float, - default=os.getenv("FAILURE_THRESHOLD", 1.e-5), - help="Threshold at which differences cause failure " - "(for continuous integration)") + description="Compares example output to file in reference directory" + ) + parser.add_argument( + "--ref_dir", type=str, help="Directory containing reference files" + ) + parser.add_argument( + "--tst_dir", type=str, default=".", help="Directory contining test values" + ) + parser.add_argument( + "--file_names", + type=str, + nargs="+", + default=[], + help="Name[s] of files to compare", + ) + parser.add_argument( + "--variables", + type=str, + nargs="+", + default=None, + help="Name[s] of files to compare", + ) + parser.add_argument( + "--report_threshold", + type=float, + default=os.getenv("REPORTING_THRESHOLD", 0.0), + help="Threshold for reporting differences", + ) + parser.add_argument( + "--failure_threshold", + type=float, + default=os.getenv("FAILURE_THRESHOLD", 1.0e-5), + help="Threshold at which differences cause failure " + "(for continuous integration)", + ) args = parser.parse_args() tst = xr.open_mfdataset( - [os.path.join(args.tst_dir, f) for f in args.file_names], - combine='by_coords') + [os.path.join(args.tst_dir, f) for f in args.file_names], combine="by_coords" + ) ref = xr.open_mfdataset( - [os.path.join(args.ref_dir, f) for f in args.file_names], - combine='by_coords') + [os.path.join(args.ref_dir, f) for f in args.file_names], combine="by_coords" + ) variables = args.variables if args.variables is not None else tst.variables failed = False for v in variables: if np.any(np.isnan(ref.variables[v].values)): - raise Exception( - v + ": some ref values are missing. That's not right.") + raise Exception(v + ": some ref values are missing. That's not right.") if np.all(np.isnan(tst.variables[v].values)): - raise Exception( - v + ": all test values are missing. Were the tests run?") + raise Exception(v + ": all test values are missing. Were the tests run?") if np.any(np.isnan(tst.variables[v].values)): - raise Exception( - v + ": some test values are missing. Now that is strange.") + raise Exception(v + ": some test values are missing. Now that is strange.") # # Reporting # @@ -69,19 +85,20 @@ avg = 0.5 * (tst + ref)[v].values # Division raises a runtime warning when we divide by zero even if # the values in those locations will be ignored. - with np.errstate(divide='ignore', invalid='ignore'): - frac_diff = np.where( - (avg > 2. * np.finfo(float).eps), diff / avg, 0) - print('Variable %s differs (max abs difference: %e; ' - 'max percent difference: %e%%)' % ( - v, diff.max(), 100.0 * frac_diff.max())) + with np.errstate(divide="ignore", invalid="ignore"): + frac_diff = np.where((avg > 2.0 * np.finfo(float).eps), diff / avg, 0) + print( + "Variable %s differs (max abs difference: %e; " + "max percent difference: %e%%)" + % (v, diff.max(), 100.0 * frac_diff.max()) + ) else: - print('Variable %s: No diffs' % v) + print("Variable %s: No diffs" % v) # # Failure # - if not np.allclose(tst[v], ref[v], atol=args.failure_threshold, - rtol=0): failed = True + if not np.allclose(tst[v], ref[v], atol=args.failure_threshold, rtol=0): + failed = True if failed: print("Tests failed") diff --git a/examples/rfmip-clear-sky/run-rfmip-examples.py b/examples/rfmip-clear-sky/run-rfmip-examples.py index c747cefd1..3520317c1 100755 --- a/examples/rfmip-clear-sky/run-rfmip-examples.py +++ b/examples/rfmip-clear-sky/run-rfmip-examples.py @@ -15,39 +15,54 @@ # Code should be run in the rfmip_dir directory conds_file = os.path.join( - os.environ["RRTMGP_DATA"], "examples", "rfmip-clear-sky", "inputs", - "multiple_input4MIPs_radiation_RFMIP_UColorado-RFMIP-1-2_none.nc") -lw_gas_coeffs_file = os.path.join( - os.environ["RRTMGP_DATA"], "rrtmgp-gas-lw-g256.nc") -sw_gas_coeffs_file = os.path.join( - os.environ["RRTMGP_DATA"], "rrtmgp-gas-sw-g224.nc") + os.environ["RRTMGP_DATA"], + "examples", + "rfmip-clear-sky", + "inputs", + "multiple_input4MIPs_radiation_RFMIP_UColorado-RFMIP-1-2_none.nc", +) +lw_gas_coeffs_file = os.path.join(os.environ["RRTMGP_DATA"], "rrtmgp-gas-lw-g256.nc") +sw_gas_coeffs_file = os.path.join(os.environ["RRTMGP_DATA"], "rrtmgp-gas-sw-g224.nc") rfmip_lw_exe_name = "./rrtmgp_rfmip_lw" rfmip_sw_exe_name = "./rrtmgp_rfmip_sw" -if __name__ == '__main__': +if __name__ == "__main__": parser = argparse.ArgumentParser( - description="Runs all-sky examples, resetting output.") - parser.add_argument("--run_command", type=str, default="", - help="Prefix ('jsrun' etc.) for running commands. Use " - "quote marks to enclose multi-part commands.") - parser.add_argument("--block_size", type=int, default=8, - help="Number of columns to compute at a time. Must be " - "a factor of 1800 (ncol*nexp)") + description="Runs all-sky examples, resetting output." + ) + parser.add_argument( + "--run_command", + type=str, + default="", + help="Prefix ('jsrun' etc.) for running commands. Use " + "quote marks to enclose multi-part commands.", + ) + parser.add_argument( + "--block_size", + type=int, + default=8, + help="Number of columns to compute at a time. Must be " + "a factor of 1800 (ncol*nexp)", + ) args = parser.parse_args() - block_size_str = '{0:4d}'.format(args.block_size) + block_size_str = "{0:4d}".format(args.block_size) if args.run_command: print("using the run command") rfmip_lw_exe_name = args.run_command + " " + rfmip_lw_exe_name rfmip_sw_exe_name = args.run_command + " " + rfmip_sw_exe_name print( - "Running " + rfmip_lw_exe_name + block_size_str + " " + - conds_file + " " + lw_gas_coeffs_file) + "Running " + + rfmip_lw_exe_name + + block_size_str + + " " + + conds_file + + " " + + lw_gas_coeffs_file + ) # arguments are block size, input conditions, coefficient files, forcing # index, physics index - subprocess.run( - [rfmip_lw_exe_name, block_size_str, conds_file, lw_gas_coeffs_file]) - subprocess.run( - [rfmip_sw_exe_name, block_size_str, conds_file, sw_gas_coeffs_file]) + subprocess.run([rfmip_lw_exe_name, block_size_str, conds_file, lw_gas_coeffs_file]) + subprocess.run([rfmip_sw_exe_name, block_size_str, conds_file, sw_gas_coeffs_file]) diff --git a/tests/validation-plots.py b/tests/validation-plots.py index 01f35b2fb..3a9fd40a8 100755 --- a/tests/validation-plots.py +++ b/tests/validation-plots.py @@ -1,15 +1,16 @@ #! /usr/bin/env python +import argparse +import urllib.request +import warnings + import colorcet as cc import matplotlib as mpl import matplotlib.pyplot as plt import numpy as np -import urllib.request -import warnings import xarray as xr from matplotlib.backends.backend_pdf import PdfPages -import argparse def mae(diff, col_dim): # @@ -25,23 +26,17 @@ def rms(diff, col_dim): return np.sqrt(np.square(diff).mean(dim=col_dim)) -def make_comparison_plot(variants, labels, reference, vscale, colors, - col_dim="site"): +def make_comparison_plot(variants, labels, reference, vscale, colors, col_dim="site"): # # Make a plot comparing differences with respect to reference # if type(variants) is not list: - make_comparison_plot([variants], [labels], reference, vscale, colors, - col_dim) + make_comparison_plot([variants], [labels], reference, vscale, colors, col_dim) else: for i in np.arange(0, len(variants)): delta = variants[i] - reference - plt.plot(mae(delta, col_dim), - vscale, '-', - color=colors[i], label=labels[i]) - plt.plot(rms(delta, col_dim), - vscale, '--', - color=colors[i]) + plt.plot(mae(delta, col_dim), vscale, "-", color=colors[i], label=labels[i]) + plt.plot(rms(delta, col_dim), vscale, "--", color=colors[i]) # Reverse vertical ordering plt.legend() # Reverse vertical ordering @@ -54,13 +49,25 @@ def construct_lbl_esgf_root(var, esgf_node="llnl"): # line-by-line results # model = "LBLRTM-12-8" - prefix = ("http://esgf3.dkrz.de/thredds/fileServer/cmip6/RFMIP/AER/" + - model + "/rad-irf/r1i1p1f1/Efx/") + prefix = ( + "http://esgf3.dkrz.de/thredds/fileServer/cmip6/RFMIP/AER/" + + model + + "/rad-irf/r1i1p1f1/Efx/" + ) if esgf_node == "llnl": - prefix = ("http://esgf-data1.llnl.gov/thredds/fileServer/css03_data/" - "CMIP6/RFMIP/AER/" + model + "/rad-irf/r1i1p1f1/Efx/") - return (prefix + var + "/gn/v20190514/" + var + - "_Efx_" + model + "_rad-irf_r1i1p1f1_gn.nc") + prefix = ( + "http://esgf-data1.llnl.gov/thredds/fileServer/css03_data/" + "CMIP6/RFMIP/AER/" + model + "/rad-irf/r1i1p1f1/Efx/" + ) + return ( + prefix + + var + + "/gn/v20190514/" + + var + + "_Efx_" + + model + + "_rad-irf_r1i1p1f1_gn.nc" + ) ######################################################################## @@ -69,24 +76,22 @@ def main(): parser.add_argument( "--state_file", help="Path to the state information NetCDF file.", - default="multiple_input4MIPs_radiation_RFMIP_UColorado-RFMIP-1-2_none.nc" + default="multiple_input4MIPs_radiation_RFMIP_UColorado-RFMIP-1-2_none.nc", ) parser.add_argument( "--lw_vars_file", help="Path to the LW results file", - default="lw_flux_variants.nc" - + default="lw_flux_variants.nc", ) parser.add_argument( "--sw_vars_file", help="Path to the SW results file", - default="sw_flux_variants.nc" - + default="sw_flux_variants.nc", ) parser.add_argument( "--output_pdf", help="Path to the output PDF file for validation plots.", - default="validation-figures.pdf" + default="validation-figures.pdf", ) args = parser.parse_args() @@ -105,17 +110,17 @@ def main(): for v in fluxes: try: try: - urllib.request.urlretrieve(construct_lbl_esgf_root(v), - v + lbl_suffix) + urllib.request.urlretrieve(construct_lbl_esgf_root(v), v + lbl_suffix) except: urllib.request.urlretrieve( - construct_lbl_esgf_root(v, esgf_node="dkrz"), - v + lbl_suffix) + construct_lbl_esgf_root(v, esgf_node="dkrz"), v + lbl_suffix + ) except: raise Exception("Failed to download {0}".format(v + lbl_suffix)) - lbl = xr.open_mfdataset([v + lbl_suffix for v in fluxes], - combine="by_coords").sel(expt=0) + lbl = xr.open_mfdataset([v + lbl_suffix for v in fluxes], combine="by_coords").sel( + expt=0 + ) # # Open the test results @@ -125,9 +130,10 @@ def main(): # Does the flux plus the Jacobian equal a calculation with perturbed surface # temperature? # - gp['lw_flux_up_from_deriv'] = gp.lw_flux_up + gp.lw_jaco_up + gp["lw_flux_up_from_deriv"] = gp.lw_flux_up + gp.lw_jaco_up gp.lw_flux_up_from_deriv.attrs = { - "description": "LW flux up, surface T+1K, computed from Jacobian"} + "description": "LW flux up, surface T+1K, computed from Jacobian" + } ######################################################################## # # The RFMIP cases are on an irregular pressure grid so we can't compute @@ -137,16 +143,28 @@ def main(): # plevs = np.arange(0, lbl.plev.max(), 1000) plevs[0] = lbl.plev.min().values - plev = xr.DataArray(plevs, dims=['plev'], coords={'plev': plevs}) - lbli = xr.concat([lbl.sel(site=i).reset_coords().swap_dims( - {"level": "plev"}).interp(plev=plev) for i in - np.arange(0, lbl.site.size)], dim='site') + plev = xr.DataArray(plevs, dims=["plev"], coords={"plev": plevs}) + lbli = xr.concat( + [ + lbl.sel(site=i) + .reset_coords() + .swap_dims({"level": "plev"}) + .interp(plev=plev) + for i in np.arange(0, lbl.site.size) + ], + dim="site", + ) gpi = xr.concat( - [gp.sel(site=i).rename( - {"pres_level": "plev"}).reset_coords().swap_dims( - {"level": "plev"}).interp(plev=plev) - for i in np.arange(0, gp.site.size)], - dim='site') + [ + gp.sel(site=i) + .rename({"pres_level": "plev"}) + .reset_coords() + .swap_dims({"level": "plev"}) + .interp(plev=plev) + for i in np.arange(0, gp.site.size) + ], + dim="site", + ) cols = cc.glasbey_dark mpl.rcParams["legend.frameon"] = False @@ -160,31 +178,57 @@ def main(): # # Accuracy - 3-angle and single-angle # - variants = [[gpi.lw_flux_dn, gpi.lw_flux_dn_alt, gpi.lw_flux_dn_optang, - gpi.lw_flux_dn_alt_oa, gpi.lw_flux_dn_3ang, - gpi.lw_flux_dn_2str, gpi.lw_flux_dn_1rescl], - [gpi.lw_flux_up, gpi.lw_flux_up_alt, gpi.lw_flux_up_optang, - gpi.lw_flux_up_alt_oa, gpi.lw_flux_up_3ang, - gpi.lw_flux_up_2str, gpi.lw_flux_up_1rescl], - [gpi.lw_flux_net, - gpi.lw_flux_net_alt, - gpi.lw_flux_dn_optang - gpi.lw_flux_up_optang, - gpi.lw_flux_dn_alt_oa - gpi.lw_flux_up_alt_oa, - gpi.lw_flux_dn_3ang - gpi.lw_flux_up_3ang, - gpi.lw_flux_dn_2str - gpi.lw_flux_up_2str, - gpi.lw_flux_dn_1rescl - gpi.lw_flux_up_1rescl]] + variants = [ + [ + gpi.lw_flux_dn, + gpi.lw_flux_dn_alt, + gpi.lw_flux_dn_optang, + gpi.lw_flux_dn_alt_oa, + gpi.lw_flux_dn_3ang, + gpi.lw_flux_dn_2str, + gpi.lw_flux_dn_1rescl, + ], + [ + gpi.lw_flux_up, + gpi.lw_flux_up_alt, + gpi.lw_flux_up_optang, + gpi.lw_flux_up_alt_oa, + gpi.lw_flux_up_3ang, + gpi.lw_flux_up_2str, + gpi.lw_flux_up_1rescl, + ], + [ + gpi.lw_flux_net, + gpi.lw_flux_net_alt, + gpi.lw_flux_dn_optang - gpi.lw_flux_up_optang, + gpi.lw_flux_dn_alt_oa - gpi.lw_flux_up_alt_oa, + gpi.lw_flux_dn_3ang - gpi.lw_flux_up_3ang, + gpi.lw_flux_dn_2str - gpi.lw_flux_up_2str, + gpi.lw_flux_dn_1rescl - gpi.lw_flux_up_1rescl, + ], + ] refs = [lbli.rld, lbli.rlu, lbli.rld - lbli.rlu] - titles = ["Accuracy wrt LBLRTM: LW down", "Accuracy wrt LBLRTM: LW up", - "Accuracy: LW net"] + titles = [ + "Accuracy wrt LBLRTM: LW down", + "Accuracy wrt LBLRTM: LW up", + "Accuracy: LW net", + ] for v, r, t in zip(variants, refs, titles): - make_comparison_plot(v, - labels=["default", "fewer-g-points", - "optimal-angle", - "fewer points + optimal-angle", - "3-angle", "2-stream", "rescaled"], - reference=r, - vscale=plev / 100., - colors=cols) + make_comparison_plot( + v, + labels=[ + "default", + "fewer-g-points", + "optimal-angle", + "fewer points + optimal-angle", + "3-angle", + "2-stream", + "rescaled", + ], + reference=r, + vscale=plev / 100.0, + colors=cols, + ) plt.ylabel("Pressure (Pa)") plt.xlabel("Error (W/m$^2$), solid=mean, dash=RMS") plt.title(t) @@ -198,16 +242,20 @@ def main(): # Using scattering representations of optical properties to do a # no-scattering calculation # - variants = [[gpi.lw_flux_dn_notlev, gpi.lw_flux_dn_2str], - [gpi.lw_flux_up_notlev, gpi.lw_flux_up_2str]] + variants = [ + [gpi.lw_flux_dn_notlev, gpi.lw_flux_dn_2str], + [gpi.lw_flux_up_notlev, gpi.lw_flux_up_2str], + ] refs = [gpi.lw_flux_dn, gpi.lw_flux_up] titles = ["Variants: LW down", "Variants: LW up"] for v, r, t in zip(variants, refs, titles): - make_comparison_plot(v, - labels=["no-tlev", "2str"], - reference=r, - vscale=plev / 100., - colors=cols) + make_comparison_plot( + v, + labels=["no-tlev", "2str"], + reference=r, + vscale=plev / 100.0, + colors=cols, + ) plt.ylabel("Pressure (Pa)") plt.xlabel("Difference (W/m$^2$), solid=mean, dash=RMS") plt.title(t) @@ -221,16 +269,20 @@ def main(): # # Accuracy comparison # - variants = [[gpi.sw_flux_dn, gpi.sw_flux_dn_alt], - [gpi.sw_flux_up, gpi.sw_flux_up_alt]] + variants = [ + [gpi.sw_flux_dn, gpi.sw_flux_dn_alt], + [gpi.sw_flux_up, gpi.sw_flux_up_alt], + ] refs = [lbli.rsd, lbli.rsu] titles = ["Accuracy: SW down", "Accuracy: SW up"] for v, r, t in zip(variants, refs, titles): - make_comparison_plot(v, - labels=["default", "fewer-g-points"], - reference=r, - vscale=plev / 100., - colors=cols) + make_comparison_plot( + v, + labels=["default", "fewer-g-points"], + reference=r, + vscale=plev / 100.0, + colors=cols, + ) plt.ylabel("Pressure (Pa)") plt.xlabel("Error (W/m$^2$), solid=mean, dash=RMS") plt.title(t) @@ -238,5 +290,5 @@ def main(): plt.close() -if __name__ == '__main__': +if __name__ == "__main__": main()