diff --git a/.github/workflows/Dockerfile.gnu b/.github/workflows/Dockerfile.gnu
new file mode 100644
index 0000000000..3506c2b9ee
--- /dev/null
+++ b/.github/workflows/Dockerfile.gnu
@@ -0,0 +1,68 @@
+#***********************************************************************
+#* GNU Lesser General Public License
+#*
+#* This file is part of the GFDL Flexible Modeling System (FMS).
+#*
+#* FMS is free software: you can redistribute it and/or modify it under
+#* the terms of the GNU Lesser General Public License as published by
+#* the Free Software Foundation, either version 3 of the License, or (at
+#* your option) any later version.
+#*
+#* FMS is distributed in the hope that it will be useful, but WITHOUT
+#* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+#* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+#* for more details.
+#*
+#* You should have received a copy of the GNU Lesser General Public
+#* License along with FMS. If not, see .
+#***********************************************************************
+# FMS CI image recipefile for GNU
+# Runs on centos stream (builder has same base from redhat registry)
+#
+# arguments to specify versions to build can be given to docker or changed here (--build-arg name=val)
+FROM spack/rockylinux9:latest as builder
+
+ARG gcc_version=12.3.0
+ARG netcdfc_version=4.9.0
+ARG netcdff_version=4.6.0
+ARG libyaml_version=0.2.5
+ARG mpich_version=4.0.2
+
+COPY spack.env /opt/deps/spack.env
+
+# perl's download kept timing out
+RUN sed -i 's/connect_timeout: 10/connect_timeout: 600/' /opt/spack/etc/spack/defaults/config.yaml && \
+ spack install gcc@${gcc_version} && \
+ source /opt/spack/share/spack/setup-env.sh && \
+ spack load gcc@${gcc_version} && \
+ spack compiler find && \
+ sed "s/COMPILER/gcc@$gcc_version/" /opt/deps/spack.env > spack.yaml && \
+ sed -i "s/NETCDF_C_VERSION/$netcdfc_version/" spack.yaml && \
+ sed -i "s/NETCDF_F_VERSION/$netcdff_version/" spack.yaml && \
+ sed -i "s/LIBYAML_VERSION/$libyaml_version/" spack.yaml && \
+ sed -i "s/MPI_LIB/mpich@$mpich_version/" spack.yaml && \
+ spack env activate -d . && \
+ spack -e . concretize -f > /opt/deps/deps.log && \
+ spack install --fail-fast
+
+# copy built software to base from first image
+FROM rockylinux:9
+
+COPY --from=builder /opt/view/ /opt/view/
+COPY --from=builder /opt/deps/ /opt/deps/
+
+# input files used with --enable-input-tests
+# need to be on the dev boxes if building
+COPY ./fms_test_input /home/unit_tests_input
+
+RUN dnf install -y autoconf make automake m4 libtool pkg-config zip
+
+ENV FC="mpifort"
+ENV CC="mpicc"
+ENV MPICH_FC="/opt/view/bin/gfortran"
+ENV MPICH_CC="/opt/view/bin/gcc"
+ENV FCFLAGS="-I/opt/view/include"
+ENV CFLAGS="-I/opt/view/include"
+ENV LDFLAGS="-L/opt/view/lib"
+ENV LD_LIBRARY_PATH="/opt/view/lib:/opt/view/lib64:/usr/local/lib:/usr/local/lib64"
+ENV PATH="/opt/view/bin:/usr/local/bin:/usr/bin:/usr/local/sbin:/usr/sbin"
diff --git a/.github/workflows/github_autotools_gnu.yml b/.github/workflows/github_autotools_gnu.yml
new file mode 100644
index 0000000000..b7008a0814
--- /dev/null
+++ b/.github/workflows/github_autotools_gnu.yml
@@ -0,0 +1,39 @@
+# 'main' required ci, does a distcheck (builds, tests, check install)
+# image created off dockerfile in repo, compile/link flags are set there
+name: Build libFMS test with autotools
+
+on: [push, pull_request]
+
+jobs:
+ build:
+ runs-on: ubuntu-latest
+ strategy:
+ matrix:
+ conf-flag: [ --disable-openmp, --enable-mixed-mode, --disable-setting-flags, --with-mpi=no]
+ input-flag: [--with-yaml, --enable-test-input=/home/unit_tests_input]
+ exclude:
+ - conf-flag: --with-mpi=no
+ input-flag: --enable-test-input=/home/unit_tests_input
+ container:
+ image: noaagfdl/fms-ci-rocky-gnu:12.3.0
+ env:
+ TEST_VERBOSE: 1
+ DISTCHECK_CONFIGURE_FLAGS: "${{ matrix.conf-flag }} ${{ matrix.input-flag }} ${{ matrix.io-flag }}"
+ SKIP_TESTS: "test_yaml_parser.5" # temporary till fixes are in
+ steps:
+ - name: Checkout code
+ uses: actions/checkout@v2
+ - name: Prepare GNU autoconf for build
+ run: autoreconf -if
+ - name: Configure the build
+ if: ${{ matrix.conf-flag != '--disable-setting-flags' }}
+ run: ./configure ${DISTCHECK_CONFIGURE_FLAGS} || cat config.log
+ - name: Configure the build with compiler flags
+ if: ${{ matrix.conf-flag == '--disable-setting-flags' }}
+ run: ./configure ${DISTCHECK_CONFIGURE_FLAGS} FCFLAGS="-fdefault-real-8 -fdefault-double-8 -fcray-pointer -ffree-line-length-none -I/usr/include $FCFLAGS" || cat config.log
+ - name: Build the library
+ run: make distcheck
+ if: ${{ matrix.conf-flag != '--with-mpi=no' }}
+ - name: Build the library (without test suite for serial build)
+ run: make
+ if: ${{ matrix.conf-flag == '--with-mpi=no' }}
diff --git a/.github/workflows/intel_pr.yml b/.github/workflows/github_autotools_intel.yml
similarity index 92%
rename from .github/workflows/intel_pr.yml
rename to .github/workflows/github_autotools_intel.yml
index 62a15361ea..851ad71520 100644
--- a/.github/workflows/intel_pr.yml
+++ b/.github/workflows/github_autotools_intel.yml
@@ -2,13 +2,16 @@ on: pull_request
jobs:
intel-autotools:
runs-on: ubuntu-latest
+ strategy:
+ matrix:
+ io-flag: ["--disable-deprecated-io", "--enable-deprecated-io"]
container:
image: intel/oneapi-hpckit:2023.1.0-devel-ubuntu20.04
env:
CC: mpiicc
FC: mpiifort
CFLAGS: "-I/libs/include"
- FCFLAGS: "-I/libs/include -g -traceback"
+ FCFLAGS: "-I/libs/include -g -traceback ${{ matrix.io-flag }}"
LDFLAGS: "-L/libs/lib"
TEST_VERBOSE: 1
I_MPI_FABRICS: "shm" # needed for mpi in image
diff --git a/.github/workflows/build_cmake_gnu.yml b/.github/workflows/github_cmake_gnu.yml
similarity index 100%
rename from .github/workflows/build_cmake_gnu.yml
rename to .github/workflows/github_cmake_gnu.yml
diff --git a/.github/workflows/coupler.yml b/.github/workflows/github_coupler_gnu.yml
similarity index 100%
rename from .github/workflows/coupler.yml
rename to .github/workflows/github_coupler_gnu.yml
diff --git a/.github/workflows/update_docs.yml b/.github/workflows/github_doc_site.yml
similarity index 100%
rename from .github/workflows/update_docs.yml
rename to .github/workflows/github_doc_site.yml
diff --git a/.github/workflows/lint_fms.yml b/.github/workflows/github_linter.yml
similarity index 100%
rename from .github/workflows/lint_fms.yml
rename to .github/workflows/github_linter.yml
diff --git a/.github/workflows/am4_regression_parallelWorks_intel_tag.yml b/.github/workflows/parallelworks_am4_intel.yml
similarity index 100%
rename from .github/workflows/am4_regression_parallelWorks_intel_tag.yml
rename to .github/workflows/parallelworks_am4_intel.yml
diff --git a/.github/workflows/spack.env b/.github/workflows/spack.env
new file mode 100644
index 0000000000..69a3bdcbd0
--- /dev/null
+++ b/.github/workflows/spack.env
@@ -0,0 +1,17 @@
+# template for spack environment yaml
+# uppercase words get replaced before activating
+spack:
+ specs:
+ - COMPILER
+ - MPI_LIB
+ - netcdf-c@NETCDF_C_VERSION ^MPI_LIB
+ - netcdf-fortran@NETCDF_F_VERSION
+ - libyaml@LIBYAML_VERSION
+ concretizer:
+ unify: true
+ packages:
+ all:
+ compiler: [ COMPILER ]
+ config:
+ install_tree: /opt/deps
+ view: /opt/view
diff --git a/.github/workflows/version.yml b/.github/workflows/version.yml
index 1e71236386..31641aea69 100644
--- a/.github/workflows/version.yml
+++ b/.github/workflows/version.yml
@@ -1,3 +1,5 @@
+# appends -dev to the version upon release and opens pr
+# CI won't run on generated PR, easiest workaround is to close + reopen
on:
release:
types: [published]
@@ -16,4 +18,4 @@ jobs:
branch-suffix: timestamp # add a timestamp to branch name
delete-branch: true # delete afer merge
title: Append dev to version number post-release
- body: automated change, adds '-dev' to the version number upon releases
+ body: automated change, adds '-dev' to the version number upon releases. This PR will need to be closed and reopened to run CI testing.
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 0cc9802f8f..2c616e647d 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -6,6 +6,45 @@ and this project uses `yyyy.rr[.pp]`, where `yyyy` is the year a patch is releas
`rr` is a sequential release number (starting from `01`), and an optional two-digit
sequential patch number (starting from `01`).
+## [2023.02] - 2023-07-27
+### Known Issues
+- GCC 11.1.0 is unsupported due to compilation issues with select type. The issue is resolved in later GCC releases.
+- When outputting sub-region diagnostics, the current diag_manager does not add "tileX" to the filename when using a cube sphere. This leads to trouble when trying to combine the files and regrid them (if the region is in two different tiles)
+- GCC 10 and greater causing io issues when compiled using O2 optimization flags
+- GNU compilers prior to the GCC 9.0 release are unsupported for this release due to lack of support for the findloc intrinsic function. This will result in an error saying 'findloc' has no IMPLICIT type and can be resolved by compiling with gcc version 9.0 or greater.
+
+### Added
+- MPP/EXCHANGE: Adds association checks before pointer deallocations in mpp includes and xgrid
+
+### Changed
+- LIBFMS: The libFMS.F90 file (module name `fms`) meant to provide global access has been updated to include 'fms' and it's module/subdirectory name as prefixes for all names. This will only affect external codes that are already using the global module (via `use fms`) and not individual modules.
+- MIXED PRECISION: Updates the axis_utils2, horiz_interp, sat_vapor_pressure, and axis_utils subdirectories to support mixed precision real values.
+- FMS2_IO: Added in mpp_scatter and mpp_gather performance changes from the 2023.01.01 patch. See below for more details.
+- FMS2_IO: Improved error messages to give more debugging information
+- FMS_MOD: Changed fms_init to include a system call to set the stack size to unlimited, removed previously added stack size fixes
+- MONIN_OBUKHOV: Restructures the subroutines in `stable_mix` interface so that 1d calls the underlying implementation, and 2 and 3d call it on 1d slices of the data as opposed to passing in mismatched arrays.
+- MPP: Updates from JEDI for ajoint version the mpp halo filling (mpp_do_update_ad.fh), adds checkpoint for forward buffer information.
+
+### Fixed
+- MPP: mpp_broadcast causing an unintended error message due to checking the wrong pe value
+- MPP: Added workaround for GCC 12 issues causing errors with string lengths in fms2_io
+- FMS2_IO: Fixed support for 'packed' data when using NF_SHORT variables. Scale_factor and add_offset attributes will now be applied if present.
+- DOCS: Improved doxygen comments for tranlon, updated deployment action for site
+- TESTS: Workaround added for ICE coming from mpp_alltoall test with intel 2022.3, and fixes for any test scripts missing input.nml creation. Fixes for mpp/test_global_array failures.
+- TIME_INTERP: Fixes crashes when calling with a non-existant field
+- DIAG_MANAGER: Fixes a module dependency issue causing failures during parallel builds
+- AXIS_UTILS2: Fixes an out of bounds memory index
+
+### Removed
+- FMS_IO/MPP_IO: The two older io modules, fms_io_mod and mpp_io_mod, have been deprecated and will not be compiled by default. If you wish to compile these modules, you must use the -Duse_deprecated_io CPP flag or the --enable-deprecated-io configure option if building with autotools.
+
+### Tag Commit Hashes
+- 2023.02-beta1 2be8aa452ad3e5f43e92c38a64f12d1ae6c43fb8
+- 2023.02-alpha3 8c73bd18dc1d580f2ee524c37cf903ff54d40501
+- 2023.02-alpha2 783019fdec89a8db2b26247c2f63d4782e1495c0
+- 2023.02-alpga1 419c66be31f82ebb13a91ea5e837c707eb54473b
+
+
## [2023.01.01] - 2023-06-06
### Changed
- FMS2_IO: Performance changes for domain_reads_2d and domain_reads_3d:
diff --git a/CI.md b/CI.md
index 225b25129b..89b4db256e 100644
--- a/CI.md
+++ b/CI.md
@@ -8,24 +8,30 @@ Required CI for pull requests are listed first.
## Pull Request CI and checks
### Build libFMS with autotools
+
Required GNU build test for all pull requests/pushes.
Runs `make distcheck` after configuring via GNU autotools.
+Runs on a container image with spack installed dependencies, on top a rocky linux base.
+
+Dockerfile for image is stored at .github/workflows/Dockerfile.gnu for more specific information on the CI environment.
+
Container environment:
-gcc v7.3.0
-mpich v3.3a2
-netcdf v4.6.0
-netcdf-fortran v4.4.4
+gcc v12.3.0
+mpich v4.0.2
+netcdf v4.9.0
+netcdf-fortran v4.6.0
autoconf v2.69
+libyaml v0.2.5
-container hosted at [noaagfdl/ubuntu_libfms_gnu:latest](https://hub.docker.com/r/noaagfdl/ubuntu_libfms_gnu)
-
-`./configure` flags:
-- `--enable-openmp`
+`./configure` flags tested:
- `--disable-openmp`
- `--enable-mixed-mode`
+- `--with-mpi=no` (disables unit testing)
- `--disable-setting-flags`
- `--with-yaml`
+- `--enable-test-input=/home/unit_tests_input`
+
### Build libfms with cmake
Required GNU build test for all pull requests/pushes.
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 61fccf9c62..338dfb8f18 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -26,7 +26,7 @@ set(CMAKE_Fortran_FLAGS_DEBUG)
# Define the CMake project
project(FMS
- VERSION 2023.01.0
+ VERSION 2023.02.0
DESCRIPTION "GFDL FMS Library"
HOMEPAGE_URL "https://www.gfdl.noaa.gov/fms"
LANGUAGES C Fortran)
@@ -196,6 +196,7 @@ list(APPEND fms_fortran_src_files
# Collect FMS C source files
list(APPEND fms_c_src_files
affinity/affinity.c
+ fms/fms_stacksize.c
mosaic/create_xgrid.c
mosaic/gradient_c2l.c
mosaic/interp.c
@@ -314,9 +315,12 @@ foreach(kind ${kinds})
constants
topography/include
axis_utils/include
+ constants
+ astronomy/include
field_manager/include
time_interp/include
- tracer_manager/include)
+ tracer_manager/include
+ interpolator/include)
target_compile_definitions(${libTgt}_f PRIVATE "${fms_defs}")
target_compile_definitions(${libTgt}_f PRIVATE "${${kind}_defs}")
@@ -353,6 +357,7 @@ foreach(kind ${kinds})
$
$
$
+ $
$
$
$
@@ -363,7 +368,10 @@ foreach(kind ${kinds})
$
$
$
- $)
+ $
+ $
+ $
+ $)
target_include_directories(${libTgt} INTERFACE
diff --git a/astronomy/Makefile.am b/astronomy/Makefile.am
index 33d4acaf3c..50a1449412 100644
--- a/astronomy/Makefile.am
+++ b/astronomy/Makefile.am
@@ -23,14 +23,23 @@
# Ed Hartnett 2/22/19
# Include .h and .mod files.
-AM_CPPFLAGS = -I$(top_srcdir)/include
+AM_CPPFLAGS = -I$(top_srcdir)/include -I$(top_srcdir)/astronomy/include
AM_FCFLAGS = $(FC_MODINC). $(FC_MODOUT)$(MODDIR)
# Build this uninstalled convenience library.
noinst_LTLIBRARIES = libastronomy.la
# The convenience library depends on its source.
-libastronomy_la_SOURCES = astronomy.F90
+libastronomy_la_SOURCES = \
+ astronomy.F90 \
+ include/astronomy_r4.fh \
+ include/astronomy_r8.fh \
+ include/astronomy.inc
+
+astronomy.$(FC_MODEXT): \
+include/astronomy_r4.fh \
+include/astronomy_r8.fh \
+include/astronomy.inc
BUILT_SOURCES = astronomy_mod.$(FC_MODEXT)
nodist_include_HEADERS = astronomy_mod.$(FC_MODEXT)
diff --git a/astronomy/astronomy.F90 b/astronomy/astronomy.F90
index 192fc8f22a..0a7af2cc6e 100644
--- a/astronomy/astronomy.F90
+++ b/astronomy/astronomy.F90
@@ -29,296 +29,320 @@
module astronomy_mod
-use fms_mod, only: fms_init, &
- mpp_pe, mpp_root_pe, stdlog, &
- write_version_number, &
- check_nml_error, error_mesg, &
- FATAL, NOTE, WARNING
-use time_manager_mod, only: time_type, set_time, get_time, &
- get_date_julian, set_date_julian, &
- set_date, length_of_year, &
- time_manager_init, &
- operator(-), operator(+), &
- operator( // ), operator(<)
-use constants_mod, only: constants_init, PI
-use mpp_mod, only: input_nml_file
-
-!--------------------------------------------------------------------
-
-implicit none
-private
-
-!---------------------------------------------------------------------
-!----------- version number for this module --------------------------
+ use fms_mod, only: fms_init, &
+ mpp_pe, mpp_root_pe, stdlog, &
+ write_version_number, &
+ check_nml_error, error_mesg, &
+ FATAL, NOTE, WARNING
+ use time_manager_mod, only: time_type, set_time, get_time, &
+ get_date_julian, set_date_julian, &
+ set_date, length_of_year, &
+ time_manager_init, &
+ operator(-), operator(+), &
+ operator( // ), operator(<)
+ use constants_mod, only: constants_init, PI
+ use mpp_mod, only: input_nml_file
+ use platform_mod, only: r4_kind, r8_kind, i4_kind, i8_kind
+
+ !--------------------------------------------------------------------
+
+ implicit none
+ private
+
+ !---------------------------------------------------------------------
+ !----------- version number for this module --------------------------
! Include variable "version" to be written to log file.
#include
-!---------------------------------------------------------------------
-!------- interfaces --------
-
-public &
- astronomy_init, get_period, set_period, &
- set_orbital_parameters, get_orbital_parameters, &
- set_ref_date_of_ae, get_ref_date_of_ae, &
- diurnal_solar, daily_mean_solar, annual_mean_solar, &
- astronomy_end, universal_time, orbital_time
-
-!> @}
-
-!> @brief Calculates solar information for the given location(lat & lon) and time
-!!
-!> ~~~~~~~~~~{.f90}
-!! call diurnal_solar (lat, lon, time, cosz, fracday, rrsun, dt_time)
-!! call diurnal_solar (lat, lon, gmt, time_since_ae, cosz, fracday, rrsun, dt)
-!! ~~~~~~~~~~
-!!
-!! The first option (used in conjunction with time_manager_mod)
-!! generates the real variables gmt and time_since_ae from the
-!! time_type input, and then calls diurnal_solar with these real inputs.
-!!
-!! The time of day is set by
-!! ~~~~~~~~~~{.f90}
-!! real, intent(in) :: gmt
-!! ~~~~~~~~~~
-!! The time of year is set by
-!! ~~~~~~~~~~{.f90}
-!! real, intent(in) :: time_since_ae
-!! ~~~~~~~~~~
-!! with time_type input, both of these are extracted from
-!! ~~~~~~~~~~{.f90}
-!! type(time_type), intent(in) :: time
-!! ~~~~~~~~~~
-!!
-!! Separate routines exist within this interface for scalar,
-!! 1D or 2D input and output fields:
-!!
-!! ~~~~~~~~~~{.f90}
-!! real, intent(in), dimension(:,:) :: lat, lon
-!! real, intent(in), dimension(:) :: lat, lon
-!! real, intent(in) :: lat, lon
-!!
-!! real, intent(out), dimension(:,:) :: cosz, fracday
-!! real, intent(out), dimension(:) :: cosz, fracday
-!! real, intent(out) :: cosz, fracday
-!! ~~~~~~~~~~
-!!
-!! One may also average the output fields over the time interval
-!! between gmt and gmt + dt by including the optional argument dt (or
-!! dt_time). dt is measured in radians and must be less than pi
-!! (1/2 day). This average is computed analytically, and should be
-!! exact except for the fact that changes in earth-sun distance over
-!! the time interval dt are ignored. In the context of a diurnal GCM,
-!! this option should always be employed to insure that the total flux
-!! at the top of the atmosphere is not modified by time truncation error.
-!!
-!! ~~~~~~~~~~{.f90}
-!! real, intent(in), optional :: dt
-!! type(time_type), optional :: dt_time
-!! ~~~~~~~~~~
-!!
-!! @param [in] Latitudes of model grid points [radians]
-!! @param [in] Longitudes of model grid points [radians]
-!! @param [in] Time of day at longitude 0.0; midnight = 0.0, one day = 2 * pi [radians]
-!! @param [in] Time of year; autumnal equinox = 0.0, one year = 2 * pi [radians]
-!! @param [in] Time at which astronomical values are desired (time_type variable) [seconds, days]
-!! @param [out] Cosine of solar zenith angle, set to zero when entire period is in darkness [dimensionless]
-!! @param [out] Daylight fraction of time interval [dimensionless]
-!! @param [out] Earth-Sun distance (r) relative to semi-major axis of orbital ellipse
-!! (a):(a/r)**2 [dimensionless]
-!! @param [in] OPTIONAL: Time interval after gmt over which the astronomical variables are to be
-!! averaged. this produces averaged output rather than instantaneous. [radians], (1 day = 2 * pi)
-!! @param [in] OPTIONAL: Time interval after gmt over which the astronomical variables are to be
-!! averaged. this produces averaged output rather than instantaneous. time_type, [days, seconds]
-!! @param [in] Allow negative values for cosz?
-!! @param [out] half_day_out
-!> @ingroup astronomy_mod
-interface diurnal_solar
- module procedure diurnal_solar_2d
- module procedure diurnal_solar_1d
- module procedure diurnal_solar_0d
- module procedure diurnal_solar_cal_2d
- module procedure diurnal_solar_cal_1d
- module procedure diurnal_solar_cal_0d
-end interface
-
-!> @brief Calculates the daily mean solar information for a given time and latitude.
-!!
-!> ~~~~~~~~~~{.f90}
-!! call daily_mean_solar (lat, time, cosz, fracday, rrsun)
-!! call daily_mean_solar (lat, time_since_ae, cosz, fracday, rrsun)
-!! call daily_mean_solar (lat, time, cosz, solar)
-!! call daily_mean_solar (lat, time_since_ae, cosz, solar)
-!! ~~~~~~~~~~
-!!
-!! The first option (used in conjunction with time_manager_mod)
-!! generates the real variable time_since_ae from the time_type
-!! input time, and then calls daily_mean_solar with this real input
-!! (option 2). The third and fourth options correspond to the first
-!! and second and are used with then spectral 2-layer model, where
-!! only cosz and solar are desired as output. These routines generate
-!! dummy arguments and then call option 2, where the calculation is done.
-!!
-!! The time of year is set by
-!! ~~~~~~~~~~{.f90}
-!! real, intent(in) :: time_since_ae
-!! ~~~~~~~~~~
-!! With time_type input, the time of year is extracted from
-!! ~~~~~~~~~~{.f90}
-!! type(time_type), intent(in) :: time
-!! ~~~~~~~~~~
-!!
-!! Separate routines exist within this interface for scalar,
-!! 1D or 2D input and output fields:
-!!
-!! ~~~~~~~~~~{.f90}
-!! real, intent(in), dimension(:,:) :: lat
-!! real, intent(in), dimension(:) :: lat
-!! real, intent(in) :: lat
-!!
-!! real, intent(out), dimension(:,:) :: cosz, fracday
-!! real, intent(out), dimension(:) :: cosz, fracday
-!! real, intent(out) :: cosz, fracday
-!! ~~~~~~~~~~
-!!
-!! @param [in] Latitudes of model grid points [radians]
-!! @param [in] Time of year; autumnal equinox = 0.0, one year = 2 * pi [radians]
-!! @param [in] Time at which astronomical values are desired (time_type variable) [seconds, days]
-!! @param [out] Cosine of solar zenith angle, set to zero when entire period is in darkness [dimensionless]
-!! @param [out] Daylight fraction of time interval [dimensionless]
-!! @param [out] Earth-Sun distance (r) relative to semi-major axis of orbital ellipse
-!! (a):(a/r)**2 [dimensionless]
-!! @param [out] shortwave flux factor: cosine of zenith angle * daylight fraction /
-!! (earth-sun distance squared) [dimensionless]
-!> @ingroup astronomy_mod
-interface daily_mean_solar
- module procedure daily_mean_solar_2d
- module procedure daily_mean_solar_1d
- module procedure daily_mean_solar_2level
- module procedure daily_mean_solar_0d
- module procedure daily_mean_solar_cal_2d
- module procedure daily_mean_solar_cal_1d
- module procedure daily_mean_solar_cal_2level
- module procedure daily_mean_solar_cal_0d
-end interface
-
-!> Calculates the annual mean of solar information for a given latitude and time.
-!!
-!> ~~~~~~~~~~{.f90}
-!! call annual_mean_solar (js, je, lat, cosz, solar, fracday, rrsun)
-!! call annual_mean_solar (lat, cosz, solar)
-!! ~~~~~~~~~~
-!!
-!! The second interface above is used by the spectral 2-layer model,
-!! which requires only cosz and solar as output arguments, and which
-!! makes this call during the initialization phase of the model.
-!! Separate routines exist within this interface for 1D or 2D input
-!! and output fields:
-!!
-!! ~~~~~~~~~~{.f90}
-!! real, intent(in), dimension(:,:) :: lat
-!! real, intent(in), dimension(:) :: lat
-!!
-!! real, intent(out), dimension(:,:) :: cosz, solar, fracday
-!! real, intent(out), dimension(:) :: cosz, solar, fracday
-!! ~~~~~~~~~~
-!!
-!! @param [in] Starting subdomain j indices of data in the physics wiondow being integrated
-!! @param [in] Ending subdomain j indices of data in the physics wiondow being integrated
-!! @param [in] Latitudes of model grid points [radians]
-!! @param [out] cosz is the average over the year of the cosine of an effective zenith angle
-!! that would produce the correct daily solar flux if the sun were fixed at that
-!! single position for the period of daylight on the given day. in this average,
-!! the daily mean effective cosz is weighted by the daily mean solar flux. [dimensionless]
-!! @param [out] Normalized solar flux, averaged over the year, equal to the product of
-!! fracday*cosz*rrsun [dimensionless]
-!! @param [out] Daylight fraction calculated so as to make the average flux (solar) equal to the
-!! product of the flux-weighted avg cosz * this fracday * assumed annual mean avg
-!! Earth-Sun distance of 1.0. [dimensionless]
-!! @param [out] Annual mean Earth-Sun distance (r) relative to semi-major axis of orbital ellipse
-!! (a):(a/r)**2 [dimensionless]
-!> @ingroup astronomy_mod
-interface annual_mean_solar
- module procedure annual_mean_solar_2d
- module procedure annual_mean_solar_1d
- module procedure annual_mean_solar_2level
-end interface
-
-!> Gets the length of year for current calendar
-!!
-!> ~~~~~~~~~~{.f90}
-!! call get_period (period)
-!! ~~~~~~~~~~
-!!
-!! Separate routines exist within this interface for integer
-!! and time_type output:
-!!
-!! ~~~~~~~~~~{.f90}
-!! integer, intent(out) :: period
-!! type(time_type), intent(out) :: period
-!! ~~~~~~~~~~
-!!
-!! @param [out] Length of year for calendar in use
-!> @ingroup astronomy_mod
-interface get_period
- module procedure get_period_time_type, get_period_integer
-end interface
-
-!> Sets the length of a year for the calendar in use
-!!
-!> ~~~~~~~~~~{.f90}
-!! call set_period (period_in)
-!! ~~~~~~~~~~
-!!
-!! Separate routines exist within this interface for integer
-!! and time_type output:
-!!
-!! ~~~~~~~~~~{.f90}
-!! integer, intent(out) :: period_in
-!! type(time_type), intent(out) :: period_in
-!! ~~~~~~~~~~
-!!
-!! @param [in] Length of year for calendar in use
-!> @ingroup astronomy_mod
-interface set_period
- module procedure set_period_time_type, set_period_integer
-end interface
-
-
-private &
- orbit, & ! Called from astronomy_init and set_orbital_parameters
- r_inv_squared, & ! Called from diurnal_solar, daily_mean_solar and orbit
- angle, declination, half_day ! called from diurnal_solar and daily_mean_solar
-! half_day, orbital_time, & ! called from diurnal_solar and daily_mean_solar
-! universal_time ! called from diurnal_solar:
+ !---------------------------------------------------------------------
+ !------- interfaces --------
+
+ public &
+ astronomy_init, get_period, set_period, &
+ set_orbital_parameters, get_orbital_parameters, &
+ set_ref_date_of_ae, get_ref_date_of_ae, &
+ diurnal_solar, daily_mean_solar, annual_mean_solar, &
+ astronomy_end, universal_time, orbital_time
+
+ interface set_orbital_parameters
+ module procedure set_orbital_parameters_r4, set_orbital_parameters_r8
+ end interface set_orbital_parameters
+
+ interface get_orbital_parameters
+ module procedure get_orbital_parameters_r4, get_orbital_parameters_r8
+ end interface get_orbital_parameters
+
+ !> @}
+
+ !> @brief Calculates solar information for the given location(lat & lon) and time
+ !!
+ !> ~~~~~~~~~~{.f90}
+ !! call diurnal_solar (lat, lon, time, cosz, fracday, rrsun, dt_time)
+ !! call diurnal_solar (lat, lon, gmt, time_since_ae, cosz, fracday, rrsun, dt)
+ !! ~~~~~~~~~~
+ !!
+ !! The first option (used in conjunction with time_manager_mod)
+ !! generates the real variables gmt and time_since_ae from the
+ !! time_type input, and then calls diurnal_solar with these real inputs.
+ !!
+ !! The time of day is set by
+ !! ~~~~~~~~~~{.f90}
+ !! real, intent(in) :: gmt
+ !! ~~~~~~~~~~
+ !! The time of year is set by
+ !! ~~~~~~~~~~{.f90}
+ !! real, intent(in) :: time_since_ae
+ !! ~~~~~~~~~~
+ !! with time_type input, both of these are extracted from
+ !! ~~~~~~~~~~{.f90}
+ !! type(time_type), intent(in) :: time
+ !! ~~~~~~~~~~
+ !!
+ !! Separate routines exist within this interface for scalar,
+ !! 1D or 2D input and output fields:
+ !!
+ !! ~~~~~~~~~~{.f90}
+ !! real, intent(in), dimension(:,:) :: lat, lon
+ !! real, intent(in), dimension(:) :: lat, lon
+ !! real, intent(in) :: lat, lon
+ !!
+ !! real, intent(out), dimension(:,:) :: cosz, fracday
+ !! real, intent(out), dimension(:) :: cosz, fracday
+ !! real, intent(out) :: cosz, fracday
+ !! ~~~~~~~~~~
+ !!
+ !! One may also average the output fields over the time interval
+ !! between gmt and gmt + dt by including the optional argument dt (or
+ !! dt_time). dt is measured in radians and must be less than pi
+ !! (1/2 day). This average is computed analytically, and should be
+ !! exact except for the fact that changes in earth-sun distance over
+ !! the time interval dt are ignored. In the context of a diurnal GCM,
+ !! this option should always be employed to insure that the total flux
+ !! at the top of the atmosphere is not modified by time truncation error.
+ !!
+ !! ~~~~~~~~~~{.f90}
+ !! real, intent(in), optional :: dt
+ !! type(time_type), optional :: dt_time
+ !! ~~~~~~~~~~
+ !!
+ !! @param [in] Latitudes of model grid points [radians]
+ !! @param [in] Longitudes of model grid points [radians]
+ !! @param [in] Time of day at longitude 0.0; midnight = 0.0, one day = 2 * pi [radians]
+ !! @param [in] Time of year; autumnal equinox = 0.0, one year = 2 * pi [radians]
+ !! @param [in] Time at which astronomical values are desired (time_type variable) [seconds, days]
+ !! @param [out] Cosine of solar zenith angle, set to zero when entire period is in darkness [dimensionless]
+ !! @param [out] Daylight fraction of time interval [dimensionless]
+ !! @param [out] Earth-Sun distance (r) relative to semi-major axis of orbital ellipse
+ !! (a):(a/r)**2 [dimensionless]
+ !! @param [in] OPTIONAL: Time interval after gmt over which the astronomical variables are to be
+ !! averaged. this produces averaged output rather than instantaneous. [radians], (1 day = 2 * pi)
+ !! @param [in] OPTIONAL: Time interval after gmt over which the astronomical variables are to be
+ !! averaged. this produces averaged output rather than instantaneous. time_type,
+ !! [days, seconds]
+ !! @param [in] Allow negative values for cosz?
+ !! @param [out] half_day_out
+ !> @ingroup astronomy_mod
+ interface diurnal_solar
+ module procedure diurnal_solar_2d_r4, diurnal_solar_2d_r8
+ module procedure diurnal_solar_1d_r4, diurnal_solar_1d_r8
+ module procedure diurnal_solar_0d_r4, diurnal_solar_0d_r8
+ module procedure diurnal_solar_cal_2d_r4, diurnal_solar_cal_2d_r8
+ module procedure diurnal_solar_cal_1d_r4, diurnal_solar_cal_1d_r8
+ module procedure diurnal_solar_cal_0d_r4, diurnal_solar_cal_0d_r8
+ end interface diurnal_solar
+
+ !> @brief Calculates the daily mean solar information for a given time and latitude.
+ !!
+ !> ~~~~~~~~~~{.f90}
+ !! call daily_mean_solar (lat, time, cosz, fracday, rrsun)
+ !! call daily_mean_solar (lat, time_since_ae, cosz, fracday, rrsun)
+ !! call daily_mean_solar (lat, time, cosz, solar)
+ !! call daily_mean_solar (lat, time_since_ae, cosz, solar)
+ !! ~~~~~~~~~~
+ !!
+ !! The first option (used in conjunction with time_manager_mod)
+ !! generates the real variable time_since_ae from the time_type
+ !! input time, and then calls daily_mean_solar with this real input
+ !! (option 2). The third and fourth options correspond to the first
+ !! and second and are used with then spectral 2-layer model, where
+ !! only cosz and solar are desired as output. These routines generate
+ !! dummy arguments and then call option 2, where the calculation is done.
+ !!
+ !! The time of year is set by
+ !! ~~~~~~~~~~{.f90}
+ !! real, intent(in) :: time_since_ae
+ !! ~~~~~~~~~~
+ !! With time_type input, the time of year is extracted from
+ !! ~~~~~~~~~~{.f90}
+ !! type(time_type), intent(in) :: time
+ !! ~~~~~~~~~~
+ !!
+ !! Separate routines exist within this interface for scalar,
+ !! 1D or 2D input and output fields:
+ !!
+ !! ~~~~~~~~~~{.f90}
+ !! real, intent(in), dimension(:,:) :: lat
+ !! real, intent(in), dimension(:) :: lat
+ !! real, intent(in) :: lat
+ !!
+ !! real, intent(out), dimension(:,:) :: cosz, fracday
+ !! real, intent(out), dimension(:) :: cosz, fracday
+ !! real, intent(out) :: cosz, fracday
+ !! ~~~~~~~~~~
+ !!
+ !! @param [in] Latitudes of model grid points [radians]
+ !! @param [in] Time of year; autumnal equinox = 0.0, one year = 2 * pi [radians]
+ !! @param [in] Time at which astronomical values are desired (time_type variable) [seconds, days]
+ !! @param [out] Cosine of solar zenith angle, set to zero when entire period is in darkness [dimensionless]
+ !! @param [out] Daylight fraction of time interval [dimensionless]
+ !! @param [out] Earth-Sun distance (r) relative to semi-major axis of orbital ellipse
+ !! (a):(a/r)**2 [dimensionless]
+ !! @param [out] shortwave flux factor: cosine of zenith angle * daylight fraction /
+ !! (earth-sun distance squared) [dimensionless]
+ !> @ingroup astronomy_mod
+ interface daily_mean_solar
+ module procedure daily_mean_solar_2d_r4, daily_mean_solar_2d_r8
+ module procedure daily_mean_solar_1d_r4, daily_mean_solar_1d_r8
+ module procedure daily_mean_solar_2level_r4, daily_mean_solar_2level_r8
+ module procedure daily_mean_solar_0d_r4, daily_mean_solar_0d_r8
+ module procedure daily_mean_solar_cal_2d_r4, daily_mean_solar_cal_2d_r8
+ module procedure daily_mean_solar_cal_1d_r4, daily_mean_solar_cal_1d_r8
+ module procedure daily_mean_solar_cal_2level_r4, daily_mean_solar_cal_2level_r8
+ module procedure daily_mean_solar_cal_0d_r4, daily_mean_solar_cal_0d_r8
+ end interface daily_mean_solar
+
+ !> Calculates the annual mean of solar information for a given latitude and time.
+ !!
+ !> ~~~~~~~~~~{.f90}
+ !! call annual_mean_solar (js, je, lat, cosz, solar, fracday, rrsun)
+ !! call annual_mean_solar (lat, cosz, solar)
+ !! ~~~~~~~~~~
+ !!
+ !! The second interface above is used by the spectral 2-layer model,
+ !! which requires only cosz and solar as output arguments, and which
+ !! makes this call during the initialization phase of the model.
+ !! Separate routines exist within this interface for 1D or 2D input
+ !! and output fields:
+ !!
+ !! ~~~~~~~~~~{.f90}
+ !! real, intent(in), dimension(:,:) :: lat
+ !! real, intent(in), dimension(:) :: lat
+ !!
+ !! real, intent(out), dimension(:,:) :: cosz, solar, fracday
+ !! real, intent(out), dimension(:) :: cosz, solar, fracday
+ !! ~~~~~~~~~~
+ !!
+ !! @param [in] Starting subdomain j indices of data in the physics wiondow being integrated
+ !! @param [in] Ending subdomain j indices of data in the physics wiondow being integrated
+ !! @param [in] Latitudes of model grid points [radians]
+ !! @param [out] cosz is the average over the year of the cosine of an effective zenith angle
+ !! that would produce the correct daily solar flux if the sun were fixed at that
+ !! single position for the period of daylight on the given day. in this average,
+ !! the daily mean effective cosz is weighted by the daily mean solar flux. [dimensionless]
+ !! @param [out] Normalized solar flux, averaged over the year, equal to the product of
+ !! fracday*cosz*rrsun [dimensionless]
+ !! @param [out] Daylight fraction calculated so as to make the average flux (solar) equal to the
+ !! product of the flux-weighted avg cosz * this fracday * assumed annual mean avg
+ !! Earth-Sun distance of 1.0. [dimensionless]
+ !! @param [out] Annual mean Earth-Sun distance (r) relative to semi-major axis of orbital ellipse
+ !! (a):(a/r)**2 [dimensionless]
+ !> @ingroup astronomy_mod
+ interface annual_mean_solar
+ module procedure annual_mean_solar_2d_r4, annual_mean_solar_2d_r8
+ module procedure annual_mean_solar_1d_r4, annual_mean_solar_1d_r8
+ module procedure annual_mean_solar_2level_r4, annual_mean_solar_2level_r8
+ end interface annual_mean_solar
+
+ !> Gets the length of year for current calendar
+ !!
+ !> ~~~~~~~~~~{.f90}
+ !! call get_period (period)
+ !! ~~~~~~~~~~
+ !!
+ !! Separate routines exist within this interface for integer
+ !! and time_type output:
+ !!
+ !! ~~~~~~~~~~{.f90}
+ !! integer, intent(out) :: period
+ !! type(time_type), intent(out) :: period
+ !! ~~~~~~~~~~
+ !!
+ !! @param [out] Length of year for calendar in use
+ !> @ingroup astronomy_mod
+ interface get_period
+ module procedure get_period_time_type, get_period_integer
+ end interface
+
+ !> Sets the length of a year for the calendar in use
+ !!
+ !> ~~~~~~~~~~{.f90}
+ !! call set_period (period_in)
+ !! ~~~~~~~~~~
+ !!
+ !! Separate routines exist within this interface for integer
+ !! and time_type output:
+ !!
+ !! ~~~~~~~~~~{.f90}
+ !! integer, intent(out) :: period_in
+ !! type(time_type), intent(out) :: period_in
+ !! ~~~~~~~~~~
+ !!
+ !! @param [in] Length of year for calendar in use
+ !> @ingroup astronomy_mod
+ interface set_period
+ module procedure set_period_time_type, set_period_integer
+ end interface
+
+
+ private &
+ orbit, & ! Called from astronomy_init and set_orbital_parameters
+ r_inv_squared, & ! Called from diurnal_solar, daily_mean_solar and orbit
+ angle, declination, half_day ! called from diurnal_solar and daily_mean_solar
+ ! half_day, orbital_time, & ! called from diurnal_solar and daily_mean_solar
+ ! universal_time ! called from diurnal_solar:
+
+ interface r_inv_squared
+ module procedure r_inv_squared_r4, r_inv_squared_r8
+ end interface r_inv_squared
+
+ interface angle
+ module procedure angle_r4, angle_r8
+ end interface angle
+
+ interface declination
+ module procedure declination_r4, declination_r8
+ end interface declination
+
+ !> Private interface for internal use by dirunal_solar and daily_mean_solar.
+ !!
+ !> Example usage:
+ !! ~~~~~~~~~~{.f90}
+ !! half_day (latitude, dec) result (h)
+ !! ~~~~~~~~~~
+ !!
+ !! Separate routines exist within this interface for scalar,
+ !! or 2D input and output fields:
+ !!
+ !! ~~~~~~~~~~{.f90}
+ !! real, intent(in), dimension(:,:) :: latitude
+ !! real, intent(in) :: latitude
+ !!
+ !! real, dimension(size(latitude,1),size(latitude,2)) :: h
+ !! real :: h
+ !! ~~~~~~~~~~
+ !!
+ !! @param [in] Latitudes of model grid points [radians]
+ !! @param [in] Solar declination [radians]
+ !! @param [out] Half of the length of daylight at the given latitude and orbital position (dec); value
+ !! ranges between 0 (all darkness) and pi (all daylight) [dimensionless]
+ !> @ingroup astronomy_mod
+ interface half_day
+ module procedure half_day_2d_r4, half_day_2d_r8
+ module procedure half_day_0d_r4, half_day_0d_r8
+ end interface half_day
-!> Private interface for internal use by dirunal_solar and daily_mean_solar.
-!!
-!> Example usage:
-!! ~~~~~~~~~~{.f90}
-!! half_day (latitude, dec) result (h)
-!! ~~~~~~~~~~
-!!
-!! Separate routines exist within this interface for scalar,
-!! or 2D input and output fields:
-!!
-!! ~~~~~~~~~~{.f90}
-!! real, intent(in), dimension(:,:) :: latitude
-!! real, intent(in) :: latitude
-!!
-!! real, dimension(size(latitude,1),size(latitude,2)) :: h
-!! real :: h
-!! ~~~~~~~~~~
-!!
-!! @param [in] Latitudes of model grid points [radians]
-!! @param [in] Solar declination [radians]
-!! @param [out] Half of the length of daylight at the given latitude and orbital position (dec); value
-!! ranges between 0 (all darkness) and pi (all daylight) [dimensionless]
-!> @ingroup astronomy_mod
-interface half_day
- module procedure half_day_2d, half_day_0d
-end interface
!> @addtogroup astronomy_mod
!> @{
@@ -326,22 +350,22 @@ module astronomy_mod
!---------------------------------------------------------------------
!-------- namelist ---------
-real :: ecc = 0.01671 !< Eccentricity of Earth's orbit [dimensionless]
-real :: obliq = 23.439 !< Obliquity [degrees]
-real :: per = 102.932 !< Longitude of perihelion with respect
- !! to autumnal equinox in NH [degrees]
-integer :: period = 0 !< Specified length of year [seconds];
- !! must be specified to override default
- !! value given by length_of_year in
- !! time_manager_mod
-integer :: day_ae = 23 !< Day of specified autumnal equinox
-integer :: month_ae = 9 !< Month of specified autumnal equinox
-integer :: year_ae = 1998 !< Year of specified autumnal equinox
-integer :: hour_ae = 5 !< Hour of specified autumnal equinox
-integer :: minute_ae = 37 !< Minute of specified autumnal equinox
-integer :: second_ae = 0 !< Second of specified autumnal equinox
-integer :: num_angles = 3600 !< Number of intervals into which the year
- !! is divided to compute orbital positions
+real(r8_kind) :: ecc = 0.01671_r8_kind !< Eccentricity of Earth's orbit [dimensionless]
+real(r8_kind) :: obliq = 23.439_r8_kind !< Obliquity [degrees]
+real(r8_kind) :: per = 102.932_r8_kind !< Longitude of perihelion with respect
+ !! to autumnal equinox in NH [degrees]
+integer :: period = 0 !< Specified length of year [seconds];
+ !! must be specified to override default
+ !! value given by length_of_year in
+ !! time_manager_mod
+integer :: day_ae = 23 !< Day of specified autumnal equinox
+integer :: month_ae = 9 !< Month of specified autumnal equinox
+integer :: year_ae = 1998 !< Year of specified autumnal equinox
+integer :: hour_ae = 5 !< Hour of specified autumnal equinox
+integer :: minute_ae = 37 !< Minute of specified autumnal equinox
+integer :: second_ae = 0 !< Second of specified autumnal equinox
+integer :: num_angles = 3600 !< Number of intervals into which the year
+ !! is divided to compute orbital positions
namelist /astronomy_nml/ ecc, obliq, per, period, &
@@ -363,32 +387,32 @@ module astronomy_mod
type(time_type) :: period_time_type !< time_type variable containing
!! period of one orbit
-real, dimension(:), allocatable :: orb_angle !< table of orbital positions (0 to
+real(r8_kind), dimension(:), allocatable :: orb_angle !< table of orbital positions (0 to
!! 2*pi) as a function of time used
!! to find actual orbital position
!! via interpolation
-real :: seconds_per_day=86400. !< seconds in a day
-real :: deg_to_rad !< conversion from degrees to radians
-real :: twopi !< 2 *PI
-logical :: module_is_initialized=.false. !< has the module been initialized ?
+real(r8_kind) :: seconds_per_day = 86400.0_r8_kind !< seconds in a day
+real(r8_kind) :: deg_to_rad !< conversion from degrees to radians
+real(r8_kind) :: twopi !< 2 *PI
+logical :: module_is_initialized=.false. !< has the module been initialized ?
-real, dimension(:,:), allocatable :: &
+real(r8_kind), dimension(:,:), allocatable :: &
cosz_ann, & !< annual mean cos of zenith angle
solar_ann, & !< annual mean solar factor
fracday_ann !< annual mean daylight fraction
-real :: rrsun_ann !< annual mean earth-sun distance
-logical :: annual_mean_calculated=.false. !< have the annual mean values been calculated?
-integer :: num_pts = 0 !< count of grid_boxes for which
- !! annual mean astronomy values have
- !! been calculated
-integer :: total_pts !< number of grid boxes owned by the processor
+real(r8_kind) :: rrsun_ann !< annual mean earth-sun distance
+logical :: annual_mean_calculated=.false. !< have the annual mean values been calculated?
+integer :: num_pts = 0 !< count of grid_boxes for which
+ !! annual mean astronomy values have
+ !! been calculated
+integer :: total_pts !< number of grid boxes owned by the processor
!--------------------------------------------------------------------
- contains
+contains
@@ -405,8 +429,9 @@ module astronomy_mod
!! @throw FATAL, "astronomy_mod perihelion must be between 0 and 360 degrees"
subroutine astronomy_init (latb, lonb)
-real, dimension(:,:), intent(in), optional :: latb !< 2d array of model latitudes at cell corners [radians]
-real, dimension(:,:), intent(in), optional :: lonb !< 2d array of model longitudes at cell corners [radians]
+class(*), dimension(:,:), intent(in), optional :: latb !< 2d array of model latitudes at cell corners [radians]
+class(*), dimension(:,:), intent(in), optional :: lonb !< 2d array of model longitudes at cell corners [radians]
+logical :: is_valid
!-------------------------------------------------------------------
! local variables:
@@ -416,92 +441,115 @@ subroutine astronomy_init (latb, lonb)
!-------------------------------------------------------------------
! if module has already been initialized, exit.
!-------------------------------------------------------------------
- if (module_is_initialized) return
+ if (module_is_initialized) return
!-------------------------------------------------------------------
!> This routine will:
!> Verify that modules used by this module have been initialized.
!-------------------------------------------------------------------
- call fms_init
- call time_manager_init
- call constants_init
+ call fms_init
+ call time_manager_init
+ call constants_init
!-----------------------------------------------------------------------
!> Read namelist.
!-----------------------------------------------------------------------
- read (input_nml_file, astronomy_nml, iostat=io)
- ierr = check_nml_error(io,'astronomy_nml')
+ read (input_nml_file, astronomy_nml, iostat=io)
+ ierr = check_nml_error(io,'astronomy_nml')
!---------------------------------------------------------------------
!> Write version number and namelist to logfile.
!---------------------------------------------------------------------
- call write_version_number("ASTRONOMY_MOD", version)
- if (mpp_pe() == mpp_root_pe() ) then
- unit = stdlog()
- write (unit, nml=astronomy_nml)
- endif
+ call write_version_number("ASTRONOMY_MOD", version)
+ if (mpp_pe() == mpp_root_pe() ) then
+ unit = stdlog()
+ write (unit, nml=astronomy_nml)
+ endif
!--------------------------------------------------------------------
!> Be sure input values are within valid ranges.
! QUESTION : ARE THESE THE RIGHT LIMITS ???
!---------------------------------------------------------------------
- if (ecc < 0.0 .or. ecc > 0.99) &
- call error_mesg ('astronomy_mod', &
- 'ecc must be between 0 and 0.99', FATAL)
- if (obliq < -90. .or. obliq > 90.) &
+ if (ecc < 0.0_r8_kind .or. ecc > 0.99_r8_kind) &
+ call error_mesg ('astronomy_mod', &
+ 'ecc must be between 0 and 0.99', FATAL)
+ if (obliq < -90.0_r8_kind .or. obliq > 90.0_r8_kind) &
call error_mesg ('astronomy_mod', &
- 'obliquity must be between -90 and 90 degrees', FATAL)
- if (per < 0.0 .or. per > 360.0) &
+ 'obliquity must be between -90 and 90 degrees', FATAL)
+ if (per < 0.0_r8_kind .or. per > 360.0_r8_kind) &
call error_mesg ('astronomy_mod', &
- 'perihelion must be between 0 and 360 degrees', FATAL)
+ 'perihelion must be between 0 and 360 degrees', FATAL)
!----------------------------------------------------------------------
!> Set up time-type variable defining specified time of autumnal equinox.
!----------------------------------------------------------------------
- autumnal_eq_ref = set_date (year_ae,month_ae,day_ae, &
- hour_ae,minute_ae,second_ae)
+ autumnal_eq_ref = set_date (year_ae,month_ae,day_ae, &
+ hour_ae,minute_ae,second_ae)
!---------------------------------------------------------------------
!> Set up time-type variable defining length of year.
!----------------------------------------------------------------------
- if (period == 0) then
+ if (period == 0) then
period_time_type = length_of_year()
call get_time (period_time_type, seconds, days)
- period = int(seconds_per_day*days + seconds)
- else
+ period = int(seconds_per_day) * days + seconds
+ else
period_time_type = set_time(period,0)
- endif
+ endif
!---------------------------------------------------------------------
!> Define useful module variables.
!----------------------------------------------------------------------
- twopi = 2.*PI
- deg_to_rad = twopi/360.
+ twopi = 2.0_r8_kind * PI
+ deg_to_rad = twopi/360.0_r8_kind
!---------------------------------------------------------------------
!> Call orbit to define table of orbital angles as function of
!! orbital time.
!----------------------------------------------------------------------
! wfc moved here from orbit
- allocate ( orb_angle(0:num_angles) )
- call orbit
+ allocate ( orb_angle(0:num_angles) )
+ call orbit
!--------------------------------------------------------------------
!> If annual mean radiation is desired, then latb will be present.
!! allocate arrays to hold the needed astronomical factors. define
!! the total number of points that the processor is responsible for.
!--------------------------------------------------------------------
- if (present(latb)) then
- jd = size(latb,2) - 1
- id = size(lonb,1) - 1
- allocate (cosz_ann(id, jd))
- allocate (solar_ann(id, jd))
- allocate (fracday_ann(id, jd))
- total_pts = jd*id
- endif
+ ! check that no invalid types (integers or characters) are given as optional arg
+
+ is_valid = .false.
+ if (present(latb) .and. present(lonb)) then
+ select type (latb)
+ type is (real(r4_kind))
+ select type (lonb)
+ type is (real(r4_kind))
+ is_valid = .true.
+ end select
+ type is (real(r8_kind))
+ select type (lonb)
+ type is (real(r8_kind))
+ is_valid = .true.
+ end select
+ end select
+ if( is_valid ) then
+ jd = size(latb,2) - 1
+ id = size(lonb,1) - 1
+ allocate (cosz_ann(id, jd))
+ allocate (solar_ann(id, jd))
+ allocate (fracday_ann(id, jd))
+ total_pts = jd*id
+ else
+ call error_mesg('astronomy_mod', 'latb and/or lonb are not valid types', FATAL)
+ end if
+ elseif (present(latb) .and. .not. present(lonb)) then
+ call error_mesg ('astronomy_mod', 'lat and lon must both be present', FATAL)
+ elseif (present(lonb) .and. .not. present(latb)) then
+ call error_mesg ('astronomy_mod', 'lat and lon must both be present', FATAL)
+ endif
!---------------------------------------------------------------------
!> Mark the module as initialized.
!---------------------------------------------------------------------
- module_is_initialized=.true.
+ module_is_initialized=.true.
!---------------------------------------------------------------------
@@ -519,20 +567,19 @@ subroutine get_period_integer (period_out)
!--------------------------------------------------------------------
! local variables:
- integer :: seconds, days
+integer :: seconds, days
!---------------------------------------------------------------------
! exit if module has not been initialized.
!---------------------------------------------------------------------
- if (.not. module_is_initialized) &
- call error_mesg ( 'astronomy_mod', &
- ' module has not been initialized', FATAL)
+ if (.not. module_is_initialized) &
+ call error_mesg ( 'astronomy_mod', ' module has not been initialized', FATAL)
!--------------------------------------------------------------------
! define length of year in seconds.
!--------------------------------------------------------------------
- call get_time (period_time_type, seconds, days)
- period_out = int(seconds_per_day*days + seconds)
+ call get_time (period_time_type, seconds, days)
+ period_out = int(seconds_per_day) * days + seconds
end subroutine get_period_integer
@@ -548,14 +595,13 @@ subroutine get_period_time_type (period_out)
!---------------------------------------------------------------------
! exit if module has not been initialized.
!---------------------------------------------------------------------
- if (.not. module_is_initialized) &
- call error_mesg ( 'astronomy_mod', &
- ' module has not been initialized', FATAL)
+ if (.not. module_is_initialized) &
+ call error_mesg ('astronomy_mod', 'module has not been initialized', FATAL)
!--------------------------------------------------------------------
! define length of year as a time_type variable.
!--------------------------------------------------------------------
- period_out = period_time_type
+ period_out = period_time_type
end subroutine get_period_time_type
@@ -570,15 +616,14 @@ subroutine set_period_integer (period_in)
!---------------------------------------------------------------------
! exit if module has not been initialized.
!---------------------------------------------------------------------
- if (.not. module_is_initialized) &
- call error_mesg ( 'astronomy_mod', &
- ' module has not been initialized', FATAL)
+ if (.not. module_is_initialized) &
+ call error_mesg ('astronomy_mod', 'module has not been initialized', FATAL)
!---------------------------------------------------------------------
! define time_type variable defining the length of year from the
! input value (integer seconds).
!---------------------------------------------------------------------
- period_time_type = set_time(period_in, 0)
+ period_time_type = set_time(period_in, 0)
end subroutine set_period_integer
@@ -594,107 +639,18 @@ subroutine set_period_time_type(period_in)
!---------------------------------------------------------------------
! exit if module has not been initialized.
!---------------------------------------------------------------------
- if (.not. module_is_initialized) &
- call error_mesg ( 'astronomy_mod', &
- ' module has not been initialized', FATAL)
+ if (.not. module_is_initialized) &
+ call error_mesg ('astronomy_mod', 'module has not been initialized', FATAL)
!---------------------------------------------------------------------
! define time_type variable defining the length of year from the
! input value (time_type).
!---------------------------------------------------------------------
- period_time_type = period_in
+ period_time_type = period_in
end subroutine set_period_time_type
-!> @brief set_orbital_parameters saves the input values of eccentricity,
-!! obliquity and perihelion time as module variables for use by
-!! astronomy_mod.
-!!
-!! @throw FATAL, "astronomy_mod module has not been initialized"
-!! @throw FATAL, "astronomy_mod ecc must be between 0 and 0.99"
-!! @throw FATAL, "astronomy_mod obliquity must be between -90. and 90. degrees"
-!! @throw FATAL, "astronomy_mod perihelion must be between 0.0 and 360. degrees"
-subroutine set_orbital_parameters (ecc_in, obliq_in, per_in)
-
-real, intent(in) :: ecc_in !< Eccentricity of orbital ellipse [dimensionless]
-real, intent(in) :: obliq_in !< Obliquity [degrees]
-real, intent(in) :: per_in !< Longitude of perihelion with respect to autumnal
- !! equinox in northern hemisphere [degrees]
-
-!---------------------------------------------------------------------
-! exit if module has not been initialized.
-!---------------------------------------------------------------------
- if (.not. module_is_initialized) &
- call error_mesg ( 'astronomy_mod', &
- ' module has not been initialized', FATAL)
-
-!--------------------------------------------------------------------
-! be sure input values are within valid ranges.
-! QUESTION : ARE THESE THE RIGHT LIMITS ???
-!---------------------------------------------------------------------
- if (ecc_in < 0.0 .or. ecc_in > 0.99) &
- call error_mesg ('astronomy_mod', &
- 'ecc must be between 0 and 0.99', FATAL)
- if (obliq_in < -90.0 .or. obliq > 90.0) &
- call error_mesg ('astronomy_mod', &
- 'obliquity must be between -90. and 90. degrees', FATAL)
- if (per_in < 0.0 .or. per_in > 360.0) &
- call error_mesg ('astronomy_mod', &
- 'perihelion must be between 0.0 and 360. degrees', FATAL)
-
-!---------------------------------------------------------------------
-! save input values into module variables.
-!---------------------------------------------------------------------
- ecc = ecc_in
- obliq = obliq_in
- per = per_in
-
-!---------------------------------------------------------------------
-! call orbit to define table of orbital angles as function of
-! orbital time using the input values of parameters just supplied.
-!----------------------------------------------------------------------
- call orbit
-
-!----------------------------------------------------------------------
-
-end subroutine set_orbital_parameters
-
-!> @brief get_orbital_parameters retrieves the orbital parameters for use
-!! by another module.
-!!
-!! @throw FATAL, "astronomy_mod module has not been initialized"
-subroutine get_orbital_parameters (ecc_out, obliq_out, per_out)
-
-!-------------------------------------------------------------------
-! get_orbital_parameters retrieves the orbital parameters for use
-! by another module.
-!--------------------------------------------------------------------
-
-real, intent(out) :: ecc_out !< Eccentricity of orbital ellipse [dimensionless]
-real, intent(out) :: obliq_out !< Obliquity [degrees]
-real, intent(out) :: per_out !< Longitude of perihelion with respect to autumnal
- !! equinox in northern hemisphere [degrees]
-
-!---------------------------------------------------------------------
-! exit if module has not been initialized.
-!---------------------------------------------------------------------
- if (.not. module_is_initialized) &
- call error_mesg ( 'astronomy_mod', &
- ' module has not been initialized', FATAL)
-
-!--------------------------------------------------------------------
-! fill the output arguments with the eccentricity, obliquity and
-! perihelion angle.
-!--------------------------------------------------------------------
- ecc_out = ecc
- obliq_out = obliq
- per_out = per
-
-
-end subroutine get_orbital_parameters
-
-
!> @brief set_ref_date_of_ae provides a means of specifying the reference
!! date of the NH autumnal equinox for a particular year.
!!
@@ -714,6 +670,7 @@ end subroutine get_orbital_parameters
!! @param [out] OPTIONAL: Hour of reference autumnal equinox
!!
!! @throw FATAL, "astronomy_mod module has not been initialized"
+
subroutine set_ref_date_of_ae (day_in,month_in,year_in, &
second_in,minute_in,hour_in)
@@ -723,30 +680,29 @@ subroutine set_ref_date_of_ae (day_in,month_in,year_in, &
!---------------------------------------------------------------------
! exit if module has not been initialized.
!---------------------------------------------------------------------
- if (.not. module_is_initialized) &
- call error_mesg ( 'astronomy_mod', &
- ' module has not been initialized', FATAL)
+ if (.not. module_is_initialized) &
+ call error_mesg ('astronomy_mod', 'module has not been initialized', FATAL)
!--------------------------------------------------------------------
! save the input time of ae specification into a time_type module
! variable autumnal_eq_ref.
!--------------------------------------------------------------------
- day_ae = day_in
- month_ae = month_in
- year_ae = year_in
+ day_ae = day_in
+ month_ae = month_in
+ year_ae = year_in
- if (present(second_in)) then
+ if (present(second_in)) then
second_ae = second_in
minute_ae = minute_in
hour_ae = hour_in
- else
+ else
second_ae = 0
minute_ae = 0
hour_ae = 0
- endif
+ endif
- autumnal_eq_ref = set_date (year_ae,month_ae,day_ae, &
- hour_ae,minute_ae,second_ae)
+ autumnal_eq_ref = set_date (year_ae,month_ae,day_ae, &
+ hour_ae,minute_ae,second_ae)
!---------------------------------------------------------------------
@@ -774,1476 +730,132 @@ subroutine get_ref_date_of_ae (day_out,month_out,year_out,&
!---------------------------------------------------------------------
! exit if module has not been initialized.
!---------------------------------------------------------------------
- if (.not. module_is_initialized) &
- call error_mesg ( 'astronomy_mod', &
- ' module has not been initialized', FATAL)
+ if (.not. module_is_initialized) &
+ call error_mesg ('astronomy_mod', 'module has not been initialized', FATAL)
!---------------------------------------------------------------------
! fill the output fields with the proper module data.
!---------------------------------------------------------------------
- day_out = day_ae
- month_out = month_ae
- year_out = year_ae
- second_out = second_ae
- minute_out = minute_ae
- hour_out = hour_ae
+ day_out = day_ae
+ month_out = month_ae
+ year_out = year_ae
+ second_out = second_ae
+ minute_out = minute_ae
+ hour_out = hour_ae
end subroutine get_ref_date_of_ae
+!> @brief astronomy_end is the destructor for astronomy_mod.
+subroutine astronomy_end
-!> @brief diurnal_solar_2d returns 2d fields of cosine of zenith angle,
-!! daylight fraction and earth-sun distance at the specified latitudes,
-!! longitudes and time. These values may be instantaneous or averaged
-!! over a specified time interval.
-!!
-!! @param [in] Latitudes of model grid points
-!! @param [in] Longitudes of model grid points
-!! @param [in] Time of day at longitude 0.0; midnight = 0.0, one day = 2 * pi
-!! @param [in] Time of year; autumnal equinox = 0.0, one year = 2 * pi
-!! @param [out] Cosine of solar zenith angle
-!! @param [out] Daylight fraction of time interval
-!! @param [out] earth-sun distance (r) relative to semi-major axis of orbital ellipse (a):(a/r)**2
-!! @param [in] OPTIONAL: time interval after gmt over which the astronomical variables are to be
-!! averaged. this produces averaged output rather than instantaneous.
-!! @param [in] Allow negative values for cosz?
-!! @param [out] half_day_out
-!!
-!! @throw FATAL, "astronomy_mod time_since_ae not between 0 and 2pi"
-!! @throw FATAL, "astronomy_mod gmt not between 0 and 2pi"
-subroutine diurnal_solar_2d (lat, lon, gmt, time_since_ae, cosz, &
- fracday, rrsun, dt, allow_negative_cosz, &
- half_day_out)
-
-real, dimension(:,:), intent(in) :: lat, lon
-real, intent(in) :: gmt, time_since_ae
-real, dimension(:,:), intent(out) :: cosz, fracday
-real, intent(out) :: rrsun
-real, intent(in), optional :: dt
-logical, intent(in), optional :: allow_negative_cosz
-real, dimension(:,:), intent(out), optional :: half_day_out
-
+ !----------------------------------------------------------------------
+ !> check if the module has been initialized.
+ !----------------------------------------------------------------------
+ if (.not. module_is_initialized) return
+ ! call error_mesg ( 'astronomy_mod', &
+ ! ' module has not been initialized', FATAL)
+
+ !----------------------------------------------------------------------
+ !> deallocate module variables.
+ !----------------------------------------------------------------------
+ deallocate (orb_angle)
+ if (allocated(cosz_ann) ) then
+ deallocate (cosz_ann)
+ deallocate (fracday_ann)
+ deallocate (solar_ann)
+ endif
-!---------------------------------------------------------------------
-! local variables
+ !----------------------------------------------------------------------
+ !> mark the module as uninitialized.
+ !----------------------------------------------------------------------
+ module_is_initialized = .false.
- real, dimension(size(lat,1),size(lat,2)) :: t, tt, h, aa, bb, &
- st, stt, sh
- real :: ang, dec
- logical :: Lallow_negative
+end subroutine astronomy_end
-!---------------------------------------------------------------------
-! local variables
+!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!
-! t time of day with respect to local noon (2 pi = 1 day)
-! [ radians ]
-! tt end of averaging period [ radians ]
-! h half of the daily period of daylight, centered at noon
-! [ radians, -pi --> pi ]
-! aa sin(lat) * sin(declination)
-! bb cos(lat) * cos(declination)
-! st sine of time of day
-! stt sine of time of day at end of averaging period
-! sh sine of half-day period
-! ang position of earth in its orbit wrt autumnal equinox
-! [ radians ]
-! dec earth's declination [ radians ]
+! PRIVATE SUBROUTINES
!
-!--------------------------------------------------------------------
-
-!--------------------------------------------------------------------
-! be sure the time in the annual cycle is legitimate.
-!---------------------------------------------------------------------
- if (time_since_ae < 0.0 .or. time_since_ae > twopi) &
- call error_mesg('astronomy_mod', &
- 'time_since_ae not between 0 and 2pi', FATAL)
-
-!--------------------------------------------------------------------
-! be sure the time at longitude = 0.0 is legitimate.
-!---------------------------------------------------------------------
- if (gmt < 0.0 .or. gmt > twopi) &
- call error_mesg('astronomy_mod', &
- 'gmt not between 0 and 2pi', FATAL)
-
-!---------------------------------------------------------------------
-!> define the orbital angle (location in year), solar declination and
-!! earth sun distance factor. use functions contained in this module.
-!---------------------------------------------------------------------
- ang = angle(time_since_ae)
- dec = declination(ang)
- rrsun = r_inv_squared(ang)
-
-!---------------------------------------------------------------------
-!> define terms needed in the cosine zenith angle equation.
-!--------------------------------------------------------------------
- aa = sin(lat)*sin(dec)
- bb = cos(lat)*cos(dec)
-
-!---------------------------------------------------------------------
-!> define local time. force it to be between -pi and pi.
-!--------------------------------------------------------------------
- t = gmt + lon - PI
- where(t >= PI) t = t - twopi
- where(t < -PI) t = t + twopi
-
- Lallow_negative = .false.
- if (present(allow_negative_cosz)) then
- if (allow_negative_cosz) Lallow_negative = .true.
- end if
-
-!---------------------------------------------------------------------
-!> perform a time integration to obtain cosz and fracday if desired.
-!! output is valid over the period from t to t + dt.
-!--------------------------------------------------------------------
- h = half_day (lat,dec)
-
- if ( present(half_day_out) ) then
- half_day_out = h
- end if
-
- if ( present(dt) ) then ! (perform time averaging)
- tt = t + dt
- st = sin(t)
- stt = sin(tt)
- sh = sin(h)
- cosz = 0.0
-
- if (.not. Lallow_negative) then
-!-------------------------------------------------------------------
-! case 1: entire averaging period is before sunrise.
-!-------------------------------------------------------------------
- where (t < -h .and. tt < -h) cosz = 0.0
-
-!-------------------------------------------------------------------
-! case 2: averaging period begins before sunrise, ends after sunrise
-! but before sunset
-!-------------------------------------------------------------------
- where ( (tt+h) /= 0.0 .and. t < -h .and. abs(tt) <= h) &
- cosz = aa + bb*(stt + sh)/ (tt + h)
-
-!-------------------------------------------------------------------
-! case 3: averaging period begins before sunrise, ends after sunset,
-! but before the next sunrise. modify if averaging period extends
-! past the next day's sunrise, but if averaging period is less than
-! a half- day (pi) that circumstance will never occur.
-!-------------------------------------------------------------------
- where (t < -h .and. h /= 0.0 .and. h < tt) &
- cosz = aa + bb*( sh + sh)/(h+h)
-
-!-------------------------------------------------------------------
-! case 4: averaging period begins after sunrise, ends before sunset.
-!-------------------------------------------------------------------
- where ( abs(t) <= h .and. abs(tt) <= h) &
- cosz = aa + bb*(stt - st)/ (tt - t)
-
-!-------------------------------------------------------------------
-! case 5: averaging period begins after sunrise, ends after sunset.
-! modify when averaging period extends past the next day's sunrise.
-!-------------------------------------------------------------------
- where ((h-t) /= 0.0 .and. abs(t) <= h .and. h < tt) &
- cosz = aa + bb*(sh - st)/(h-t)
-
-!-------------------------------------------------------------------
-! case 6: averaging period begins after sunrise , ends after the
-! next day's sunrise. note that this includes the case when the
-! day length is one day (h = pi).
-!-------------------------------------------------------------------
- where (twopi - h < tt .and. (tt+h-twopi) /= 0.0 .and. t <= h ) &
- cosz = (cosz*(h - t) + (aa*(tt + h - twopi) + &
- bb*(stt + sh))) / ((h - t) + (tt + h - twopi))
-
-!-------------------------------------------------------------------
-! case 7: averaging period begins after sunset and ends before the
-! next day's sunrise
-!-------------------------------------------------------------------
- where( h < t .and. twopi - h >= tt ) cosz = 0.0
-
-!-------------------------------------------------------------------
-! case 8: averaging period begins after sunset and ends after the
-! next day's sunrise but before the next day's sunset. if the
-! averaging period is less than a half-day (pi) the latter
-! circumstance will never occur.
-!-----------------------------------------------------------------
- where( h < t .and. twopi - h < tt )
- cosz = aa + bb*(stt + sh) / (tt + h - twopi)
- end where
-
- else
- cosz = aa + bb*(stt - st)/ (tt - t)
- end if
-
-
-
-!-------------------------------------------------------------------
-! day fraction is the fraction of the averaging period contained
-! within the (-h,h) period.
-!-------------------------------------------------------------------
- where (t < -h .and. tt < -h) fracday = 0.0
- where (t < -h .and. abs(tt) <= h) fracday = (tt + h )/dt
- where (t < -h .and. h < tt) fracday = ( h + h )/dt
- where (abs(t) <= h .and. abs(tt) <= h) fracday = (tt - t )/dt
- where (abs(t) <= h .and. h < tt) fracday = ( h - t )/dt
- where ( h < t ) fracday = 0.0
- where (twopi - h < tt) fracday = fracday + &
- (tt + h - &
- twopi)/dt
-!----------------------------------------------------------------------
-!> if instantaneous values are desired, define cosz at time t.
-!----------------------------------------------------------------------
- else ! (no time averaging)
- if (.not. Lallow_negative) then
- where (abs(t) < h)
- cosz = aa + bb*cos(t)
- fracday = 1.0
- elsewhere
- cosz = 0.0
- fracday = 0.0
- end where
- else
- cosz = aa + bb*cos(t)
- where (abs(t) < h)
- fracday = 1.0
- elsewhere
- fracday = 0.0
- end where
- end if
- end if
-
-!----------------------------------------------------------------------
-!> Check that cosz is not negative, if desired.
-!----------------------------------------------------------------------
- if (.not. Lallow_negative) then
- cosz = max(0.0, cosz)
- end if
-
-end subroutine diurnal_solar_2d
-
-
-!> @brief diurnal_solar_1d takes 1-d input fields, adds a second dimension
-!! and calls diurnal_solar_2d. on return, the 2d fields are returned
-!! to the original 1d fields.
-!!
-!! @param [in] Latitudes of model grid points
-!! @param [in] Longitudes of model grid points
-!! @param [in] Time of day at longitude 0.0; midnight = 0.0, one day = 2 * pi
-!! @param [in] Time of year; autumnal equinox = 0.0, one year = 2 * pi
-!! @param [out] Cosine of solar zenith angle
-!! @param [out] Daylight fraction of time interval
-!! @param [out] earth-sun distance (r) relative to semi-major axis of orbital ellipse (a):(a/r)**2
-!! @param [in] OPTIONAL: time interval after gmt over which the astronomical variables are to be
-!! averaged. this produces averaged output rather than instantaneous.
-!! @param [in] Allow negative values for cosz?
-!! @param [out] half_day_out
-subroutine diurnal_solar_1d (lat, lon, gmt, time_since_ae, cosz, &
- fracday, rrsun, dt, allow_negative_cosz, &
- half_day_out)
-
-!---------------------------------------------------------------------
-real, dimension(:), intent(in) :: lat, lon
-real, intent(in) :: gmt, time_since_ae
-real, dimension(:), intent(out) :: cosz, fracday
-real, intent(out) :: rrsun
-real, intent(in), optional :: dt
-logical, intent(in), optional :: allow_negative_cosz
-real, dimension(:), intent(out), optional :: half_day_out
+!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+!> @brief Orbit computes and stores a table of value of orbital angles as a
+!! function of orbital time (both the angle and time are zero at
+!! autumnal equinox in the NH, and range from 0 to 2*pi).
+subroutine orbit
!---------------------------------------------------------------------
-! local variables
-!---------------------------------------------------------------------
- real, dimension(size(lat),1) :: lat_2d, lon_2d, cosz_2d, &
- fracday_2d,halfday_2d
-
-!--------------------------------------------------------------------
-!> define 2-d versions of input data arrays.
-!--------------------------------------------------------------------
- lat_2d(:,1) = lat
- lon_2d(:,1) = lon
-
-!--------------------------------------------------------------------
-!> call diurnal_solar_2d to calculate astronomy fields.
-!--------------------------------------------------------------------
-! if (present(dt)) then
- call diurnal_solar_2d (lat_2d, lon_2d, gmt, time_since_ae,&
- cosz_2d, fracday_2d, rrsun, dt=dt, &
- allow_negative_cosz=allow_negative_cosz, &
- half_day_out=halfday_2d)
-! else
-! call diurnal_solar_2d (lat_2d, lon_2d, gmt, time_since_ae, &
-! cosz_2d, fracday_2d, rrsun)
-! endif
-
-!-------------------------------------------------------------------
-!> place output fields into 1-d arguments for return to
-!! calling routine.
-!-------------------------------------------------------------------
- fracday = fracday_2d(:,1)
- cosz = cosz_2d (:,1)
- if (present(half_day_out)) then
- half_day_out = halfday_2d(:,1)
- end if
-
-end subroutine diurnal_solar_1d
-
-
-!> @brief diurnal_solar_0d takes scalar input fields, makes them into 2d
-!! arrays dimensioned (1,1), and calls diurnal_solar_2d. on return,
-!! the 2d fields are converted back to the desired scalar output.
-!!
-!! @param [in] Latitudes of model grid points
-!! @param [in] Longitudes of model grid points
-!! @param [in] Time of day at longitude 0.0; midnight = 0.0, one day = 2 * pi
-!! @param [in] Time of year; autumnal equinox = 0.0, one year = 2 * pi
-!! @param [out] Cosine of solar zenith angle
-!! @param [out] Daylight fraction of time interval
-!! @param [out] earth-sun distance (r) relative to semi-major axis of orbital ellipse (a):(a/r)**2
-!! @param [in] OPTIONAL: time interval after gmt over which the astronomical variables are to be
-!! averaged. this produces averaged output rather than instantaneous.
-!! @param [in] Allow negative values for cosz?
-!! @param [out] half_day_out
-subroutine diurnal_solar_0d (lat, lon, gmt, time_since_ae, cosz, &
- fracday, rrsun, dt, allow_negative_cosz, &
- half_day_out)
-
-real, intent(in) :: lat, lon, gmt, time_since_ae
-real, intent(out) :: cosz, fracday, rrsun
-real, intent(in), optional :: dt
-logical,intent(in),optional :: allow_negative_cosz
-real, intent(out), optional :: half_day_out
-
-!--------------------------------------------------------------------
-! local variables:
-!--------------------------------------------------------------------
- real, dimension(1,1) :: lat_2d, lon_2d, cosz_2d, fracday_2d, halfday_2d
+! local variables
-!---------------------------------------------------------------------
-!> create 2d arrays from the scalar input fields.
-!---------------------------------------------------------------------
- lat_2d = lat
- lon_2d = lon
+integer :: n
+real(kind=r8_kind) :: d1, d2, d3, d4, d5, dt, norm
!--------------------------------------------------------------------
-!> call diurnal_solar_2d to calculate astronomy fields.
+!> allocate the orbital angle array, sized by the namelist parameter
+!! num_angles, defining the annual cycle resolution of the earth's
+!! orbit. define some constants to be used.
!--------------------------------------------------------------------
-! if (present(dt)) then
- call diurnal_solar_2d (lat_2d, lon_2d, gmt, time_since_ae, &
- cosz_2d, fracday_2d, rrsun, dt=dt, &
- allow_negative_cosz=allow_negative_cosz, &
- half_day_out=halfday_2d)
-! else
-! call diurnal_solar_2d (lat_2d, lon_2d, gmt, time_since_ae, &
-! cosz_2d, fracday_2d, rrsun)
-! end if
-
-!-------------------------------------------------------------------
-!> place output fields into scalars for return to calling routine.
-!-------------------------------------------------------------------
- fracday = fracday_2d(1,1)
- cosz = cosz_2d(1,1)
- if (present(half_day_out)) then
- half_day_out = halfday_2d(1,1)
- end if
-
-end subroutine diurnal_solar_0d
-
-
-!> @brief diurnal_solar_cal_2d receives time_type inputs, converts
-!! them to real variables and then calls diurnal_solar_2d to
-!! compute desired astronomical variables.
-!!
-!! @param [in] Latitudes of model grid points
-!! @param [in] Longitudes of model grid points
-!! @param [in] Time of day at longitude 0.0; midnight = 0.0, one day = 2 * pi
-!! @param [in] Time of year (time_type)
-!! @param [out] Cosine of solar zenith angle
-!! @param [out] Daylight fraction of time interval
-!! @param [out] earth-sun distance (r) relative to semi-major axis of orbital ellipse (a):(a/r)**2
-!! @param [in] OPTIONAL: time interval after gmt over which the astronomical variables are to be
-!! averaged. this produces averaged output rather than instantaneous.
-!! @param [in] Allow negative values for cosz?
-!! @param [out] half_day_out
-!!
-!! @throw FATAL, "astronomy_mod radiation time step must be no longer than 12 hrs"
-!! @throw FATAL, "astronomy_mod radiation time step must not be an integral number of days"
-subroutine diurnal_solar_cal_2d (lat, lon, time, cosz, fracday, &
- rrsun, dt_time, allow_negative_cosz, &
- half_day_out)
-
-!-------------------------------------------------------------------
-real, dimension(:,:), intent(in) :: lat, lon
-type(time_type), intent(in) :: time
-real, dimension(:,:), intent(out) :: cosz, fracday
-real, intent(out) :: rrsun
-type(time_type), intent(in), optional :: dt_time
-logical, intent(in), optional :: allow_negative_cosz
-real, dimension(:,:), intent(out), optional :: half_day_out
-!---------------------------------------------------------------------
-
-!---------------------------------------------------------------------
-! local variables
-!---------------------------------------------------------------------
- real :: dt
- real :: gmt, time_since_ae
+! wfc moving to astronomy_init
+! allocate ( orb_angle(0:num_angles) )
+orb_angle(0) = 0.0_r8_kind
+dt = twopi/real(num_angles, r8_kind)
+norm = sqrt(1.0_r8_kind - ecc**2)
+dt = dt*norm
!---------------------------------------------------------------------
-!> Extract time of day (gmt) from time_type variable time with
-!! function universal_time.
+!> define the orbital angle at each of the num_angles locations in
+!! the orbit.
!---------------------------------------------------------------------
- gmt = universal_time(time)
+ do n = 1, num_angles
+ d1 = dt*r_inv_squared(orb_angle(n-1))
+ d2 = dt*r_inv_squared(orb_angle(n-1) + 0.5_r8_kind * d1)
+ d3 = dt*r_inv_squared(orb_angle(n-1) + 0.5_r8_kind * d2)
+ d4 = dt*r_inv_squared(orb_angle(n-1) + d3)
+ d5 = d1/6.0_r8_kind + d2/3.0_r8_kind + d3/3.0_r8_kind + d4/6.0_r8_kind
+ orb_angle(n) = orb_angle(n-1) + d5
+ end do
-!---------------------------------------------------------------------
-!> Extract the time of year (time_since_ae) from time_type variable
-!! time using the function orbital_time.
-!---------------------------------------------------------------------
- time_since_ae = orbital_time(time)
+end subroutine orbit
-!---------------------------------------------------------------------
-!> Convert optional time_type variable dt_time (length of averaging
-!! period) to a real variable dt with the function universal_time.
-!---------------------------------------------------------------------
- if (present(dt_time)) then
- dt = universal_time(dt_time)
- if (dt > PI) then
- call error_mesg ( 'astronomy_mod', &
- 'radiation time step must be no longer than 12 hrs', &
- FATAL)
- endif
- if (dt == 0.0) then
- call error_mesg ( 'astronomy_mod', &
- 'radiation time step must not be an integral &
- &number of days', FATAL)
- endif
-!--------------------------------------------------------------------
-!> Call diurnal_solar_2d to calculate astronomy fields, with or
-!! without the optional argument dt.
-!--------------------------------------------------------------------
- call diurnal_solar_2d (lat, lon, gmt, time_since_ae, cosz, &
- fracday, rrsun, dt=dt, &
- allow_negative_cosz=allow_negative_cosz, &
- half_day_out=half_day_out)
- else
- call diurnal_solar_2d (lat, lon, gmt, time_since_ae, cosz, &
- fracday, rrsun, &
- allow_negative_cosz=allow_negative_cosz, &
- half_day_out=half_day_out)
- end if
-
-end subroutine diurnal_solar_cal_2d
-
-
-!> @brief diurnal_solar_cal_1d receives time_type inputs, converts
-!! them to real variables and then calls diurnal_solar_2d to
-!! compute desired astronomical variables.
+!> @brief Orbital time returns the time (1 year = 2*pi) since autumnal
+!! equinox
!!
-!! @param [in] Latitudes of model grid points
-!! @param [in] Longitudes of model grid points
-!! @param [in] Time of day at longitude 0.0; midnight = 0.0, one day = 2 * pi
-!! @param [in] Time of year (time_type)
-!! @param [out] Cosine of solar zenith angle
-!! @param [out] Daylight fraction of time interval
-!! @param [out] earth-sun distance (r) relative to semi-major axis of orbital ellipse (a):(a/r)**2
-!! @param [in] OPTIONAL: time interval after gmt over which the astronomical variables are to be
-!! averaged. this produces averaged output rather than instantaneous.
-!! @param [in] Allow negative values for cosz?
-!! @param [out] half_day_out
-subroutine diurnal_solar_cal_1d (lat, lon, time, cosz, fracday, &
- rrsun, dt_time, allow_negative_cosz, &
- half_day_out)
+!! @details Orbital time returns the time (1 year = 2*pi) since autumnal
+!! equinox; autumnal_eq_ref is a module variable of time_type and
+!! will have been defined by default or by a call to
+!! set_ref_date_of_ae; length_of_year is available through the time
+!! manager and is set at the value approriate for the calandar being used
+function orbital_time(time) result(t)
-!--------------------------------------------------------------------
-real, dimension(:), intent(in) :: lat, lon
-type(time_type), intent(in) :: time
-real, dimension(:), intent(out) :: cosz, fracday
-real, intent(out) :: rrsun
-type(time_type), intent(in), optional :: dt_time
-logical, intent(in), optional :: allow_negative_cosz
-real, dimension(:), intent(out), optional :: half_day_out
-!--------------------------------------------------------------------
+type(time_type), intent(in) :: time !< time (1 year = 2*pi) since autumnal equinox
+real(kind=r8_kind) :: t
-!-------------------------------------------------------------------
-! local variables
-!-------------------------------------------------------------------
- real, dimension(size(lat),1) :: lat_2d, lon_2d, cosz_2d, &
- fracday_2d, halfday_2d
+ t = (time - autumnal_eq_ref)//period_time_type
+ t = twopi*(t - real(floor(t), r8_kind))
+ if (time < autumnal_eq_ref) t = twopi - t
-!--------------------------------------------------------------------
-!> Define 2-d versions of input data arrays.
-!--------------------------------------------------------------------
- lat_2d(:,1) = lat
- lon_2d(:,1) = lon
+end function orbital_time
-!--------------------------------------------------------------------
-!> Call diurnal_solar_cal_2d to convert the time_types to reals and
-!! then calculate the astronomy fields.
-!--------------------------------------------------------------------
- if (present(dt_time)) then
- call diurnal_solar_cal_2d (lat_2d, lon_2d, time, cosz_2d, &
- fracday_2d, rrsun, dt_time=dt_time, &
- allow_negative_cosz=allow_negative_cosz, &
- half_day_out=halfday_2d)
- else
- call diurnal_solar_cal_2d (lat_2d, lon_2d, time, cosz_2d, &
- fracday_2d, rrsun, &
- allow_negative_cosz=allow_negative_cosz, &
- half_day_out=halfday_2d)
- end if
-!-------------------------------------------------------------------
-!> Place output fields into 1-d arguments for return to
-!! calling routine.
-!-------------------------------------------------------------------
- fracday = fracday_2d(:,1)
- cosz = cosz_2d (:,1)
- if (present(half_day_out)) then
- half_day_out = halfday_2d(:,1)
- end if
+!> @brief universal_time returns the time of day at longitude = 0.0
+!! (1 day = 2*pi)
+function universal_time(time) result(t)
-end subroutine diurnal_solar_cal_1d
+type(time_type), intent(in) :: time !< Time (1 year = 2*pi) since autumnal equinox
+real(kind=r8_kind) :: t
+ !--------------------------------------------------------------------
+ ! local variables
+ !--------------------------------------------------------------------
+ integer :: seconds, days
-!> @brief diurnal_solar_cal_0d receives time_type inputs, converts them to real variables
-!! and then calls diurnal_solar_2d to compute desired astronomical variables.
-!!
-!! @param [in] Latitudes of model grid points
-!! @param [in] Longitudes of model grid points
-!! @param [in] Time of year (time_type)
-!! @param [out] Cosine of solar zenith angle
-!! @param [out] Daylight fraction of time interval
-!! @param [out] Earth-Sun distance (r) relative to semi-major axis of orbital ellipse (a) : (a/r)**2
-!! @param [out] OPTIONAL: time interval after gmt over which the astronomical variables are
-!! to be averaged. this produces averaged output rather than instantaneous.
-!! @param [in] allow_negative_cosz
-!! @param [out] half_day_out
-subroutine diurnal_solar_cal_0d (lat, lon, time, cosz, fracday, &
- rrsun, dt_time, allow_negative_cosz, &
- half_day_out)
+ call get_time(time, seconds, days)
+ t = twopi * real(seconds, r8_kind)/real(seconds_per_day, r8_kind)
-!---------------------------------------------------------------------
-real, intent(in) :: lat, lon
-type(time_type), intent(in) :: time
-real, intent(out) :: cosz, fracday, rrsun
-type(time_type), intent(in), optional :: dt_time
-logical, intent(in), optional :: allow_negative_cosz
-real, intent(out), optional :: half_day_out
-!---------------------------------------------------------------------
+end function universal_time
-!---------------------------------------------------------------------
-! local variables
-!---------------------------------------------------------------------
- real, dimension(1,1) :: lat_2d, lon_2d, cosz_2d, fracday_2d, halfday_2d
-
-!--------------------------------------------------------------------
-!> Define 2-d versions of input data arrays.
-!--------------------------------------------------------------------
- lat_2d = lat
- lon_2d = lon
-
-!--------------------------------------------------------------------
-!> Call diurnal_solar_cal_2d to convert the time_types to reals and
-!! then calculate the astronomy fields.
-!--------------------------------------------------------------------
- if (present(dt_time)) then
- call diurnal_solar_cal_2d (lat_2d, lon_2d, time, cosz_2d, &
- fracday_2d, rrsun, dt_time=dt_time, &
- allow_negative_cosz=allow_negative_cosz, &
- half_day_out=halfday_2d)
- else
- call diurnal_solar_cal_2d (lat_2d, lon_2d, time, cosz_2d, &
- fracday_2d, rrsun, &
- allow_negative_cosz=allow_negative_cosz, &
- half_day_out=halfday_2d)
- end if
-
-!-------------------------------------------------------------------
-!> Place output fields into 1-d arguments for return to
-!! calling routine.
-!-------------------------------------------------------------------
- fracday= fracday_2d(1,1)
- cosz = cosz_2d(1,1)
- if (present(half_day_out)) then
- half_day_out = halfday_2d(1,1)
- end if
-
-end subroutine diurnal_solar_cal_0d
-
-
-!> @brief daily_mean_solar_2d computes the daily mean astronomical parameters for
-!! the input points at latitude lat and time of year time_since_ae.
-!!
-!! @param [in] Latitudes of model grid points
-!! @param [in] Time of year; autumnal equinox = 0.0, one year = 2 * pi
-!! @param [out] Cosine of solar zenith angle
-!! @param [out] 2-d array of half-day lengths at the latitudes
-!! @param [out] the inverse of the square of the earth-sun distance relative
-!! to the mean distance at angle ang in the earth's orbit.
-!!
-!! @throw FATAL, "astronomy_mod time_since_ae not between 0 and 2pi"
-subroutine daily_mean_solar_2d (lat, time_since_ae, cosz, h_out, rr_out)
-
-!----------------------------------------------------------------------
-real, dimension(:,:), intent(in) :: lat
-real, intent(in) :: time_since_ae
-real, dimension(:,:), intent(out) :: cosz, h_out
-real, intent(out) :: rr_out
-!----------------------------------------------------------------------
-
-!--------------------------------------------------------------------
-! local variables
-!--------------------------------------------------------------------
- real, dimension(size(lat,1),size(lat,2)) :: h
- real :: ang, dec, rr
-
-!--------------------------------------------------------------------
-! be sure the time in the annual cycle is legitimate.
-!---------------------------------------------------------------------
- if (time_since_ae < 0.0 .or. time_since_ae > twopi) &
- call error_mesg('astronomy_mod', &
- 'time_since_ae not between 0 and 2pi', FATAL)
-
-!---------------------------------------------------------------------
-!> Define the orbital angle (location in year), solar declination,
-!! half-day length and earth sun distance factor. Use functions
-!! contained in this module.
-!---------------------------------------------------------------------
- ang = angle (time_since_ae)
- dec = declination(ang)
- h = half_day (lat, dec)
- rr = r_inv_squared (ang)
-
-!---------------------------------------------------------------------
-!> Where the entire day is dark, define cosz to be zero. otherwise
-!! use the standard formula. Define the daylight fraction and earth-
-!! sun distance.
-!---------------------------------------------------------------------
- where (h == 0.0)
- cosz = 0.0
- elsewhere
- cosz = sin(lat)*sin(dec) + cos(lat)*cos(dec)*sin(h)/h
- end where
- h_out = h/PI
- rr_out = rr
-
-end subroutine daily_mean_solar_2d
-
-
-!> @brief daily_mean_solar_1d takes 1-d input fields, adds a second dimension
-!! and calls daily_mean_solar_2d. on return, the 2d fields are
-!! returned to the original 1d fields.
-!!
-!! @param [in] Latitudes of model grid points
-!! @param [in] Time of year; autumnal equinox = 0.0, one year = 2 * pi
-!! @param [out] Cosine of solar zenith angle
-!! @param [out] 2-d array of half-day lengths at the latitudes
-!! @param [out] the inverse of the square of the earth-sun distance relative
-!! to the mean distance at angle ang in the earth's orbit.
-subroutine daily_mean_solar_1d (lat, time_since_ae, cosz, h_out, rr_out)
-
-!----------------------------------------------------------------------
-real, intent(in), dimension(:) :: lat
-real, intent(in) :: time_since_ae
-real, intent(out), dimension(size(lat(:))) :: cosz
-real, intent(out), dimension(size(lat(:))) :: h_out
-real, intent(out) :: rr_out
-!----------------------------------------------------------------------
-
-!----------------------------------------------------------------------
-! local variables
-!----------------------------------------------------------------------
- real, dimension(size(lat),1) :: lat_2d, cosz_2d, hout_2d
-
-!--------------------------------------------------------------------
-!> define 2-d versions of input data array.
-!--------------------------------------------------------------------
- lat_2d(:,1) = lat
-
-!--------------------------------------------------------------------
-!> call daily_mean_solar_2d to calculate astronomy fields.
-!--------------------------------------------------------------------
- call daily_mean_solar_2d (lat_2d, time_since_ae, cosz_2d, &
- hout_2d, rr_out)
-
-!-------------------------------------------------------------------
-!> place output fields into 1-d arguments for return to
-!! calling routine.
-!-------------------------------------------------------------------
- h_out = hout_2d(:,1)
- cosz = cosz_2d(:,1)
-
-end subroutine daily_mean_solar_1d
-
-
-!> @brief daily_mean_solar_2level takes 1-d input fields, adds a second
-!! dimension and calls daily_mean_solar_2d. on return, the 2d fields
-!! are returned to the original 1d fields.
-!!
-!! @param [in] Latitudes of model grid points
-!! @param [in] Time of year; autumnal equinox = 0.0, one year = 2 * pi
-!! @param [out] Cosine of solar zenith angle
-!! @param [out] Shortwave flux factor: cosine of zenith angle * daylight fraction /
-!! (earth-sun distance squared)
-subroutine daily_mean_solar_2level (lat, time_since_ae, cosz, solar)
-
-!----------------------------------------------------------------------
-real, intent(in), dimension(:) :: lat
-real, intent(in) :: time_since_ae
-real, intent(out), dimension(size(lat(:))) :: cosz, solar
-!----------------------------------------------------------------------
-
-!----------------------------------------------------------------------
-! local variables
-!----------------------------------------------------------------------
- real, dimension(size(lat),1) :: lat_2d, cosz_2d, hout_2d
- real :: rr_out
-
-!--------------------------------------------------------------------
-!> define 2-d versions of input data array.
-!--------------------------------------------------------------------
- lat_2d(:,1) = lat
-
-!--------------------------------------------------------------------
-!> call daily_mean_solar_2d to calculate astronomy fields.
-!--------------------------------------------------------------------
- call daily_mean_solar_2d (lat_2d, time_since_ae, cosz_2d, &
- hout_2d, rr_out)
-
-!-------------------------------------------------------------------
-!> place output fields into 1-d arguments for return to
-!! calling routine.
-!-------------------------------------------------------------------
- solar = cosz_2d(:,1)*hout_2d(:,1)*rr_out
- cosz = cosz_2d(:,1)
-
-end subroutine daily_mean_solar_2level
-
-
-!> @brief daily_mean_solar_1d takes 1-d input fields, adds a second dimension
-!! and calls daily_mean_solar_2d. on return, the 2d fields are
-!! returned to the original 1d fields.
-!!
-!! @param [in] Latitudes of model grid points
-!! @param [in] Time of year; autumnal equinox = 0.0, one year = 2 * pi
-!! @param [out] Cosine of solar zenith angle
-!! @param [out] 2-d array of half-day lengths at the latitudes
-!! @param [out] the inverse of the square of the earth-sun distance relative to
-!! the mean distance at angle ang in the earth's orbit.
-subroutine daily_mean_solar_0d (lat, time_since_ae, cosz, h_out, rr_out)
-
-real, intent(in) :: lat, time_since_ae
-real, intent(out) :: cosz, h_out, rr_out
-
-!--------------------------------------------------------------------
-! local variables
-!--------------------------------------------------------------------
- real, dimension(1,1) :: lat_2d, cosz_2d, hout_2d
-
-!--------------------------------------------------------------------
-!> define 2-d versions of input data array.
-!--------------------------------------------------------------------
- lat_2d = lat
-
-!--------------------------------------------------------------------
-!> call daily_mean_solar_2d to calculate astronomy fields.
-!--------------------------------------------------------------------
- call daily_mean_solar_2d (lat_2d, time_since_ae, cosz_2d, &
- hout_2d, rr_out)
-
-!-------------------------------------------------------------------
-!> return output fields to scalars for return to calling routine.
-!-------------------------------------------------------------------
- h_out = hout_2d(1,1)
- cosz = cosz_2d(1,1)
-
-end subroutine daily_mean_solar_0d
-
-
-!> @brief daily_mean_solar_cal_2d receives time_type inputs, converts
-!! them to real variables and then calls daily_mean_solar_2d to
-!! compute desired astronomical variables.
-!!
-!! @param [in] Latitudes of model grid points
-!! @param [in] Time of year (time_type)
-!! @param [out] Cosine of solar zenith angle
-!! @param [out] Daylight fraction of time interval
-!! @param [out] Earth-Sun distance (r) relative to semi-major axis of orbital ellipse (a):(a/r)**2
-!!
-!! @throw FATAL, "astronomy_mod time_since_ae not between 0 and 2pi"
-subroutine daily_mean_solar_cal_2d (lat, time, cosz, fracday, rrsun)
-
-!-------------------------------------------------------------------
-real, dimension(:,:), intent(in) :: lat
-type(time_type), intent(in) :: time
-real, dimension(:,:), intent(out) :: cosz, fracday
-real, intent(out) :: rrsun
-!-------------------------------------------------------------------
-
-!-------------------------------------------------------------------
-! local variables
-!-------------------------------------------------------------------
- real :: time_since_ae
-
-!--------------------------------------------------------------------
-! be sure the time in the annual cycle is legitimate.
-!---------------------------------------------------------------------
- time_since_ae = orbital_time(time)
- if (time_since_ae < 0.0 .or. time_since_ae > twopi) &
- call error_mesg ('astronomy_mod', &
- 'time_since_ae not between 0 and 2pi', FATAL)
-
-!--------------------------------------------------------------------
-! call daily_mean_solar_2d to calculate astronomy fields.
-!--------------------------------------------------------------------
- call daily_mean_solar_2d (lat, time_since_ae, cosz, &
- fracday, rrsun)
-
-end subroutine daily_mean_solar_cal_2d
-
-
-!> @brief daily_mean_solar_cal_1d receives time_type inputs, converts
-!! them to real, 2d variables and then calls daily_mean_solar_2d to
-!! compute desired astronomical variables.
-!!
-!! @param [in] Latitudes of model grid points
-!! @param [in] Time of year (time_type)
-!! @param [out] Cosine of solar zenith angle
-!! @param [out] Daylight fraction of time interval
-!! @param [out] Earth-Sun distance (r) relative to semi-major axis of orbital ellipse (a):(a/r)**2
-subroutine daily_mean_solar_cal_1d (lat, time, cosz, fracday, rrsun)
-
-real, dimension(:), intent(in) :: lat
-type(time_type), intent(in) :: time
-real, dimension(:), intent(out) :: cosz, fracday
-real, intent(out) :: rrsun
-
-!---------------------------------------------------------------------
-! local variables
-!---------------------------------------------------------------------
- real, dimension(size(lat),1) :: lat_2d, cosz_2d, fracday_2d
-
-
-!--------------------------------------------------------------------
-!> define 2-d versions of input data array.
-!--------------------------------------------------------------------
- lat_2d(:,1) = lat
-
-!--------------------------------------------------------------------
-!> call daily_mean_solar_cal_2d to convert the time_types to reals and
-!! then calculate the astronomy fields.
-!--------------------------------------------------------------------
- call daily_mean_solar_cal_2d (lat_2d, time, cosz_2d, &
- fracday_2d, rrsun)
-
-!-------------------------------------------------------------------
-!> place output fields into 1-d arguments for return to
-!! calling routine.
-!-------------------------------------------------------------------
- fracday = fracday_2d(:,1)
- cosz = cosz_2d(:,1)
-
-end subroutine daily_mean_solar_cal_1d
-
-
-!> @brief daily_mean_solar_cal_2level receives 1d arrays and time_type input,
-!! converts them to real, 2d variables and then calls
-!! daily_mean_solar_2d to compute desired astronomical variables.
-!!
-!! @param [in] Latitudes of model grid points
-!! @param [in] Time of year (time_type)
-!! @param [out] Cosine of solar zenith angle
-!! @param [out] Shortwave flux factor: cosine of zenith angle * daylight fraction /
-!! (earth-sun distance squared)
-subroutine daily_mean_solar_cal_2level (lat, time, cosz, solar)
-
-real, dimension(:), intent(in) :: lat
-type(time_type), intent(in) :: time
-real, dimension(:), intent(out) :: cosz, solar
-
-!---------------------------------------------------------------------
-! local variables
-!---------------------------------------------------------------------
- real, dimension(size(lat),1) :: lat_2d, cosz_2d, fracday_2d
- real :: rrsun
-
-
-!--------------------------------------------------------------------
-!> define 2-d versions of input data array.
-!--------------------------------------------------------------------
- lat_2d(:,1) = lat
-
-!--------------------------------------------------------------------
-!> call daily_mean_solar_cal_2d to convert the time_types to reals and
-!! then calculate the astronomy fields.
-!--------------------------------------------------------------------
- call daily_mean_solar_cal_2d (lat_2d, time, cosz_2d, &
- fracday_2d, rrsun)
-
-!-------------------------------------------------------------------
-!> place output fields into 1-d arguments for return to
-!! calling routine.
-!-------------------------------------------------------------------
- solar = cosz_2d(:,1)*fracday_2d(:,1)*rrsun
- cosz = cosz_2d(:,1)
-
-end subroutine daily_mean_solar_cal_2level
-
-
-!> @brief daily_mean_solar_cal_0d converts scalar input fields to real, 2d variables and
-!! then calls daily_mean_solar_2d to compute desired astronomical variables.
-!!
-!! @param [in] Latitudes of model grid points
-!! @param [in] Time of year (time_type)
-!! @param [out] Cosine of solar zenith angle
-!! @param [out] Daylight fraction of time interval
-!! @param [out] Earth-Sun distance (r) relative to semi-major axis of orbital ellipse (a):(a/r)**2
-subroutine daily_mean_solar_cal_0d (lat, time, cosz, fracday, rrsun)
-
-!--------------------------------------------------------------------
-real, intent(in) :: lat
-type(time_type), intent(in) :: time
-real, intent(out) :: cosz, fracday, rrsun
-!--------------------------------------------------------------------
-
-!--------------------------------------------------------------------
-! local variables
-!--------------------------------------------------------------------
- real, dimension(1,1) :: lat_2d, cosz_2d, fracday_2d
-
-!--------------------------------------------------------------------
-!> define 2-d versions of input data array.
-!--------------------------------------------------------------------
- lat_2d = lat
-
-!--------------------------------------------------------------------
-!> call daily_mean_solar_cal_2d to convert the time_types to reals and
-!! then calculate the astronomy fields.
-!--------------------------------------------------------------------
- call daily_mean_solar_cal_2d (lat_2d, time, cosz_2d, &
- fracday_2d, rrsun)
-
-!-------------------------------------------------------------------
-!> place output fields into scalar arguments for return to
-!! calling routine.
-!-------------------------------------------------------------------
- fracday = fracday_2d(1,1)
- cosz = cosz_2d(1,1)
-
-end subroutine daily_mean_solar_cal_0d
-
-
-!> @brief annual_mean_solar_2d returns 2d fields of annual mean values of the cosine of
-!! zenith angle, daylight fraction and earth-sun distance at the specified latitude.
-!!
-!! @param [in] Starting index of latitude window
-!! @param [in] Ending index of latitude window
-!! @param [in] Latitudes of model grid points
-!! @param [out] Cosine of solar zenith angle
-!! @param [out] Shortwave flux factor: cosine of zenith angle * daylight fraction /
-!! (earth-sun distance squared)
-!! @param [out] Daylight fraction of time interval
-!! @param [out] Earth-Sun distance (r) relative to semi-major axis of orbital ellipse (a):(a/r)**2
-subroutine annual_mean_solar_2d (js, je, lat, cosz, solar, fracday, &
- rrsun)
-
-!--------------------------------------------------------------------
-integer, intent(in) :: js, je
-real, dimension(:,:), intent(in) :: lat
-real, dimension(:,:), intent(out) :: solar, cosz, fracday
-real, intent(out) :: rrsun
-!--------------------------------------------------------------------
-
-!--------------------------------------------------------------------
-! local variables
-!--------------------------------------------------------------------
- real, dimension(size(lat,1),size(lat,2)) :: s,z
- real :: t
- integer :: n
-
-!--------------------------------------------------------------------
-! if the calculation has not yet been done, do it here.
-!--------------------------------------------------------------------
- if (.not. annual_mean_calculated) then
-
-!----------------------------------------------------------------------
-!> determine annual mean values of solar flux and product of cosz
-!! and solar flux by integrating the annual cycle in num_angles
-!! orbital increments.
-!----------------------------------------------------------------------
- solar = 0.0
- cosz = 0.0
- do n =1, num_angles
- t = float(n-1)*twopi/float(num_angles)
- call daily_mean_solar(lat,t, z, fracday, rrsun)
- s = z*rrsun*fracday
- solar = solar + s
- cosz = cosz + z*s
- end do
- solar = solar/float(num_angles)
- cosz = cosz/float(num_angles)
-
-!--------------------------------------------------------------------
-!> define the flux-weighted annual mean cosine of the zenith angle.
-!--------------------------------------------------------------------
- where(solar.eq.0.0)
- cosz = 0.0
- elsewhere
- cosz = cosz/solar
- end where
-
-!-------------------------------------------------------------------
-!> define avg fracday such as to make the avg flux (solar) equal to
-!! the product of the avg cosz * avg fracday * assumed mean avg
-!! radius of 1.0. it is unlikely that these avg fracday and avg rr
-!! will ever be used.
-!--------------------------------------------------------------------
- where(solar .eq.0.0)
- fracday = 0.0
- elsewhere
- fracday = solar/cosz
- end where
- rrsun = 1.00
-
-!---------------------------------------------------------------------
-!> save the values that have been calculated as module variables, if
-!! those variables are present; i.e., not the spectral 2-layer model.
-!---------------------------------------------------------------------
- if (allocated (cosz_ann)) then
- cosz_ann = cosz
- solar_ann = solar
- fracday_ann = fracday
- rrsun_ann = rrsun
-
-!--------------------------------------------------------------------
-!> increment the points computed counter. set flag to end execution
-!! once values have been calculated for all points owned by the
-!! processor.
-!--------------------------------------------------------------------
- num_pts = num_pts + size(lat,1)*size(lat,2)
- if ( num_pts == total_pts) annual_mean_calculated = .true.
- endif
-
-!--------------------------------------------------------------------
-!> if the calculation has been done, return the appropriate module
-!! variables.
-!--------------------------------------------------------------------
- else
- if (allocated (cosz_ann)) then
- cosz = cosz_ann
- solar = solar_ann
- fracday = fracday_ann
- rrsun = rrsun_ann
- endif
- endif
-
-end subroutine annual_mean_solar_2d
-
-
-!> @brief annual_mean_solar_1d creates 2-d input fields from 1-d input fields and then calls
-!! annual_mean_solar_2d to obtain 2-d output fields which are then stored into 1-d
-!! fields for return to the calling subroutine.
-!!
-!! @param [in] Starting index of latitude window
-!! @param [in] Ending index of latitude window
-!! @param [in] Latitudes of model grid points
-!! @param [out] Cosine of solar zenith angle
-!! @param [out] Shortwave flux factor: cosine of zenith angle * daylight fraction /
-!! (earth-sun distance squared)
-!! @param [out] Daylight fraction of time interval
-!! @param [out] Earth-Sun distance (r) relative to semi-major axis of orbital ellipse (a):(a/r)**2
-subroutine annual_mean_solar_1d (jst, jnd, lat, cosz, solar, &
- fracday, rrsun_out)
-
-!---------------------------------------------------------------------
-integer, intent(in) :: jst, jnd
-real, dimension(:), intent(in) :: lat(:)
-real, dimension(:), intent(out) :: cosz, solar, fracday
-real, intent(out) :: rrsun_out
-!---------------------------------------------------------------------
-
-!---------------------------------------------------------------------
-! local variables
-
- real, dimension(size(lat),1) :: lat_2d, solar_2d, cosz_2d, &
- fracday_2d
- real :: rrsun
-
-!--------------------------------------------------------------------
-! if the calculation has not been done, do it here.
-!--------------------------------------------------------------------
- if ( .not. annual_mean_calculated) then
-
-!--------------------------------------------------------------------
-!> define 2-d versions of input data array.
-!--------------------------------------------------------------------
- lat_2d(:,1) = lat
-
-!--------------------------------------------------------------------
-!> call annual_mean_solar_2d to calculate the astronomy fields.
-!--------------------------------------------------------------------
- call annual_mean_solar_2d (jst, jnd, lat_2d, cosz_2d, &
- solar_2d, fracday_2d, rrsun)
-
-!-------------------------------------------------------------------
-!> place output fields into 1-D arrays for return to calling routine.
-!-------------------------------------------------------------------
- fracday = fracday_2d(:,1)
- rrsun_out = rrsun
- solar = solar_2d(:,1)
- cosz = cosz_2d(:,1)
-
-!--------------------------------------------------------------------
-!> if the calculation has been done, simply return the module
-!! variables contain the results at the desired latitudes.
-!--------------------------------------------------------------------
- else
- cosz(:) = cosz_ann(1,jst:jnd)
- solar(:) = solar_ann(1,jst:jnd)
- fracday(:) = fracday_ann(1,jst:jnd)
- rrsun = rrsun_ann
- endif
-
-end subroutine annual_mean_solar_1d
-
-
-!> @brief annual_mean_solar_2level creates 2-d input fields from 1-d input fields
-!! and then calls annual_mean_solar_2d to obtain 2-d output fields which are
-!! then stored into 1-d fields for return to the calling subroutine. This
-!! subroutine will be called during model initialization.
-!!
-!! @throw FATAL, "astronomy_mod annual_mean_solar_2level should be called only once"
-subroutine annual_mean_solar_2level (lat, cosz, solar)
-
-!---------------------------------------------------------------------
-real, dimension(:), intent(in) :: lat !< Latitudes of model grid points
-real, dimension(:), intent(out) :: cosz !< Cosine of solar zenith angle
-real, dimension(:), intent(out) :: solar !< shortwave flux factor: cosine of zenith angle *
- !! daylight fraction / (earth-sun distance squared)
-!---------------------------------------------------------------------
-
-!---------------------------------------------------------------------
-! local variables
-
- real, dimension(size(lat),1) :: lat_2d, solar_2d, cosz_2d, &
- fracday_2d
- integer :: jst, jnd
- real :: rrsun
-
-!--------------------------------------------------------------------
-! if the calculation has not been done, do it here.
-!--------------------------------------------------------------------
- if ( .not. annual_mean_calculated) then
-
-!--------------------------------------------------------------------
-!> define 2-d versions of input data array.
-!--------------------------------------------------------------------
- lat_2d(:,1) = lat
- jst = 1
- jnd = size(lat(:))
-
-!--------------------------------------------------------------------
-!> call annual_mean_solar_2d to calculate the astronomy fields.
-!--------------------------------------------------------------------
- call annual_mean_solar_2d (jst, jnd, lat_2d, cosz_2d, &
- solar_2d, fracday_2d, rrsun)
-
-!-------------------------------------------------------------------
-!> place output fields into 1-D arrays for return to calling routine.
-!-------------------------------------------------------------------
- solar = solar_2d(:,1)
- cosz = cosz_2d(:,1)
-
-!--------------------------------------------------------------------
-!> if the calculation has been done, print an error message since
-!! this subroutine should be called only once.
-!--------------------------------------------------------------------
- else
- call error_mesg ('astronomy_mod', &
- 'annual_mean_solar_2level should be called only once', &
- FATAL)
- endif
- annual_mean_calculated = .true.
-
-end subroutine annual_mean_solar_2level
-
-
-!> @brief astronomy_end is the destructor for astronomy_mod.
-subroutine astronomy_end
-
-!----------------------------------------------------------------------
-!> check if the module has been initialized.
-!----------------------------------------------------------------------
- if (.not. module_is_initialized) return
-! call error_mesg ( 'astronomy_mod', &
-! ' module has not been initialized', FATAL)
-
-!----------------------------------------------------------------------
-!> deallocate module variables.
-!----------------------------------------------------------------------
- deallocate (orb_angle)
- if (allocated(cosz_ann) ) then
- deallocate (cosz_ann)
- deallocate (fracday_ann)
- deallocate (solar_ann)
- endif
-
-!----------------------------------------------------------------------
-!> mark the module as uninitialized.
-!----------------------------------------------------------------------
- module_is_initialized = .false.
-
-end subroutine astronomy_end
-
-
-!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-!
-! PRIVATE SUBROUTINES
-!
-!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-!> @brief Orbit computes and stores a table of value of orbital angles as a
-!! function of orbital time (both the angle and time are zero at
-!! autumnal equinox in the NH, and range from 0 to 2*pi).
-subroutine orbit
-
-!---------------------------------------------------------------------
-! local variables
-
- integer :: n
- real :: d1, d2, d3, d4, d5, dt, norm
-
-!--------------------------------------------------------------------
-!> allocate the orbital angle array, sized by the namelist parameter
-!! num_angles, defining the annual cycle resolution of the earth's
-!! orbit. define some constants to be used.
-!--------------------------------------------------------------------
-! wfc moving to astronomy_init
-! allocate ( orb_angle(0:num_angles) )
- orb_angle(0) = 0.0
- dt = twopi/float(num_angles)
- norm = sqrt(1.0 - ecc**2)
- dt = dt*norm
-
-!---------------------------------------------------------------------
-!> define the orbital angle at each of the num_angles locations in
-!! the orbit.
-!---------------------------------------------------------------------
- do n = 1, num_angles
- d1 = dt*r_inv_squared(orb_angle(n-1))
- d2 = dt*r_inv_squared(orb_angle(n-1)+0.5*d1)
- d3 = dt*r_inv_squared(orb_angle(n-1)+0.5*d2)
- d4 = dt*r_inv_squared(orb_angle(n-1)+d3)
- d5 = d1/6.0 + d2/3.0 +d3/3.0 +d4/6.0
- orb_angle(n) = orb_angle(n-1) + d5
- end do
-
-end subroutine orbit
-
-
-!> @brief r_inv_squared returns the inverse of the square of the earth-sun
-!! distance relative to the mean distance at angle ang in the Earth's orbit.
-function r_inv_squared (ang)
-
-!--------------------------------------------------------------------
-real, intent(in) :: ang !< angular position of earth in its orbit, relative to a
- !! value of 0.0 at the NH autumnal equinox, value between
- !! 0.0 and 2 * pi [radians]
-!--------------------------------------------------------------------
-
-!---------------------------------------------------------------------
-! local variables
-
-real :: r_inv_squared !< The inverse of the square of the earth-sun distance relative
- !! to the mean distance [dimensionless]
-real :: r !< Earth-Sun distance relative to mean distance [dimensionless]
-real :: rad_per !< Angular position of perihelion [radians]
-
-!--------------------------------------------------------------------
-!> define the earth-sun distance (r) and then return the inverse of
-!! its square (r_inv_squared) to the calling routine.
-!--------------------------------------------------------------------
- rad_per = per*deg_to_rad
- r = (1. - ecc**2)/(1. + ecc*cos(ang - rad_per))
- r_inv_squared = r**(-2)
-
-
-end function r_inv_squared
-
-
-!> @brief angle determines the position within the earth's orbit at time t
-!! in the year (t = 0 at NH autumnal equinox) by interpolating
-!! into the orbital position table.
-function angle (t)
-
-!--------------------------------------------------------------------
-real, intent(in) :: t !< time of year (between 0 and 2*pi; t=0 at NH autumnal equinox
-!--------------------------------------------------------------------
-
-!--------------------------------------------------------------------
-! local variables
-!--------------------------------------------------------------------
- real :: angle !< Orbital position relative to NH autumnal equinox [radians]
- real :: norm_time !< Index into orbital table corresponding to input time [dimensionless]
- real :: x !< Fractional distance between the orbital table entries bracketing the input time [dimensionless]
- integer :: int !< Table index which is lower than actual position, but closest to it [dimensionless]
- integer :: int_1 !< Next table index just larger than actual orbital position [dimensionless]
-
-!--------------------------------------------------------------------
-!> Define orbital tables indices bracketing current orbital time
-!! (int and int_1). Define table index distance between the lower
-!! table value (int) and the actual orbital time (x). Define orbital
-!! position as being x of the way between int and int_1. Renormalize
-!! angle to be within the range 0 to 2*pi.
-!--------------------------------------------------------------------
- norm_time = t*float(num_angles)/twopi
- int = floor(norm_time)
- int = modulo(int,num_angles)
- int_1 = int+1
- x = norm_time - floor(norm_time)
- angle = (1.0 -x)*orb_angle(int) + x*orb_angle(int_1)
- angle = modulo(angle, twopi)
-
-end function angle
-
-
-!> @brief Declination returns the solar declination angle at orbital
-!! position ang in earth's orbit.
-function declination (ang)
-
-!--------------------------------------------------------------------
-real, intent(in) :: ang !< solar orbital position ang in earth's orbit
-!--------------------------------------------------------------------
-
-!--------------------------------------------------------------------
-! local variables
-
- real :: declination !< Solar declination angle [radians]
- real :: rad_obliq !< Obliquity of the ecliptic [radians]
- real :: sin_dec !< Sine of the solar declination [dimensionless]
-
-!---------------------------------------------------------------------
-! compute the solar declination.
-!---------------------------------------------------------------------
- rad_obliq = obliq*deg_to_rad
- sin_dec = - sin(rad_obliq)*sin(ang)
- declination = asin(sin_dec)
-
-end function declination
-
-
-!> @brief half_day_2d returns a 2-d array of half-day lengths at the
-!! latitudes and declination provided.
-!!
-function half_day_2d (latitude, dec) result(h)
-
-!---------------------------------------------------------------------
-real, dimension(:,:), intent(in) :: latitude !< Latitutde of view point
-real, intent(in) :: dec !< Solar declination angle at view point
-real, dimension(size(latitude,1),size(latitude,2)) :: h
-!---------------------------------------------------------------------
-
-!---------------------------------------------------------------------
-! local variables
-!---------------------------------------------------------------------
- real, dimension (size(latitude,1),size(latitude,2)):: &
- cos_half_day, & !< Cosine of half-day length [dimensionless]
- lat !< Model latitude, adjusted so that it is never 0.5*pi or -0.5*pi
- real :: tan_dec !< tangent of solar declination [dimensionless]
- real :: eps = 1.0E-05 !< small increment
-
-!--------------------------------------------------------------------
-!> define tangent of the declination.
-!--------------------------------------------------------------------
- tan_dec = tan(dec)
-
-!--------------------------------------------------------------------
-!> adjust latitude so that its tangent will be defined.
-!--------------------------------------------------------------------
- lat = latitude
- where (latitude == 0.5*PI) lat= latitude - eps
- where (latitude == -0.5*PI) lat= latitude + eps
-
-!--------------------------------------------------------------------
-!> define the cosine of the half-day length. adjust for cases of
-!! all daylight or all night.
-!--------------------------------------------------------------------
- cos_half_day = -tan(lat)*tan_dec
- where (cos_half_day <= -1.0) h = PI
- where (cos_half_day >= +1.0) h = 0.0
- where(cos_half_day > -1.0 .and. cos_half_day < 1.0) &
- h = acos(cos_half_day)
-
-end function half_day_2d
-
-
-!> @brief half_day_0d takes scalar input fields, makes them into 2-d fields
-!! dimensioned (1,1), and calls half_day_2d. On return, the 2-d
-!! fields are converted to the desired scalar output.
-!!
-!! @param [in] Latitutde of view point
-!! @param [in] Solar declination angle at view point
-function half_day_0d(latitude, dec) result(h)
-
-real, intent(in) :: latitude, dec
-real :: h
-
-!----------------------------------------------------------------------
-! local variables
-!----------------------------------------------------------------------
- real, dimension(1,1) :: lat_2d, h_2d
-
-!---------------------------------------------------------------------
-! create 2d array from the input latitude field.
-!---------------------------------------------------------------------
- lat_2d = latitude
-
-!---------------------------------------------------------------------
-! call half_day with the 2d arguments to calculate half-day length.
-!---------------------------------------------------------------------
- h_2d = half_day (lat_2d, dec)
-
-!---------------------------------------------------------------------
-! create scalar from 2d array.
-!---------------------------------------------------------------------
- h = h_2d(1,1)
-
-end function half_day_0d
-
-
-!> @brief Orbital time returns the time (1 year = 2*pi) since autumnal
-!! equinox
-!!
-!! @details Orbital time returns the time (1 year = 2*pi) since autumnal
-!! equinox; autumnal_eq_ref is a module variable of time_type and
-!! will have been defined by default or by a call to
-!! set_ref_date_of_ae; length_of_year is available through the time
-!! manager and is set at the value approriate for the calandar being used
-function orbital_time(time) result(t)
-
-type(time_type), intent(in) :: time !< time (1 year = 2*pi) since autumnal equinox
-real :: t
-
- t = real ( (time - autumnal_eq_ref)//period_time_type)
- t = twopi*(t - floor(t))
- if (time < autumnal_eq_ref) t = twopi - t
-
-
-end function orbital_time
-
-
-!> @brief universal_time returns the time of day at longitude = 0.0
-!! (1 day = 2*pi)
-function universal_time(time) result(t)
-
-type(time_type), intent(in) :: time !< Time (1 year = 2*pi) since autumnal equinox
-real :: t
-
-!--------------------------------------------------------------------
-! local variables
-!--------------------------------------------------------------------
- integer :: seconds, days
-
- call get_time (time, seconds, days)
- t = twopi*real(seconds)/seconds_per_day
-
-end function universal_time
+#include "astronomy_r4.fh"
+#include "astronomy_r8.fh"
end module astronomy_mod
!> @}
diff --git a/astronomy/include/astronomy.inc b/astronomy/include/astronomy.inc
index 192fc8f22a..d26f97a7de 100644
--- a/astronomy/include/astronomy.inc
+++ b/astronomy/include/astronomy.inc
@@ -26,586 +26,6 @@
!> @addtogroup astronomy_mod
!> @{
-module astronomy_mod
-
-
-use fms_mod, only: fms_init, &
- mpp_pe, mpp_root_pe, stdlog, &
- write_version_number, &
- check_nml_error, error_mesg, &
- FATAL, NOTE, WARNING
-use time_manager_mod, only: time_type, set_time, get_time, &
- get_date_julian, set_date_julian, &
- set_date, length_of_year, &
- time_manager_init, &
- operator(-), operator(+), &
- operator( // ), operator(<)
-use constants_mod, only: constants_init, PI
-use mpp_mod, only: input_nml_file
-
-!--------------------------------------------------------------------
-
-implicit none
-private
-
-!---------------------------------------------------------------------
-!----------- version number for this module --------------------------
-
-! Include variable "version" to be written to log file.
-#include
-
-
-!---------------------------------------------------------------------
-!------- interfaces --------
-
-public &
- astronomy_init, get_period, set_period, &
- set_orbital_parameters, get_orbital_parameters, &
- set_ref_date_of_ae, get_ref_date_of_ae, &
- diurnal_solar, daily_mean_solar, annual_mean_solar, &
- astronomy_end, universal_time, orbital_time
-
-!> @}
-
-!> @brief Calculates solar information for the given location(lat & lon) and time
-!!
-!> ~~~~~~~~~~{.f90}
-!! call diurnal_solar (lat, lon, time, cosz, fracday, rrsun, dt_time)
-!! call diurnal_solar (lat, lon, gmt, time_since_ae, cosz, fracday, rrsun, dt)
-!! ~~~~~~~~~~
-!!
-!! The first option (used in conjunction with time_manager_mod)
-!! generates the real variables gmt and time_since_ae from the
-!! time_type input, and then calls diurnal_solar with these real inputs.
-!!
-!! The time of day is set by
-!! ~~~~~~~~~~{.f90}
-!! real, intent(in) :: gmt
-!! ~~~~~~~~~~
-!! The time of year is set by
-!! ~~~~~~~~~~{.f90}
-!! real, intent(in) :: time_since_ae
-!! ~~~~~~~~~~
-!! with time_type input, both of these are extracted from
-!! ~~~~~~~~~~{.f90}
-!! type(time_type), intent(in) :: time
-!! ~~~~~~~~~~
-!!
-!! Separate routines exist within this interface for scalar,
-!! 1D or 2D input and output fields:
-!!
-!! ~~~~~~~~~~{.f90}
-!! real, intent(in), dimension(:,:) :: lat, lon
-!! real, intent(in), dimension(:) :: lat, lon
-!! real, intent(in) :: lat, lon
-!!
-!! real, intent(out), dimension(:,:) :: cosz, fracday
-!! real, intent(out), dimension(:) :: cosz, fracday
-!! real, intent(out) :: cosz, fracday
-!! ~~~~~~~~~~
-!!
-!! One may also average the output fields over the time interval
-!! between gmt and gmt + dt by including the optional argument dt (or
-!! dt_time). dt is measured in radians and must be less than pi
-!! (1/2 day). This average is computed analytically, and should be
-!! exact except for the fact that changes in earth-sun distance over
-!! the time interval dt are ignored. In the context of a diurnal GCM,
-!! this option should always be employed to insure that the total flux
-!! at the top of the atmosphere is not modified by time truncation error.
-!!
-!! ~~~~~~~~~~{.f90}
-!! real, intent(in), optional :: dt
-!! type(time_type), optional :: dt_time
-!! ~~~~~~~~~~
-!!
-!! @param [in] Latitudes of model grid points [radians]
-!! @param [in] Longitudes of model grid points [radians]
-!! @param [in] Time of day at longitude 0.0; midnight = 0.0, one day = 2 * pi [radians]
-!! @param [in] Time of year; autumnal equinox = 0.0, one year = 2 * pi [radians]
-!! @param [in] Time at which astronomical values are desired (time_type variable) [seconds, days]
-!! @param [out] Cosine of solar zenith angle, set to zero when entire period is in darkness [dimensionless]
-!! @param [out] Daylight fraction of time interval [dimensionless]
-!! @param [out] Earth-Sun distance (r) relative to semi-major axis of orbital ellipse
-!! (a):(a/r)**2 [dimensionless]
-!! @param [in] OPTIONAL: Time interval after gmt over which the astronomical variables are to be
-!! averaged. this produces averaged output rather than instantaneous. [radians], (1 day = 2 * pi)
-!! @param [in] OPTIONAL: Time interval after gmt over which the astronomical variables are to be
-!! averaged. this produces averaged output rather than instantaneous. time_type, [days, seconds]
-!! @param [in] Allow negative values for cosz?
-!! @param [out] half_day_out
-!> @ingroup astronomy_mod
-interface diurnal_solar
- module procedure diurnal_solar_2d
- module procedure diurnal_solar_1d
- module procedure diurnal_solar_0d
- module procedure diurnal_solar_cal_2d
- module procedure diurnal_solar_cal_1d
- module procedure diurnal_solar_cal_0d
-end interface
-
-!> @brief Calculates the daily mean solar information for a given time and latitude.
-!!
-!> ~~~~~~~~~~{.f90}
-!! call daily_mean_solar (lat, time, cosz, fracday, rrsun)
-!! call daily_mean_solar (lat, time_since_ae, cosz, fracday, rrsun)
-!! call daily_mean_solar (lat, time, cosz, solar)
-!! call daily_mean_solar (lat, time_since_ae, cosz, solar)
-!! ~~~~~~~~~~
-!!
-!! The first option (used in conjunction with time_manager_mod)
-!! generates the real variable time_since_ae from the time_type
-!! input time, and then calls daily_mean_solar with this real input
-!! (option 2). The third and fourth options correspond to the first
-!! and second and are used with then spectral 2-layer model, where
-!! only cosz and solar are desired as output. These routines generate
-!! dummy arguments and then call option 2, where the calculation is done.
-!!
-!! The time of year is set by
-!! ~~~~~~~~~~{.f90}
-!! real, intent(in) :: time_since_ae
-!! ~~~~~~~~~~
-!! With time_type input, the time of year is extracted from
-!! ~~~~~~~~~~{.f90}
-!! type(time_type), intent(in) :: time
-!! ~~~~~~~~~~
-!!
-!! Separate routines exist within this interface for scalar,
-!! 1D or 2D input and output fields:
-!!
-!! ~~~~~~~~~~{.f90}
-!! real, intent(in), dimension(:,:) :: lat
-!! real, intent(in), dimension(:) :: lat
-!! real, intent(in) :: lat
-!!
-!! real, intent(out), dimension(:,:) :: cosz, fracday
-!! real, intent(out), dimension(:) :: cosz, fracday
-!! real, intent(out) :: cosz, fracday
-!! ~~~~~~~~~~
-!!
-!! @param [in] Latitudes of model grid points [radians]
-!! @param [in] Time of year; autumnal equinox = 0.0, one year = 2 * pi [radians]
-!! @param [in] Time at which astronomical values are desired (time_type variable) [seconds, days]
-!! @param [out] Cosine of solar zenith angle, set to zero when entire period is in darkness [dimensionless]
-!! @param [out] Daylight fraction of time interval [dimensionless]
-!! @param [out] Earth-Sun distance (r) relative to semi-major axis of orbital ellipse
-!! (a):(a/r)**2 [dimensionless]
-!! @param [out] shortwave flux factor: cosine of zenith angle * daylight fraction /
-!! (earth-sun distance squared) [dimensionless]
-!> @ingroup astronomy_mod
-interface daily_mean_solar
- module procedure daily_mean_solar_2d
- module procedure daily_mean_solar_1d
- module procedure daily_mean_solar_2level
- module procedure daily_mean_solar_0d
- module procedure daily_mean_solar_cal_2d
- module procedure daily_mean_solar_cal_1d
- module procedure daily_mean_solar_cal_2level
- module procedure daily_mean_solar_cal_0d
-end interface
-
-!> Calculates the annual mean of solar information for a given latitude and time.
-!!
-!> ~~~~~~~~~~{.f90}
-!! call annual_mean_solar (js, je, lat, cosz, solar, fracday, rrsun)
-!! call annual_mean_solar (lat, cosz, solar)
-!! ~~~~~~~~~~
-!!
-!! The second interface above is used by the spectral 2-layer model,
-!! which requires only cosz and solar as output arguments, and which
-!! makes this call during the initialization phase of the model.
-!! Separate routines exist within this interface for 1D or 2D input
-!! and output fields:
-!!
-!! ~~~~~~~~~~{.f90}
-!! real, intent(in), dimension(:,:) :: lat
-!! real, intent(in), dimension(:) :: lat
-!!
-!! real, intent(out), dimension(:,:) :: cosz, solar, fracday
-!! real, intent(out), dimension(:) :: cosz, solar, fracday
-!! ~~~~~~~~~~
-!!
-!! @param [in] Starting subdomain j indices of data in the physics wiondow being integrated
-!! @param [in] Ending subdomain j indices of data in the physics wiondow being integrated
-!! @param [in] Latitudes of model grid points [radians]
-!! @param [out] cosz is the average over the year of the cosine of an effective zenith angle
-!! that would produce the correct daily solar flux if the sun were fixed at that
-!! single position for the period of daylight on the given day. in this average,
-!! the daily mean effective cosz is weighted by the daily mean solar flux. [dimensionless]
-!! @param [out] Normalized solar flux, averaged over the year, equal to the product of
-!! fracday*cosz*rrsun [dimensionless]
-!! @param [out] Daylight fraction calculated so as to make the average flux (solar) equal to the
-!! product of the flux-weighted avg cosz * this fracday * assumed annual mean avg
-!! Earth-Sun distance of 1.0. [dimensionless]
-!! @param [out] Annual mean Earth-Sun distance (r) relative to semi-major axis of orbital ellipse
-!! (a):(a/r)**2 [dimensionless]
-!> @ingroup astronomy_mod
-interface annual_mean_solar
- module procedure annual_mean_solar_2d
- module procedure annual_mean_solar_1d
- module procedure annual_mean_solar_2level
-end interface
-
-!> Gets the length of year for current calendar
-!!
-!> ~~~~~~~~~~{.f90}
-!! call get_period (period)
-!! ~~~~~~~~~~
-!!
-!! Separate routines exist within this interface for integer
-!! and time_type output:
-!!
-!! ~~~~~~~~~~{.f90}
-!! integer, intent(out) :: period
-!! type(time_type), intent(out) :: period
-!! ~~~~~~~~~~
-!!
-!! @param [out] Length of year for calendar in use
-!> @ingroup astronomy_mod
-interface get_period
- module procedure get_period_time_type, get_period_integer
-end interface
-
-!> Sets the length of a year for the calendar in use
-!!
-!> ~~~~~~~~~~{.f90}
-!! call set_period (period_in)
-!! ~~~~~~~~~~
-!!
-!! Separate routines exist within this interface for integer
-!! and time_type output:
-!!
-!! ~~~~~~~~~~{.f90}
-!! integer, intent(out) :: period_in
-!! type(time_type), intent(out) :: period_in
-!! ~~~~~~~~~~
-!!
-!! @param [in] Length of year for calendar in use
-!> @ingroup astronomy_mod
-interface set_period
- module procedure set_period_time_type, set_period_integer
-end interface
-
-
-private &
- orbit, & ! Called from astronomy_init and set_orbital_parameters
- r_inv_squared, & ! Called from diurnal_solar, daily_mean_solar and orbit
- angle, declination, half_day ! called from diurnal_solar and daily_mean_solar
-! half_day, orbital_time, & ! called from diurnal_solar and daily_mean_solar
-! universal_time ! called from diurnal_solar:
-
-!> Private interface for internal use by dirunal_solar and daily_mean_solar.
-!!
-!> Example usage:
-!! ~~~~~~~~~~{.f90}
-!! half_day (latitude, dec) result (h)
-!! ~~~~~~~~~~
-!!
-!! Separate routines exist within this interface for scalar,
-!! or 2D input and output fields:
-!!
-!! ~~~~~~~~~~{.f90}
-!! real, intent(in), dimension(:,:) :: latitude
-!! real, intent(in) :: latitude
-!!
-!! real, dimension(size(latitude,1),size(latitude,2)) :: h
-!! real :: h
-!! ~~~~~~~~~~
-!!
-!! @param [in] Latitudes of model grid points [radians]
-!! @param [in] Solar declination [radians]
-!! @param [out] Half of the length of daylight at the given latitude and orbital position (dec); value
-!! ranges between 0 (all darkness) and pi (all daylight) [dimensionless]
-!> @ingroup astronomy_mod
-interface half_day
- module procedure half_day_2d, half_day_0d
-end interface
-
-!> @addtogroup astronomy_mod
-!> @{
-
-!---------------------------------------------------------------------
-!-------- namelist ---------
-
-real :: ecc = 0.01671 !< Eccentricity of Earth's orbit [dimensionless]
-real :: obliq = 23.439 !< Obliquity [degrees]
-real :: per = 102.932 !< Longitude of perihelion with respect
- !! to autumnal equinox in NH [degrees]
-integer :: period = 0 !< Specified length of year [seconds];
- !! must be specified to override default
- !! value given by length_of_year in
- !! time_manager_mod
-integer :: day_ae = 23 !< Day of specified autumnal equinox
-integer :: month_ae = 9 !< Month of specified autumnal equinox
-integer :: year_ae = 1998 !< Year of specified autumnal equinox
-integer :: hour_ae = 5 !< Hour of specified autumnal equinox
-integer :: minute_ae = 37 !< Minute of specified autumnal equinox
-integer :: second_ae = 0 !< Second of specified autumnal equinox
-integer :: num_angles = 3600 !< Number of intervals into which the year
- !! is divided to compute orbital positions
-
-
-namelist /astronomy_nml/ ecc, obliq, per, period, &
- year_ae, month_ae, day_ae, &
- hour_ae, minute_ae, second_ae, &
- num_angles
-
-!--------------------------------------------------------------------
-!------ public data ----------
-
-
-!--------------------------------------------------------------------
-!------ private data ----------
-
-type(time_type) :: autumnal_eq_ref !< time_type variable containing
- !! specified time of reference
- !! NH autumnal equinox
-
-type(time_type) :: period_time_type !< time_type variable containing
- !! period of one orbit
-
-real, dimension(:), allocatable :: orb_angle !< table of orbital positions (0 to
- !! 2*pi) as a function of time used
- !! to find actual orbital position
- !! via interpolation
-
-real :: seconds_per_day=86400. !< seconds in a day
-real :: deg_to_rad !< conversion from degrees to radians
-real :: twopi !< 2 *PI
-logical :: module_is_initialized=.false. !< has the module been initialized ?
-
-real, dimension(:,:), allocatable :: &
- cosz_ann, & !< annual mean cos of zenith angle
- solar_ann, & !< annual mean solar factor
- fracday_ann !< annual mean daylight fraction
-real :: rrsun_ann !< annual mean earth-sun distance
-logical :: annual_mean_calculated=.false. !< have the annual mean values been calculated?
-integer :: num_pts = 0 !< count of grid_boxes for which
- !! annual mean astronomy values have
- !! been calculated
-integer :: total_pts !< number of grid boxes owned by the processor
-
-!--------------------------------------------------------------------
-
-
-
- contains
-
-
-
-!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-!
-! PUBLIC SUBROUTINES
-!
-!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-!> @brief astronomy_init is the constructor for astronomy_mod.
-!!
-!! @throw FATAL, "astronomy_mod ecc must be between 0 and 0.99"
-!! @throw FATAL, "astronomy_mod obliquity must be between -90 and 90 degrees"
-!! @throw FATAL, "astronomy_mod perihelion must be between 0 and 360 degrees"
-subroutine astronomy_init (latb, lonb)
-
-real, dimension(:,:), intent(in), optional :: latb !< 2d array of model latitudes at cell corners [radians]
-real, dimension(:,:), intent(in), optional :: lonb !< 2d array of model longitudes at cell corners [radians]
-
-!-------------------------------------------------------------------
-! local variables:
-!-------------------------------------------------------------------
-integer :: unit, ierr, io, seconds, days, jd, id
-
-!-------------------------------------------------------------------
-! if module has already been initialized, exit.
-!-------------------------------------------------------------------
- if (module_is_initialized) return
-
-!-------------------------------------------------------------------
-!> This routine will:
-!> Verify that modules used by this module have been initialized.
-!-------------------------------------------------------------------
- call fms_init
- call time_manager_init
- call constants_init
-
-!-----------------------------------------------------------------------
-!> Read namelist.
-!-----------------------------------------------------------------------
- read (input_nml_file, astronomy_nml, iostat=io)
- ierr = check_nml_error(io,'astronomy_nml')
-!---------------------------------------------------------------------
-!> Write version number and namelist to logfile.
-!---------------------------------------------------------------------
- call write_version_number("ASTRONOMY_MOD", version)
- if (mpp_pe() == mpp_root_pe() ) then
- unit = stdlog()
- write (unit, nml=astronomy_nml)
- endif
-!--------------------------------------------------------------------
-!> Be sure input values are within valid ranges.
-! QUESTION : ARE THESE THE RIGHT LIMITS ???
-!---------------------------------------------------------------------
- if (ecc < 0.0 .or. ecc > 0.99) &
- call error_mesg ('astronomy_mod', &
- 'ecc must be between 0 and 0.99', FATAL)
- if (obliq < -90. .or. obliq > 90.) &
- call error_mesg ('astronomy_mod', &
- 'obliquity must be between -90 and 90 degrees', FATAL)
- if (per < 0.0 .or. per > 360.0) &
- call error_mesg ('astronomy_mod', &
- 'perihelion must be between 0 and 360 degrees', FATAL)
-
-!----------------------------------------------------------------------
-!> Set up time-type variable defining specified time of autumnal equinox.
-!----------------------------------------------------------------------
- autumnal_eq_ref = set_date (year_ae,month_ae,day_ae, &
- hour_ae,minute_ae,second_ae)
-
-!---------------------------------------------------------------------
-!> Set up time-type variable defining length of year.
-!----------------------------------------------------------------------
- if (period == 0) then
- period_time_type = length_of_year()
- call get_time (period_time_type, seconds, days)
- period = int(seconds_per_day*days + seconds)
- else
- period_time_type = set_time(period,0)
- endif
-
-!---------------------------------------------------------------------
-!> Define useful module variables.
-!----------------------------------------------------------------------
- twopi = 2.*PI
- deg_to_rad = twopi/360.
-
-!---------------------------------------------------------------------
-!> Call orbit to define table of orbital angles as function of
-!! orbital time.
-!----------------------------------------------------------------------
-! wfc moved here from orbit
- allocate ( orb_angle(0:num_angles) )
- call orbit
-
-!--------------------------------------------------------------------
-!> If annual mean radiation is desired, then latb will be present.
-!! allocate arrays to hold the needed astronomical factors. define
-!! the total number of points that the processor is responsible for.
-!--------------------------------------------------------------------
- if (present(latb)) then
- jd = size(latb,2) - 1
- id = size(lonb,1) - 1
- allocate (cosz_ann(id, jd))
- allocate (solar_ann(id, jd))
- allocate (fracday_ann(id, jd))
- total_pts = jd*id
- endif
-
-!---------------------------------------------------------------------
-!> Mark the module as initialized.
-!---------------------------------------------------------------------
- module_is_initialized=.true.
-
-!---------------------------------------------------------------------
-
-end subroutine astronomy_init
-
-
-!> @brief get_period_integer returns the length of the year as an
-!! integer number of seconds.
-!!
-!! @throw FATAL, "astronomy_mod module has not been initialized"
-subroutine get_period_integer (period_out)
-
-integer, intent(out) :: period_out !< Length of year [seconds]
-
-!--------------------------------------------------------------------
-! local variables:
-
- integer :: seconds, days
-
-!---------------------------------------------------------------------
-! exit if module has not been initialized.
-!---------------------------------------------------------------------
- if (.not. module_is_initialized) &
- call error_mesg ( 'astronomy_mod', &
- ' module has not been initialized', FATAL)
-
-!--------------------------------------------------------------------
-! define length of year in seconds.
-!--------------------------------------------------------------------
- call get_time (period_time_type, seconds, days)
- period_out = int(seconds_per_day*days + seconds)
-
-
-end subroutine get_period_integer
-
-!> @brief get_period_time_type returns the length of the year as a time_type
-!! variable.
-!!
-!! @throw FATAL, "astronomy_mod module has not been initialized"
-subroutine get_period_time_type (period_out)
-
-type(time_type), intent(inout) :: period_out !< Length of year as time_type variable
-
-!---------------------------------------------------------------------
-! exit if module has not been initialized.
-!---------------------------------------------------------------------
- if (.not. module_is_initialized) &
- call error_mesg ( 'astronomy_mod', &
- ' module has not been initialized', FATAL)
-
-!--------------------------------------------------------------------
-! define length of year as a time_type variable.
-!--------------------------------------------------------------------
- period_out = period_time_type
-
-end subroutine get_period_time_type
-
-!> @brief set_period_integer saves as the input length of the year (an
-!! integer) in a time_type module variable.
-!!
-!! @throw FATAL, "astronomy_mod module has not been initialized"
-subroutine set_period_integer (period_in)
-
-integer, intent(in) :: period_in !< Length of year as a time_type
-
-!---------------------------------------------------------------------
-! exit if module has not been initialized.
-!---------------------------------------------------------------------
- if (.not. module_is_initialized) &
- call error_mesg ( 'astronomy_mod', &
- ' module has not been initialized', FATAL)
-
-!---------------------------------------------------------------------
-! define time_type variable defining the length of year from the
-! input value (integer seconds).
-!---------------------------------------------------------------------
- period_time_type = set_time(period_in, 0)
-
-end subroutine set_period_integer
-
-
-!> @brief Set_period_time_type saves the length of the year (input as a
-!! time_type variable) into a time_type module variable.
-!!
-!! @throw FATAL, "astronomy_mod module has not been initialized"
-subroutine set_period_time_type(period_in)
-
-type(time_type), intent(in) :: period_in
-
-!---------------------------------------------------------------------
-! exit if module has not been initialized.
-!---------------------------------------------------------------------
- if (.not. module_is_initialized) &
- call error_mesg ( 'astronomy_mod', &
- ' module has not been initialized', FATAL)
-
-!---------------------------------------------------------------------
-! define time_type variable defining the length of year from the
-! input value (time_type).
-!---------------------------------------------------------------------
- period_time_type = period_in
-
-
-end subroutine set_period_time_type
!> @brief set_orbital_parameters saves the input values of eccentricity,
!! obliquity and perihelion time as module variables for use by
@@ -615,181 +35,88 @@ end subroutine set_period_time_type
!! @throw FATAL, "astronomy_mod ecc must be between 0 and 0.99"
!! @throw FATAL, "astronomy_mod obliquity must be between -90. and 90. degrees"
!! @throw FATAL, "astronomy_mod perihelion must be between 0.0 and 360. degrees"
-subroutine set_orbital_parameters (ecc_in, obliq_in, per_in)
-real, intent(in) :: ecc_in !< Eccentricity of orbital ellipse [dimensionless]
-real, intent(in) :: obliq_in !< Obliquity [degrees]
-real, intent(in) :: per_in !< Longitude of perihelion with respect to autumnal
+subroutine SET_ORBITAL_PARAMETERS_(ecc_in, obliq_in, per_in)
+
+real(kind=FMS_AST_KIND_), intent(in) :: ecc_in !< Eccentricity of orbital ellipse [dimensionless]
+real(kind=FMS_AST_KIND_), intent(in) :: obliq_in !< Obliquity [degrees]
+real(kind=FMS_AST_KIND_), intent(in) :: per_in !< Longitude of perihelion with respect to autumnal
!! equinox in northern hemisphere [degrees]
+integer, parameter :: lkind = FMS_AST_KIND_
+
!---------------------------------------------------------------------
! exit if module has not been initialized.
!---------------------------------------------------------------------
- if (.not. module_is_initialized) &
- call error_mesg ( 'astronomy_mod', &
- ' module has not been initialized', FATAL)
+ if (.not. module_is_initialized) &
+ call error_mesg('astronomy_mod', 'module has not been initialized', FATAL)
!--------------------------------------------------------------------
! be sure input values are within valid ranges.
! QUESTION : ARE THESE THE RIGHT LIMITS ???
!---------------------------------------------------------------------
- if (ecc_in < 0.0 .or. ecc_in > 0.99) &
- call error_mesg ('astronomy_mod', &
- 'ecc must be between 0 and 0.99', FATAL)
- if (obliq_in < -90.0 .or. obliq > 90.0) &
- call error_mesg ('astronomy_mod', &
- 'obliquity must be between -90. and 90. degrees', FATAL)
- if (per_in < 0.0 .or. per_in > 360.0) &
- call error_mesg ('astronomy_mod', &
- 'perihelion must be between 0.0 and 360. degrees', FATAL)
+
+ if (ecc_in < 0.0_lkind .or. ecc_in > 0.99_lkind) &
+ call error_mesg('astronomy_mod', 'ecc must be between 0 and 0.99', FATAL)
+ if (obliq_in < -90.0_lkind .or. real(obliq, FMS_AST_KIND_) > 90.0_lkind) &
+ call error_mesg('astronomy_mod', 'obliquity must be between -90. and 90. degrees', FATAL)
+ if (per_in < 0.0_lkind .or. per_in > 360.0_lkind) &
+ call error_mesg('astronomy_mod', 'perihelion must be between 0.0 and 360. degrees', FATAL)
!---------------------------------------------------------------------
! save input values into module variables.
!---------------------------------------------------------------------
- ecc = ecc_in
- obliq = obliq_in
- per = per_in
+
+ ecc = real(ecc_in, r8_kind)
+ obliq = real(obliq_in, r8_kind)
+ per = real(per_in, r8_kind)
!---------------------------------------------------------------------
! call orbit to define table of orbital angles as function of
! orbital time using the input values of parameters just supplied.
!----------------------------------------------------------------------
- call orbit
+
+ call orbit
!----------------------------------------------------------------------
-end subroutine set_orbital_parameters
+end subroutine SET_ORBITAL_PARAMETERS_
!> @brief get_orbital_parameters retrieves the orbital parameters for use
!! by another module.
!!
!! @throw FATAL, "astronomy_mod module has not been initialized"
-subroutine get_orbital_parameters (ecc_out, obliq_out, per_out)
+
+subroutine GET_ORBITAL_PARAMETERS_(ecc_out, obliq_out, per_out)
!-------------------------------------------------------------------
! get_orbital_parameters retrieves the orbital parameters for use
! by another module.
!--------------------------------------------------------------------
-real, intent(out) :: ecc_out !< Eccentricity of orbital ellipse [dimensionless]
-real, intent(out) :: obliq_out !< Obliquity [degrees]
-real, intent(out) :: per_out !< Longitude of perihelion with respect to autumnal
+
+real(kind=FMS_AST_KIND_), intent(out) :: ecc_out !< Eccentricity of orbital ellipse [dimensionless]
+real(kind=FMS_AST_KIND_), intent(out) :: obliq_out !< Obliquity [degrees]
+real(kind=FMS_AST_KIND_), intent(out) :: per_out !< Longitude of perihelion with respect to autumnal
!! equinox in northern hemisphere [degrees]
!---------------------------------------------------------------------
! exit if module has not been initialized.
!---------------------------------------------------------------------
- if (.not. module_is_initialized) &
- call error_mesg ( 'astronomy_mod', &
- ' module has not been initialized', FATAL)
+
+ if (.not. module_is_initialized) &
+ call error_mesg('astronomy_mod', 'module has not been initialized', FATAL)
!--------------------------------------------------------------------
! fill the output arguments with the eccentricity, obliquity and
! perihelion angle.
!--------------------------------------------------------------------
- ecc_out = ecc
- obliq_out = obliq
- per_out = per
-
-
-end subroutine get_orbital_parameters
-
-
-!> @brief set_ref_date_of_ae provides a means of specifying the reference
-!! date of the NH autumnal equinox for a particular year.
-!!
-!! @details set_ref_date_of_ae provides a means of specifying the reference
-!! date of the NH autumnal equinox for a particular year. It is only
-!! used if calls are made to the calandar versions of the routines
-!! diurnal_solar and daily_mean_solar. If the NOLEAP calendar is
-!! used, then the date of autumnal equinox will be the same every
-!! year. If JULIAN is used, then the date of autumnal equinox will
-!! return to the same value every 4th year.
-!!
-!! @param [in] Day of reference autumnal equinox
-!! @param [in] Month of reference autumnal equinox
-!! @param [in] Year of reference autumnal equinox
-!! @param [out] OPTIONAL: Second of reference autumnal equinox
-!! @param [out] OPTIONAL: Minute of reference autumnal equinox
-!! @param [out] OPTIONAL: Hour of reference autumnal equinox
-!!
-!! @throw FATAL, "astronomy_mod module has not been initialized"
-subroutine set_ref_date_of_ae (day_in,month_in,year_in, &
- second_in,minute_in,hour_in)
-
-integer, intent(in) :: day_in, month_in, year_in
-integer, intent(in), optional :: second_in, minute_in, hour_in
-
-!---------------------------------------------------------------------
-! exit if module has not been initialized.
-!---------------------------------------------------------------------
- if (.not. module_is_initialized) &
- call error_mesg ( 'astronomy_mod', &
- ' module has not been initialized', FATAL)
-
-!--------------------------------------------------------------------
-! save the input time of ae specification into a time_type module
-! variable autumnal_eq_ref.
-!--------------------------------------------------------------------
- day_ae = day_in
- month_ae = month_in
- year_ae = year_in
-
- if (present(second_in)) then
- second_ae = second_in
- minute_ae = minute_in
- hour_ae = hour_in
- else
- second_ae = 0
- minute_ae = 0
- hour_ae = 0
- endif
-
- autumnal_eq_ref = set_date (year_ae,month_ae,day_ae, &
- hour_ae,minute_ae,second_ae)
-!---------------------------------------------------------------------
+ ecc_out = real(ecc, FMS_AST_KIND_)
+ obliq_out = real(obliq, FMS_AST_KIND_)
+ per_out = real(per, FMS_AST_KIND_)
-
-end subroutine set_ref_date_of_ae
-
-
-!> @brief get_ref_date_of_ae retrieves the reference date of the autumnal
-!! equinox as integer variables.
-!!
-!! @throw FATAL, "astronomy_mod module has not been initialized"
-!!
-!! @param [out] Day of reference autumnal equinox
-!! @param [out] Month of reference autumnal equinox
-!! @param [out] Year of reference autumnal equinox
-!! @param [out] Second of reference autumnal equinox
-!! @param [out] Minute of reference autumnal equinox
-!! @param [out] Hour of reference autumnal equinox
-subroutine get_ref_date_of_ae (day_out,month_out,year_out,&
- second_out,minute_out,hour_out)
-
-integer, intent(out) :: day_out, month_out, year_out, &
- second_out, minute_out, hour_out
-
-!---------------------------------------------------------------------
-! exit if module has not been initialized.
-!---------------------------------------------------------------------
- if (.not. module_is_initialized) &
- call error_mesg ( 'astronomy_mod', &
- ' module has not been initialized', FATAL)
-
-!---------------------------------------------------------------------
-! fill the output fields with the proper module data.
-!---------------------------------------------------------------------
- day_out = day_ae
- month_out = month_ae
- year_out = year_ae
- second_out = second_ae
- minute_out = minute_ae
- hour_out = hour_ae
-
-
-end subroutine get_ref_date_of_ae
+end subroutine GET_ORBITAL_PARAMETERS_
!> @brief diurnal_solar_2d returns 2d fields of cosine of zenith angle,
@@ -811,26 +138,28 @@ end subroutine get_ref_date_of_ae
!!
!! @throw FATAL, "astronomy_mod time_since_ae not between 0 and 2pi"
!! @throw FATAL, "astronomy_mod gmt not between 0 and 2pi"
-subroutine diurnal_solar_2d (lat, lon, gmt, time_since_ae, cosz, &
- fracday, rrsun, dt, allow_negative_cosz, &
- half_day_out)
-real, dimension(:,:), intent(in) :: lat, lon
-real, intent(in) :: gmt, time_since_ae
-real, dimension(:,:), intent(out) :: cosz, fracday
-real, intent(out) :: rrsun
-real, intent(in), optional :: dt
-logical, intent(in), optional :: allow_negative_cosz
-real, dimension(:,:), intent(out), optional :: half_day_out
+subroutine DIURNAL_SOLAR_2D_(lat, lon, gmt, time_since_ae, cosz, &
+ fracday, rrsun, dt, allow_negative_cosz, &
+ half_day_out)
+real(kind=FMS_AST_KIND_), dimension(:,:), intent(in) :: lat, lon
+real(kind=FMS_AST_KIND_), intent(in) :: gmt, time_since_ae
+real(kind=FMS_AST_KIND_), dimension(:,:), intent(out) :: cosz, fracday
+real(kind=FMS_AST_KIND_), intent(out) :: rrsun
+real(kind=FMS_AST_KIND_), intent(in), optional :: dt
+logical, intent(in), optional :: allow_negative_cosz
+real(kind=FMS_AST_KIND_), dimension(:,:), intent(out), optional :: half_day_out
!---------------------------------------------------------------------
! local variables
- real, dimension(size(lat,1),size(lat,2)) :: t, tt, h, aa, bb, &
- st, stt, sh
- real :: ang, dec
- logical :: Lallow_negative
+
+real(kind=FMS_AST_KIND_), dimension(size(lat,1),size(lat,2)) :: t, tt, h, aa, bb, st, stt, sh
+real(kind=FMS_AST_KIND_) :: ang, dec
+logical :: Lallow_negative
+integer, parameter :: lkind = FMS_AST_KIND_
+
!---------------------------------------------------------------------
! local variables
@@ -854,72 +183,79 @@ real, dimension(:,:), intent(out), optional :: half_day_out
!--------------------------------------------------------------------
! be sure the time in the annual cycle is legitimate.
!---------------------------------------------------------------------
- if (time_since_ae < 0.0 .or. time_since_ae > twopi) &
- call error_mesg('astronomy_mod', &
- 'time_since_ae not between 0 and 2pi', FATAL)
+
+ if (time_since_ae < 0.0_lkind .or. time_since_ae > real(twopi, FMS_AST_KIND_)) &
+ call error_mesg('astronomy_mod', 'time_since_ae not between 0 and 2pi', FATAL)
!--------------------------------------------------------------------
! be sure the time at longitude = 0.0 is legitimate.
!---------------------------------------------------------------------
- if (gmt < 0.0 .or. gmt > twopi) &
- call error_mesg('astronomy_mod', &
- 'gmt not between 0 and 2pi', FATAL)
+
+ if (gmt < 0.0_lkind .or. gmt > real(twopi, FMS_AST_KIND_)) &
+ call error_mesg('astronomy_mod', 'gmt not between 0 and 2pi', FATAL)
!---------------------------------------------------------------------
!> define the orbital angle (location in year), solar declination and
!! earth sun distance factor. use functions contained in this module.
!---------------------------------------------------------------------
- ang = angle(time_since_ae)
- dec = declination(ang)
- rrsun = r_inv_squared(ang)
+
+ ang = angle(time_since_ae)
+ dec = declination(ang)
+ rrsun = r_inv_squared(ang)
!---------------------------------------------------------------------
!> define terms needed in the cosine zenith angle equation.
!--------------------------------------------------------------------
- aa = sin(lat)*sin(dec)
- bb = cos(lat)*cos(dec)
+
+ aa = sin(lat)*sin(dec)
+ bb = cos(lat)*cos(dec)
!---------------------------------------------------------------------
!> define local time. force it to be between -pi and pi.
!--------------------------------------------------------------------
- t = gmt + lon - PI
- where(t >= PI) t = t - twopi
- where(t < -PI) t = t + twopi
- Lallow_negative = .false.
- if (present(allow_negative_cosz)) then
- if (allow_negative_cosz) Lallow_negative = .true.
- end if
+ t = gmt + lon - real(PI, FMS_AST_KIND_)
+ where(t >= real(PI, FMS_AST_KIND_)) t = t - real(twopi, FMS_AST_KIND_)
+ where(t < real(-PI,FMS_AST_KIND_)) t = t + real(twopi, FMS_AST_KIND_)
+
+ Lallow_negative = .false.
+ if (present(allow_negative_cosz)) then
+ if (allow_negative_cosz) Lallow_negative = .true.
+ end if
+
!---------------------------------------------------------------------
!> perform a time integration to obtain cosz and fracday if desired.
!! output is valid over the period from t to t + dt.
!--------------------------------------------------------------------
- h = half_day (lat,dec)
- if ( present(half_day_out) ) then
- half_day_out = h
- end if
+ h = half_day(lat,dec)
+
+ if ( present(half_day_out) ) then
+ half_day_out = h
+ end if
- if ( present(dt) ) then ! (perform time averaging)
- tt = t + dt
- st = sin(t)
- stt = sin(tt)
- sh = sin(h)
- cosz = 0.0
+ if ( present(dt) ) then ! (perform time averaging)
+ tt = t + dt
+ st = sin(t)
+ stt = sin(tt)
+ sh = sin(h)
+ cosz = 0.0_lkind
if (.not. Lallow_negative) then
!-------------------------------------------------------------------
! case 1: entire averaging period is before sunrise.
!-------------------------------------------------------------------
- where (t < -h .and. tt < -h) cosz = 0.0
+
+ where (t < -h .and. tt < -h) cosz = 0.0_lkind
!-------------------------------------------------------------------
! case 2: averaging period begins before sunrise, ends after sunrise
! but before sunset
!-------------------------------------------------------------------
- where ( (tt+h) /= 0.0 .and. t < -h .and. abs(tt) <= h) &
- cosz = aa + bb*(stt + sh)/ (tt + h)
+
+ where ( (tt+h) /= 0.0_lkind .and. t < -h .and. abs(tt) <= h) &
+ cosz = aa + bb*(stt + sh)/ (tt + h)
!-------------------------------------------------------------------
! case 3: averaging period begins before sunrise, ends after sunset,
@@ -927,36 +263,43 @@ real, dimension(:,:), intent(out), optional :: half_day_out
! past the next day's sunrise, but if averaging period is less than
! a half- day (pi) that circumstance will never occur.
!-------------------------------------------------------------------
- where (t < -h .and. h /= 0.0 .and. h < tt) &
- cosz = aa + bb*( sh + sh)/(h+h)
+
+ where (t < -h .and. h /= 0.0_lkind .and. h < tt) &
+ cosz = aa + bb*( sh + sh)/(h+h)
!-------------------------------------------------------------------
! case 4: averaging period begins after sunrise, ends before sunset.
!-------------------------------------------------------------------
where ( abs(t) <= h .and. abs(tt) <= h) &
- cosz = aa + bb*(stt - st)/ (tt - t)
+
+ cosz = aa + bb*(stt - st)/ (tt - t)
!-------------------------------------------------------------------
! case 5: averaging period begins after sunrise, ends after sunset.
! modify when averaging period extends past the next day's sunrise.
!-------------------------------------------------------------------
- where ((h-t) /= 0.0 .and. abs(t) <= h .and. h < tt) &
- cosz = aa + bb*(sh - st)/(h-t)
+
+ where ((h-t) /= 0.0_lkind .and. abs(t) <= h .and. h < tt) &
+ cosz = aa + bb*(sh - st)/(h-t)
!-------------------------------------------------------------------
! case 6: averaging period begins after sunrise , ends after the
! next day's sunrise. note that this includes the case when the
! day length is one day (h = pi).
!-------------------------------------------------------------------
- where (twopi - h < tt .and. (tt+h-twopi) /= 0.0 .and. t <= h ) &
- cosz = (cosz*(h - t) + (aa*(tt + h - twopi) + &
- bb*(stt + sh))) / ((h - t) + (tt + h - twopi))
+
+ where (real(twopi, FMS_AST_KIND_) - h < tt .and. &
+ (tt+h-real(twopi, FMS_AST_KIND_)) /= 0.0_lkind .and. t <= h ) &
+ cosz = (cosz*(h - t) + (aa*(tt + h - real(twopi, FMS_AST_KIND_)) + bb &
+ * (stt + sh))) / ((h - t) + (tt + h - real(twopi, FMS_AST_KIND_)))
!-------------------------------------------------------------------
! case 7: averaging period begins after sunset and ends before the
! next day's sunrise
!-------------------------------------------------------------------
- where( h < t .and. twopi - h >= tt ) cosz = 0.0
+
+ where(h < t .and. real(twopi, FMS_AST_KIND_) - h >= tt) &
+ cosz = 0.0_lkind
!-------------------------------------------------------------------
! case 8: averaging period begins after sunset and ends after the
@@ -964,12 +307,13 @@ real, dimension(:,:), intent(out), optional :: half_day_out
! averaging period is less than a half-day (pi) the latter
! circumstance will never occur.
!-----------------------------------------------------------------
- where( h < t .and. twopi - h < tt )
- cosz = aa + bb*(stt + sh) / (tt + h - twopi)
+
+ where (h < t .and. real(twopi, FMS_AST_KIND_) - h < tt)
+ cosz = aa + bb*(stt + sh) / (tt + h - real(twopi, FMS_AST_KIND_))
end where
else
- cosz = aa + bb*(stt - st)/ (tt - t)
+ cosz = aa + bb*(stt - st)/ (tt - t)
end if
@@ -978,45 +322,47 @@ real, dimension(:,:), intent(out), optional :: half_day_out
! day fraction is the fraction of the averaging period contained
! within the (-h,h) period.
!-------------------------------------------------------------------
- where (t < -h .and. tt < -h) fracday = 0.0
- where (t < -h .and. abs(tt) <= h) fracday = (tt + h )/dt
- where (t < -h .and. h < tt) fracday = ( h + h )/dt
- where (abs(t) <= h .and. abs(tt) <= h) fracday = (tt - t )/dt
- where (abs(t) <= h .and. h < tt) fracday = ( h - t )/dt
- where ( h < t ) fracday = 0.0
- where (twopi - h < tt) fracday = fracday + &
- (tt + h - &
- twopi)/dt
+
+ where (t < -h .and. tt < -h) fracday = 0.0_lkind
+ where (t < -h .and. abs(tt) <= h) fracday = (tt + h )/dt
+ where (t < -h .and. h < tt) fracday = ( h + h )/dt
+ where (abs(t) <= h .and. abs(tt) <= h) fracday = (tt - t )/dt
+ where (abs(t) <= h .and. h < tt) fracday = ( h - t )/dt
+ where (h < t) fracday = 0.0_lkind
+ where (real(twopi, FMS_AST_KIND_) - h < tt) &
+ fracday = fracday + (tt + h - real(twopi, FMS_AST_KIND_))/dt
!----------------------------------------------------------------------
!> if instantaneous values are desired, define cosz at time t.
!----------------------------------------------------------------------
- else ! (no time averaging)
+ else ! (no time averaging)
if (.not. Lallow_negative) then
- where (abs(t) < h)
- cosz = aa + bb*cos(t)
- fracday = 1.0
- elsewhere
- cosz = 0.0
- fracday = 0.0
- end where
+ where (abs(t) < h)
+ cosz = aa + bb*cos(t)
+ fracday = 1.0_lkind
+ elsewhere
+ cosz = 0.0_lkind
+ fracday = 0.0_lkind
+ end where
else
- cosz = aa + bb*cos(t)
- where (abs(t) < h)
- fracday = 1.0
- elsewhere
- fracday = 0.0
- end where
+ cosz = aa + bb*cos(t)
+ where (abs(t) < h)
+ fracday = 1.0_lkind
+ elsewhere
+ fracday = 0.0_lkind
+ end where
end if
- end if
+ end if
+
!----------------------------------------------------------------------
!> Check that cosz is not negative, if desired.
!----------------------------------------------------------------------
- if (.not. Lallow_negative) then
- cosz = max(0.0, cosz)
- end if
-end subroutine diurnal_solar_2d
+ if (.not. Lallow_negative) then
+ cosz = max(0.0_lkind, cosz)
+ end if
+
+end subroutine DIURNAL_SOLAR_2D_
!> @brief diurnal_solar_1d takes 1-d input fields, adds a second dimension
@@ -1034,39 +380,43 @@ end subroutine diurnal_solar_2d
!! averaged. this produces averaged output rather than instantaneous.
!! @param [in] Allow negative values for cosz?
!! @param [out] half_day_out
-subroutine diurnal_solar_1d (lat, lon, gmt, time_since_ae, cosz, &
+
+subroutine DIURNAL_SOLAR_1D_(lat, lon, gmt, time_since_ae, cosz, &
fracday, rrsun, dt, allow_negative_cosz, &
half_day_out)
!---------------------------------------------------------------------
-real, dimension(:), intent(in) :: lat, lon
-real, intent(in) :: gmt, time_since_ae
-real, dimension(:), intent(out) :: cosz, fracday
-real, intent(out) :: rrsun
-real, intent(in), optional :: dt
-logical, intent(in), optional :: allow_negative_cosz
-real, dimension(:), intent(out), optional :: half_day_out
+
+real(kind=FMS_AST_KIND_), dimension(:), intent(in) :: lat, lon
+real(kind=FMS_AST_KIND_), intent(in) :: gmt, time_since_ae
+real(kind=FMS_AST_KIND_), dimension(:), intent(out) :: cosz, fracday
+real(kind=FMS_AST_KIND_), intent(out) :: rrsun
+real(kind=FMS_AST_KIND_), intent(in), optional :: dt
+logical, intent(in), optional :: allow_negative_cosz
+real(kind=FMS_AST_KIND_), dimension(:), intent(out), optional :: half_day_out
!---------------------------------------------------------------------
! local variables
!---------------------------------------------------------------------
- real, dimension(size(lat),1) :: lat_2d, lon_2d, cosz_2d, &
- fracday_2d,halfday_2d
+
+real(kind=FMS_AST_KIND_), dimension(size(lat),1) :: lat_2d, lon_2d, cosz_2d, fracday_2d, halfday_2d
!--------------------------------------------------------------------
!> define 2-d versions of input data arrays.
!--------------------------------------------------------------------
- lat_2d(:,1) = lat
- lon_2d(:,1) = lon
+
+ lat_2d(:,1) = lat
+ lon_2d(:,1) = lon
+
!--------------------------------------------------------------------
!> call diurnal_solar_2d to calculate astronomy fields.
!--------------------------------------------------------------------
! if (present(dt)) then
- call diurnal_solar_2d (lat_2d, lon_2d, gmt, time_since_ae,&
- cosz_2d, fracday_2d, rrsun, dt=dt, &
- allow_negative_cosz=allow_negative_cosz, &
- half_day_out=halfday_2d)
+
+ call diurnal_solar(lat_2d, lon_2d, gmt, time_since_ae, &
+ cosz_2d, fracday_2d, rrsun, dt=dt, &
+ allow_negative_cosz=allow_negative_cosz, half_day_out=halfday_2d)
! else
! call diurnal_solar_2d (lat_2d, lon_2d, gmt, time_since_ae, &
! cosz_2d, fracday_2d, rrsun)
@@ -1076,13 +426,15 @@ real, dimension(:), intent(out), optional :: half_day_out
!> place output fields into 1-d arguments for return to
!! calling routine.
!-------------------------------------------------------------------
- fracday = fracday_2d(:,1)
- cosz = cosz_2d (:,1)
- if (present(half_day_out)) then
- half_day_out = halfday_2d(:,1)
- end if
-end subroutine diurnal_solar_1d
+ fracday = fracday_2d(:,1)
+ cosz = cosz_2d (:,1)
+ if (present(half_day_out)) then
+ half_day_out = halfday_2d(:,1)
+ end if
+
+end subroutine DIURNAL_SOLAR_1D_
+
!> @brief diurnal_solar_0d takes scalar input fields, makes them into 2d
@@ -1100,35 +452,38 @@ end subroutine diurnal_solar_1d
!! averaged. this produces averaged output rather than instantaneous.
!! @param [in] Allow negative values for cosz?
!! @param [out] half_day_out
-subroutine diurnal_solar_0d (lat, lon, gmt, time_since_ae, cosz, &
+
+subroutine DIURNAL_SOLAR_0D_(lat, lon, gmt, time_since_ae, cosz, &
fracday, rrsun, dt, allow_negative_cosz, &
half_day_out)
-real, intent(in) :: lat, lon, gmt, time_since_ae
-real, intent(out) :: cosz, fracday, rrsun
-real, intent(in), optional :: dt
-logical,intent(in),optional :: allow_negative_cosz
-real, intent(out), optional :: half_day_out
+real(kind=FMS_AST_KIND_), intent(in) :: lat, lon, gmt, time_since_ae
+real(kind=FMS_AST_KIND_), intent(out) :: cosz, fracday, rrsun
+real(kind=FMS_AST_KIND_), intent(in), optional :: dt
+logical, intent(in), optional :: allow_negative_cosz
+real(kind=FMS_AST_KIND_), intent(out), optional :: half_day_out
!--------------------------------------------------------------------
! local variables:
!--------------------------------------------------------------------
- real, dimension(1,1) :: lat_2d, lon_2d, cosz_2d, fracday_2d, halfday_2d
+
+real(kind=FMS_AST_KIND_), dimension(1,1) :: lat_2d, lon_2d, cosz_2d, fracday_2d, halfday_2d
!---------------------------------------------------------------------
!> create 2d arrays from the scalar input fields.
!---------------------------------------------------------------------
- lat_2d = lat
- lon_2d = lon
+
+ lat_2d = lat
+ lon_2d = lon
!--------------------------------------------------------------------
!> call diurnal_solar_2d to calculate astronomy fields.
!--------------------------------------------------------------------
! if (present(dt)) then
- call diurnal_solar_2d (lat_2d, lon_2d, gmt, time_since_ae, &
- cosz_2d, fracday_2d, rrsun, dt=dt, &
- allow_negative_cosz=allow_negative_cosz, &
- half_day_out=halfday_2d)
+
+ call diurnal_solar(lat_2d, lon_2d, gmt, time_since_ae, &
+ cosz_2d, fracday_2d, rrsun, dt=dt, &
+ allow_negative_cosz=allow_negative_cosz, half_day_out=halfday_2d)
! else
! call diurnal_solar_2d (lat_2d, lon_2d, gmt, time_since_ae, &
! cosz_2d, fracday_2d, rrsun)
@@ -1137,13 +492,14 @@ real, intent(out), optional :: half_day_out
!-------------------------------------------------------------------
!> place output fields into scalars for return to calling routine.
!-------------------------------------------------------------------
- fracday = fracday_2d(1,1)
- cosz = cosz_2d(1,1)
- if (present(half_day_out)) then
- half_day_out = halfday_2d(1,1)
- end if
-end subroutine diurnal_solar_0d
+ fracday = fracday_2d(1,1)
+ cosz = cosz_2d(1,1)
+ if (present(half_day_out)) then
+ half_day_out = halfday_2d(1,1)
+ end if
+
+end subroutine DIURNAL_SOLAR_0D_
!> @brief diurnal_solar_cal_2d receives time_type inputs, converts
@@ -1164,71 +520,75 @@ end subroutine diurnal_solar_0d
!!
!! @throw FATAL, "astronomy_mod radiation time step must be no longer than 12 hrs"
!! @throw FATAL, "astronomy_mod radiation time step must not be an integral number of days"
-subroutine diurnal_solar_cal_2d (lat, lon, time, cosz, fracday, &
+
+subroutine DIURNAL_SOLAR_CAL_2D_(lat, lon, time, cosz, fracday, &
rrsun, dt_time, allow_negative_cosz, &
half_day_out)
!-------------------------------------------------------------------
-real, dimension(:,:), intent(in) :: lat, lon
-type(time_type), intent(in) :: time
-real, dimension(:,:), intent(out) :: cosz, fracday
-real, intent(out) :: rrsun
-type(time_type), intent(in), optional :: dt_time
-logical, intent(in), optional :: allow_negative_cosz
-real, dimension(:,:), intent(out), optional :: half_day_out
+
+real(kind=FMS_AST_KIND_), dimension(:,:), intent(in) :: lat, lon
+type(time_type), intent(in) :: time
+real(kind=FMS_AST_KIND_), dimension(:,:), intent(out) :: cosz, fracday
+real(kind=FMS_AST_KIND_), intent(out) :: rrsun
+type(time_type), intent(in), optional :: dt_time
+logical, intent(in), optional :: allow_negative_cosz
+real(kind=FMS_AST_KIND_), dimension(:,:), intent(out), optional :: half_day_out
+
!---------------------------------------------------------------------
!---------------------------------------------------------------------
! local variables
!---------------------------------------------------------------------
- real :: dt
- real :: gmt, time_since_ae
+
+real(kind=FMS_AST_KIND_) :: dt
+real(kind=FMS_AST_KIND_) :: gmt, time_since_ae
+integer, parameter :: lkind = FMS_AST_KIND_
!---------------------------------------------------------------------
!> Extract time of day (gmt) from time_type variable time with
!! function universal_time.
!---------------------------------------------------------------------
- gmt = universal_time(time)
+
+ gmt = real(universal_time(time), FMS_AST_KIND_)
!---------------------------------------------------------------------
!> Extract the time of year (time_since_ae) from time_type variable
!! time using the function orbital_time.
!---------------------------------------------------------------------
- time_since_ae = orbital_time(time)
+
+ time_since_ae = real(orbital_time(time), FMS_AST_KIND_)
!---------------------------------------------------------------------
!> Convert optional time_type variable dt_time (length of averaging
!! period) to a real variable dt with the function universal_time.
!---------------------------------------------------------------------
- if (present(dt_time)) then
- dt = universal_time(dt_time)
- if (dt > PI) then
- call error_mesg ( 'astronomy_mod', &
- 'radiation time step must be no longer than 12 hrs', &
- FATAL)
- endif
- if (dt == 0.0) then
- call error_mesg ( 'astronomy_mod', &
- 'radiation time step must not be an integral &
- &number of days', FATAL)
+
+ if (present(dt_time)) then
+ dt = real(universal_time(dt_time), FMS_AST_KIND_)
+ if (dt > real(PI, FMS_AST_KIND_)) then
+ call error_mesg('astronomy_mod', 'radiation time step must be no longer than 12 hrs', FATAL)
endif
+ if (dt == 0.0_lkind) then
+ call error_mesg('astronomy_mod', 'radiation time step must not be an integral number of days', FATAL)
+ endif
!--------------------------------------------------------------------
!> Call diurnal_solar_2d to calculate astronomy fields, with or
!! without the optional argument dt.
!--------------------------------------------------------------------
- call diurnal_solar_2d (lat, lon, gmt, time_since_ae, cosz, &
- fracday, rrsun, dt=dt, &
- allow_negative_cosz=allow_negative_cosz, &
- half_day_out=half_day_out)
- else
- call diurnal_solar_2d (lat, lon, gmt, time_since_ae, cosz, &
- fracday, rrsun, &
- allow_negative_cosz=allow_negative_cosz, &
- half_day_out=half_day_out)
- end if
-end subroutine diurnal_solar_cal_2d
+ call diurnal_solar(lat, lon, gmt, time_since_ae, cosz, &
+ fracday, rrsun, dt=dt, &
+ allow_negative_cosz=allow_negative_cosz, &
+ half_day_out=half_day_out)
+ else
+ call diurnal_solar(lat, lon, gmt, time_since_ae, cosz, &
+ fracday, rrsun, allow_negative_cosz=allow_negative_cosz, &
+ half_day_out=half_day_out)
+ end if
+
+end subroutine DIURNAL_SOLAR_CAL_2D_
!> @brief diurnal_solar_cal_1d receives time_type inputs, converts
@@ -1246,59 +606,64 @@ end subroutine diurnal_solar_cal_2d
!! averaged. this produces averaged output rather than instantaneous.
!! @param [in] Allow negative values for cosz?
!! @param [out] half_day_out
-subroutine diurnal_solar_cal_1d (lat, lon, time, cosz, fracday, &
+
+subroutine DIURNAL_SOLAR_CAL_1D_(lat, lon, time, cosz, fracday, &
rrsun, dt_time, allow_negative_cosz, &
half_day_out)
!--------------------------------------------------------------------
-real, dimension(:), intent(in) :: lat, lon
-type(time_type), intent(in) :: time
-real, dimension(:), intent(out) :: cosz, fracday
-real, intent(out) :: rrsun
-type(time_type), intent(in), optional :: dt_time
-logical, intent(in), optional :: allow_negative_cosz
-real, dimension(:), intent(out), optional :: half_day_out
+
+real(kind=FMS_AST_KIND_), dimension(:), intent(in) :: lat, lon
+type(time_type), intent(in) :: time
+real(kind=FMS_AST_KIND_), dimension(:), intent(out) :: cosz, fracday
+real(kind=FMS_AST_KIND_), intent(out) :: rrsun
+type(time_type), intent(in), optional :: dt_time
+logical, intent(in), optional :: allow_negative_cosz
+real(kind=FMS_AST_KIND_), dimension(:), intent(out), optional :: half_day_out
!--------------------------------------------------------------------
!-------------------------------------------------------------------
! local variables
!-------------------------------------------------------------------
- real, dimension(size(lat),1) :: lat_2d, lon_2d, cosz_2d, &
+
+real(kind=FMS_AST_KIND_), dimension(size(lat),1) :: lat_2d, lon_2d, cosz_2d, &
fracday_2d, halfday_2d
!--------------------------------------------------------------------
!> Define 2-d versions of input data arrays.
!--------------------------------------------------------------------
- lat_2d(:,1) = lat
- lon_2d(:,1) = lon
+
+ lat_2d(:,1) = lat
+ lon_2d(:,1) = lon
!--------------------------------------------------------------------
!> Call diurnal_solar_cal_2d to convert the time_types to reals and
!! then calculate the astronomy fields.
!--------------------------------------------------------------------
- if (present(dt_time)) then
- call diurnal_solar_cal_2d (lat_2d, lon_2d, time, cosz_2d, &
- fracday_2d, rrsun, dt_time=dt_time, &
- allow_negative_cosz=allow_negative_cosz, &
- half_day_out=halfday_2d)
- else
- call diurnal_solar_cal_2d (lat_2d, lon_2d, time, cosz_2d, &
- fracday_2d, rrsun, &
- allow_negative_cosz=allow_negative_cosz, &
- half_day_out=halfday_2d)
- end if
+
+ if (present(dt_time)) then
+ call diurnal_solar(lat_2d, lon_2d, time, cosz_2d, &
+ fracday_2d, rrsun, dt_time=dt_time, &
+ allow_negative_cosz=allow_negative_cosz, &
+ half_day_out=halfday_2d)
+ else
+ call diurnal_solar(lat_2d, lon_2d, time, cosz_2d, &
+ fracday_2d, rrsun, allow_negative_cosz=allow_negative_cosz, &
+ half_day_out=halfday_2d)
+ end if
!-------------------------------------------------------------------
!> Place output fields into 1-d arguments for return to
!! calling routine.
!-------------------------------------------------------------------
- fracday = fracday_2d(:,1)
- cosz = cosz_2d (:,1)
- if (present(half_day_out)) then
- half_day_out = halfday_2d(:,1)
- end if
-end subroutine diurnal_solar_cal_1d
+ fracday = fracday_2d(:,1)
+ cosz = cosz_2d (:,1)
+ if (present(half_day_out)) then
+ half_day_out = halfday_2d(:,1)
+ end if
+
+end subroutine DIURNAL_SOLAR_CAL_1D_
!> @brief diurnal_solar_cal_0d receives time_type inputs, converts them to real variables
@@ -1314,57 +679,63 @@ end subroutine diurnal_solar_cal_1d
!! to be averaged. this produces averaged output rather than instantaneous.
!! @param [in] allow_negative_cosz
!! @param [out] half_day_out
-subroutine diurnal_solar_cal_0d (lat, lon, time, cosz, fracday, &
+
+subroutine DIURNAL_SOLAR_CAL_0D_(lat, lon, time, cosz, fracday, &
rrsun, dt_time, allow_negative_cosz, &
half_day_out)
!---------------------------------------------------------------------
-real, intent(in) :: lat, lon
-type(time_type), intent(in) :: time
-real, intent(out) :: cosz, fracday, rrsun
-type(time_type), intent(in), optional :: dt_time
-logical, intent(in), optional :: allow_negative_cosz
-real, intent(out), optional :: half_day_out
+
+real(kind=FMS_AST_KIND_), intent(in) :: lat, lon
+type(time_type), intent(in) :: time
+real(kind=FMS_AST_KIND_), intent(out) :: cosz, fracday, rrsun
+type(time_type), intent(in), optional :: dt_time
+logical, intent(in), optional :: allow_negative_cosz
+real(kind=FMS_AST_KIND_), intent(out), optional :: half_day_out
+
!---------------------------------------------------------------------
!---------------------------------------------------------------------
! local variables
!---------------------------------------------------------------------
- real, dimension(1,1) :: lat_2d, lon_2d, cosz_2d, fracday_2d, halfday_2d
+
+real(kind=FMS_AST_KIND_), dimension(1,1) :: lat_2d, lon_2d, cosz_2d, fracday_2d, halfday_2d
!--------------------------------------------------------------------
!> Define 2-d versions of input data arrays.
!--------------------------------------------------------------------
- lat_2d = lat
- lon_2d = lon
+
+ lat_2d = lat
+ lon_2d = lon
!--------------------------------------------------------------------
!> Call diurnal_solar_cal_2d to convert the time_types to reals and
!! then calculate the astronomy fields.
!--------------------------------------------------------------------
- if (present(dt_time)) then
- call diurnal_solar_cal_2d (lat_2d, lon_2d, time, cosz_2d, &
- fracday_2d, rrsun, dt_time=dt_time, &
- allow_negative_cosz=allow_negative_cosz, &
- half_day_out=halfday_2d)
- else
- call diurnal_solar_cal_2d (lat_2d, lon_2d, time, cosz_2d, &
- fracday_2d, rrsun, &
- allow_negative_cosz=allow_negative_cosz, &
- half_day_out=halfday_2d)
- end if
+
+ if (present(dt_time)) then
+ call diurnal_solar(lat_2d, lon_2d, time, cosz_2d, &
+ fracday_2d, rrsun, dt_time=dt_time, &
+ allow_negative_cosz=allow_negative_cosz, &
+ half_day_out=halfday_2d)
+ else
+ call diurnal_solar(lat_2d, lon_2d, time, cosz_2d, &
+ fracday_2d, rrsun, allow_negative_cosz=allow_negative_cosz, &
+ half_day_out=halfday_2d)
+ end if
!-------------------------------------------------------------------
!> Place output fields into 1-d arguments for return to
!! calling routine.
!-------------------------------------------------------------------
- fracday= fracday_2d(1,1)
- cosz = cosz_2d(1,1)
- if (present(half_day_out)) then
- half_day_out = halfday_2d(1,1)
+
+ fracday = fracday_2d(1,1)
+ cosz = cosz_2d(1,1)
+ if (present(half_day_out)) then
+ half_day_out = halfday_2d(1,1)
end if
-end subroutine diurnal_solar_cal_0d
+end subroutine DIURNAL_SOLAR_CAL_0D_
!> @brief daily_mean_solar_2d computes the daily mean astronomical parameters for
@@ -1378,52 +749,57 @@ end subroutine diurnal_solar_cal_0d
!! to the mean distance at angle ang in the earth's orbit.
!!
!! @throw FATAL, "astronomy_mod time_since_ae not between 0 and 2pi"
-subroutine daily_mean_solar_2d (lat, time_since_ae, cosz, h_out, rr_out)
+
+subroutine DAILY_MEAN_SOLAR_2D_(lat, time_since_ae, cosz, h_out, rr_out)
!----------------------------------------------------------------------
-real, dimension(:,:), intent(in) :: lat
-real, intent(in) :: time_since_ae
-real, dimension(:,:), intent(out) :: cosz, h_out
-real, intent(out) :: rr_out
+real(kind=FMS_AST_KIND_), dimension(:,:), intent(in) :: lat
+real(kind=FMS_AST_KIND_), intent(in) :: time_since_ae
+real(kind=FMS_AST_KIND_), dimension(:,:), intent(out) :: cosz, h_out
+real(kind=FMS_AST_KIND_), intent(out) :: rr_out
!----------------------------------------------------------------------
!--------------------------------------------------------------------
! local variables
!--------------------------------------------------------------------
- real, dimension(size(lat,1),size(lat,2)) :: h
- real :: ang, dec, rr
+
+real(kind=FMS_AST_KIND_), dimension(size(lat,1),size(lat,2)) :: h
+real(kind=FMS_AST_KIND_) :: ang, dec, rr
+integer, parameter :: lkind = FMS_AST_KIND_
!--------------------------------------------------------------------
! be sure the time in the annual cycle is legitimate.
!---------------------------------------------------------------------
- if (time_since_ae < 0.0 .or. time_since_ae > twopi) &
- call error_mesg('astronomy_mod', &
- 'time_since_ae not between 0 and 2pi', FATAL)
+
+ if (time_since_ae < 0.0_lkind .or. time_since_ae > real(twopi, FMS_AST_KIND_)) &
+ call error_mesg('astronomy_mod', 'time_since_ae not between 0 and 2pi', FATAL)
!---------------------------------------------------------------------
!> Define the orbital angle (location in year), solar declination,
!! half-day length and earth sun distance factor. Use functions
!! contained in this module.
!---------------------------------------------------------------------
- ang = angle (time_since_ae)
- dec = declination(ang)
- h = half_day (lat, dec)
- rr = r_inv_squared (ang)
+
+ ang = angle (time_since_ae)
+ dec = declination (ang)
+ h = half_day (lat, dec)
+ rr = r_inv_squared (ang)
!---------------------------------------------------------------------
!> Where the entire day is dark, define cosz to be zero. otherwise
!! use the standard formula. Define the daylight fraction and earth-
!! sun distance.
!---------------------------------------------------------------------
- where (h == 0.0)
- cosz = 0.0
- elsewhere
- cosz = sin(lat)*sin(dec) + cos(lat)*cos(dec)*sin(h)/h
- end where
- h_out = h/PI
- rr_out = rr
-end subroutine daily_mean_solar_2d
+ where (h == 0.0_lkind)
+ cosz = 0.0_lkind
+ elsewhere
+ cosz = sin(lat)*sin(dec) + cos(lat)*cos(dec)*sin(h)/h
+ end where
+ h_out = h/real(PI, FMS_AST_KIND_)
+ rr_out = rr
+
+end subroutine DAILY_MEAN_SOLAR_2D_
!> @brief daily_mean_solar_1d takes 1-d input fields, adds a second dimension
@@ -1436,40 +812,46 @@ end subroutine daily_mean_solar_2d
!! @param [out] 2-d array of half-day lengths at the latitudes
!! @param [out] the inverse of the square of the earth-sun distance relative
!! to the mean distance at angle ang in the earth's orbit.
-subroutine daily_mean_solar_1d (lat, time_since_ae, cosz, h_out, rr_out)
+
+subroutine DAILY_MEAN_SOLAR_1D_(lat, time_since_ae, cosz, h_out, rr_out)
!----------------------------------------------------------------------
-real, intent(in), dimension(:) :: lat
-real, intent(in) :: time_since_ae
-real, intent(out), dimension(size(lat(:))) :: cosz
-real, intent(out), dimension(size(lat(:))) :: h_out
-real, intent(out) :: rr_out
+real(kind=FMS_AST_KIND_), intent(in), dimension(:) :: lat
+real(kind=FMS_AST_KIND_), intent(in) :: time_since_ae
+real(kind=FMS_AST_KIND_), intent(out), dimension(size(lat(:))) :: cosz
+real(kind=FMS_AST_KIND_), intent(out), dimension(size(lat(:))) :: h_out
+real(kind=FMS_AST_KIND_), intent(out) :: rr_out
+
!----------------------------------------------------------------------
!----------------------------------------------------------------------
! local variables
!----------------------------------------------------------------------
- real, dimension(size(lat),1) :: lat_2d, cosz_2d, hout_2d
+
+real(kind=FMS_AST_KIND_), dimension(size(lat),1) :: lat_2d, cosz_2d, hout_2d
!--------------------------------------------------------------------
!> define 2-d versions of input data array.
!--------------------------------------------------------------------
- lat_2d(:,1) = lat
+
+ lat_2d(:,1) = lat
!--------------------------------------------------------------------
!> call daily_mean_solar_2d to calculate astronomy fields.
!--------------------------------------------------------------------
- call daily_mean_solar_2d (lat_2d, time_since_ae, cosz_2d, &
+
+ call daily_mean_solar(lat_2d, time_since_ae, cosz_2d, &
hout_2d, rr_out)
!-------------------------------------------------------------------
!> place output fields into 1-d arguments for return to
!! calling routine.
!-------------------------------------------------------------------
- h_out = hout_2d(:,1)
- cosz = cosz_2d(:,1)
-end subroutine daily_mean_solar_1d
+ h_out = hout_2d(:,1)
+ cosz = cosz_2d(:,1)
+
+end subroutine DAILY_MEAN_SOLAR_1D_
!> @brief daily_mean_solar_2level takes 1-d input fields, adds a second
@@ -1481,39 +863,45 @@ end subroutine daily_mean_solar_1d
!! @param [out] Cosine of solar zenith angle
!! @param [out] Shortwave flux factor: cosine of zenith angle * daylight fraction /
!! (earth-sun distance squared)
-subroutine daily_mean_solar_2level (lat, time_since_ae, cosz, solar)
+
+subroutine DAILY_MEAN_SOLAR_2LEVEL_(lat, time_since_ae, cosz, solar)
!----------------------------------------------------------------------
-real, intent(in), dimension(:) :: lat
-real, intent(in) :: time_since_ae
-real, intent(out), dimension(size(lat(:))) :: cosz, solar
+real(kind=FMS_AST_KIND_), intent(in), dimension(:) :: lat
+real(kind=FMS_AST_KIND_), intent(in) :: time_since_ae
+real(kind=FMS_AST_KIND_), intent(out), dimension(size(lat(:))) :: cosz, solar
+
!----------------------------------------------------------------------
!----------------------------------------------------------------------
! local variables
!----------------------------------------------------------------------
- real, dimension(size(lat),1) :: lat_2d, cosz_2d, hout_2d
- real :: rr_out
+
+real(kind=FMS_AST_KIND_), dimension(size(lat),1) :: lat_2d, cosz_2d, hout_2d
+real(kind=FMS_AST_KIND_) :: rr_out
!--------------------------------------------------------------------
!> define 2-d versions of input data array.
!--------------------------------------------------------------------
- lat_2d(:,1) = lat
+
+ lat_2d(:,1) = lat
!--------------------------------------------------------------------
!> call daily_mean_solar_2d to calculate astronomy fields.
!--------------------------------------------------------------------
- call daily_mean_solar_2d (lat_2d, time_since_ae, cosz_2d, &
+
+ call daily_mean_solar(lat_2d, time_since_ae, cosz_2d, &
hout_2d, rr_out)
!-------------------------------------------------------------------
!> place output fields into 1-d arguments for return to
!! calling routine.
!-------------------------------------------------------------------
- solar = cosz_2d(:,1)*hout_2d(:,1)*rr_out
- cosz = cosz_2d(:,1)
-end subroutine daily_mean_solar_2level
+ solar = cosz_2d(:,1)*hout_2d(:,1)*rr_out
+ cosz = cosz_2d(:,1)
+
+end subroutine DAILY_MEAN_SOLAR_2LEVEL_
!> @brief daily_mean_solar_1d takes 1-d input fields, adds a second dimension
@@ -1526,34 +914,39 @@ end subroutine daily_mean_solar_2level
!! @param [out] 2-d array of half-day lengths at the latitudes
!! @param [out] the inverse of the square of the earth-sun distance relative to
!! the mean distance at angle ang in the earth's orbit.
-subroutine daily_mean_solar_0d (lat, time_since_ae, cosz, h_out, rr_out)
-real, intent(in) :: lat, time_since_ae
-real, intent(out) :: cosz, h_out, rr_out
+subroutine DAILY_MEAN_SOLAR_0D_(lat, time_since_ae, cosz, h_out, rr_out)
+
+real(kind=FMS_AST_KIND_), intent(in) :: lat, time_since_ae
+real(kind=FMS_AST_KIND_), intent(out) :: cosz, h_out, rr_out
!--------------------------------------------------------------------
! local variables
!--------------------------------------------------------------------
- real, dimension(1,1) :: lat_2d, cosz_2d, hout_2d
+
+real(kind=FMS_AST_KIND_), dimension(1,1) :: lat_2d, cosz_2d, hout_2d
!--------------------------------------------------------------------
!> define 2-d versions of input data array.
!--------------------------------------------------------------------
- lat_2d = lat
+
+ lat_2d = lat
!--------------------------------------------------------------------
!> call daily_mean_solar_2d to calculate astronomy fields.
!--------------------------------------------------------------------
- call daily_mean_solar_2d (lat_2d, time_since_ae, cosz_2d, &
+
+ call daily_mean_solar(lat_2d, time_since_ae, cosz_2d, &
hout_2d, rr_out)
!-------------------------------------------------------------------
!> return output fields to scalars for return to calling routine.
!-------------------------------------------------------------------
- h_out = hout_2d(1,1)
- cosz = cosz_2d(1,1)
-end subroutine daily_mean_solar_0d
+ h_out = hout_2d(1,1)
+ cosz = cosz_2d(1,1)
+
+end subroutine DAILY_MEAN_SOLAR_0D_
!> @brief daily_mean_solar_cal_2d receives time_type inputs, converts
@@ -1567,35 +960,38 @@ end subroutine daily_mean_solar_0d
!! @param [out] Earth-Sun distance (r) relative to semi-major axis of orbital ellipse (a):(a/r)**2
!!
!! @throw FATAL, "astronomy_mod time_since_ae not between 0 and 2pi"
-subroutine daily_mean_solar_cal_2d (lat, time, cosz, fracday, rrsun)
+
+subroutine DAILY_MEAN_SOLAR_CAL_2D_(lat, time, cosz, fracday, rrsun)
!-------------------------------------------------------------------
-real, dimension(:,:), intent(in) :: lat
-type(time_type), intent(in) :: time
-real, dimension(:,:), intent(out) :: cosz, fracday
-real, intent(out) :: rrsun
+real(kind=FMS_AST_KIND_), dimension(:,:), intent(in) :: lat
+type(time_type), intent(in) :: time
+real(kind=FMS_AST_KIND_), dimension(:,:), intent(out) :: cosz, fracday
+real(kind=FMS_AST_KIND_), intent(out) :: rrsun
!-------------------------------------------------------------------
!-------------------------------------------------------------------
! local variables
!-------------------------------------------------------------------
- real :: time_since_ae
+
+real(kind=FMS_AST_KIND_) :: time_since_ae
+integer, parameter :: lkind = FMS_AST_KIND_
!--------------------------------------------------------------------
! be sure the time in the annual cycle is legitimate.
!---------------------------------------------------------------------
- time_since_ae = orbital_time(time)
- if (time_since_ae < 0.0 .or. time_since_ae > twopi) &
- call error_mesg ('astronomy_mod', &
- 'time_since_ae not between 0 and 2pi', FATAL)
+
+ time_since_ae = real(orbital_time(time), FMS_AST_KIND_)
+ if (time_since_ae < 0.0_lkind .or. time_since_ae > real(twopi, FMS_AST_KIND_)) &
+ call error_mesg('astronomy_mod', 'time_since_ae not between 0 and 2pi', FATAL)
!--------------------------------------------------------------------
! call daily_mean_solar_2d to calculate astronomy fields.
!--------------------------------------------------------------------
- call daily_mean_solar_2d (lat, time_since_ae, cosz, &
- fracday, rrsun)
-end subroutine daily_mean_solar_cal_2d
+ call daily_mean_solar(lat, time_since_ae, cosz, fracday, rrsun)
+
+end subroutine DAILY_MEAN_SOLAR_CAL_2D_
!> @brief daily_mean_solar_cal_1d receives time_type inputs, converts
@@ -1607,39 +1003,42 @@ end subroutine daily_mean_solar_cal_2d
!! @param [out] Cosine of solar zenith angle
!! @param [out] Daylight fraction of time interval
!! @param [out] Earth-Sun distance (r) relative to semi-major axis of orbital ellipse (a):(a/r)**2
-subroutine daily_mean_solar_cal_1d (lat, time, cosz, fracday, rrsun)
-real, dimension(:), intent(in) :: lat
-type(time_type), intent(in) :: time
-real, dimension(:), intent(out) :: cosz, fracday
-real, intent(out) :: rrsun
+subroutine DAILY_MEAN_SOLAR_CAL_1D_(lat, time, cosz, fracday, rrsun)
+
+real(kind=FMS_AST_KIND_), dimension(:), intent(in) :: lat
+type(time_type), intent(in) :: time
+real(kind=FMS_AST_KIND_), dimension(:), intent(out) :: cosz, fracday
+real(kind=FMS_AST_KIND_), intent(out) :: rrsun
!---------------------------------------------------------------------
! local variables
!---------------------------------------------------------------------
- real, dimension(size(lat),1) :: lat_2d, cosz_2d, fracday_2d
+real(kind=FMS_AST_KIND_), dimension(size(lat),1) :: lat_2d, cosz_2d, fracday_2d
!--------------------------------------------------------------------
!> define 2-d versions of input data array.
!--------------------------------------------------------------------
- lat_2d(:,1) = lat
+
+ lat_2d(:,1) = lat
!--------------------------------------------------------------------
!> call daily_mean_solar_cal_2d to convert the time_types to reals and
!! then calculate the astronomy fields.
!--------------------------------------------------------------------
- call daily_mean_solar_cal_2d (lat_2d, time, cosz_2d, &
- fracday_2d, rrsun)
+
+ call daily_mean_solar(lat_2d, time, cosz_2d, fracday_2d, rrsun)
!-------------------------------------------------------------------
!> place output fields into 1-d arguments for return to
!! calling routine.
!-------------------------------------------------------------------
- fracday = fracday_2d(:,1)
- cosz = cosz_2d(:,1)
-end subroutine daily_mean_solar_cal_1d
+ fracday = fracday_2d(:,1)
+ cosz = cosz_2d(:,1)
+
+end subroutine DAILY_MEAN_SOLAR_CAL_1D_
!> @brief daily_mean_solar_cal_2level receives 1d arrays and time_type input,
@@ -1651,39 +1050,42 @@ end subroutine daily_mean_solar_cal_1d
!! @param [out]